окно Дольфа-Чебышёва

Все что касается фильтрации
Ivan Karamazov
Сообщения: 89
Зарегистрирован: 28 окт 2010, 22:31
Откуда: Москва

Re: окно Дольфа-Чебышёва

Сообщение Ivan Karamazov »

Требование одно -- чтобы наилучшим образом по тем или иным критериям (в конечном счете среднеквадратичным) приближать АЧХ КИХ фильтра к недостижимому идеалу. Плюс, разумеется, что именно с тем или иным окном "всего виднее" на спектроанализаторе.
Трудность собственно одна -- "magic numbers" очень смущают. Ну и, конечно, формальные описания окон -- но тут, кажется, не поваляешь - не поешь.
Если ваши решения вам нравятся -- это хорошие решения. И наоборот.

eiyawii
Сообщения: 60
Зарегистрирован: 20 окт 2011, 07:38

Re: окно Дольфа-Чебышёва

Сообщение eiyawii »

Ivan Karamazov писал(а):Трудность собственно одна -- "magic numbers" очень смущают.
что есть мэджик намберз? не понял суть утверждения.)
Ivan Karamazov писал(а):формальные описания окон -- но тут, кажется, не поваляешь - не поешь.
вы имеете в виду необходим практический подход для понимания?

Ivan Karamazov
Сообщения: 89
Зарегистрирован: 28 окт 2010, 22:31
Откуда: Москва

Re: окно Дольфа-Чебышёва

Сообщение Ivan Karamazov »

magic numbers (термин дословно непереводимый) -- это [в частности и примерно] коэффициенты в коде типа:

Код: Выделить всё

/* The Blackman–Harris window
*/
inline double winBlackmanHarris(int n, int M)
{
 double arg = FIRWindows_NDIV_;
 return 0.35875 -
	0.48829 * cos(2.0 * M_PI * arg) +
	0.14128 * cos(4.0 * M_PI * arg) -
	0.01168 * cos(6.0 * M_PI * arg);
}

/* The Nuttall window
*/
inline double winNuttall(int n, int M)
{
 double arg = FIRWindows_NDIV_;
 return 0.355768 -
	0.487396 * cos(2.0 * M_PI * arg) +
	0.144232 * cos(4.0 * M_PI * arg) -
	0.012604 * cos(6.0 * M_PI * arg);
}

/* The Blackman–Nuttall window
*/
inline double winBlackmanNuttall(int n, int M)
{
 double arg = FIRWindows_NDIV_;
 return 0.3635819 -
	0.4891775 * cos(2.0 * M_PI * arg) +
	0.1365995 * cos(4.0 * M_PI * arg) -
	0.0106411 * cos(6.0 * M_PI * arg);
}

/* The Flat top window
*/
inline double winFlatTop(int n, int M)
{
 double arg = FIRWindows_NDIV_;
 return 1.0 -
	1.93  * cos(2.0 * M_PI * arg) +
	1.29  * cos(4.0 * M_PI * arg) -
	0.388 * cos(6.0 * M_PI * arg) +
	0.032 * cos(8.0 * M_PI * arg);
}
(код потенциально для выкладывания, посему комментарии на ломаном английском.)
Ну а практический подход необходим. (пример -- на аудиокабели имеет свойство наводиться строчная развертка ТВ-вещания -- не то что бы сильно, но если специально слушать -- то можно услышать, во всяком случае, на фоне "цифровой тишины"). Рассмотреть таковую (сразу, а не на слух, каковой может подвести) (и принять меры) позволяло многие годы только окно Блэкмена(?)-Харриса -- хотя, "по классификации" -- оно "одно из многих". Т.е. расширяя репертуар не хочется опечаток. И не боле.
Т.е. на самом деле беда в том, что я за последний год вполне конкретно разучился верить "печатному слову". Многое чего, конечно, можно перевывести ручкой и бумажкой, но здесь немного не тот случай.
Если ваши решения вам нравятся -- это хорошие решения. И наоборот.

eiyawii
Сообщения: 60
Зарегистрирован: 20 окт 2011, 07:38

Re: окно Дольфа-Чебышёва

Сообщение eiyawii »

нашел по данной теме интересную статью: полностью описывает как провести синтез окна Дольфа-Чебышева.

PETER LYNCH, The Dolph–Chebyshev Window: A Simple Optimal Filter.

В статье формула (3) является расчетной формулой для коэффициентов окна (или я не правильно понял?). Сам сейчас рассчитываю их реккурентным способом и иногда кажется, что ошибка округления накапливается слишком большая. Эта формула могла бы пролить свет на точность вычислений моим способом, но как вычислить cosh ? при вычислениях возникает ошибка "overflow encountered in cosh"

Код: Выделить всё

wn = [] 								# empty array for coefficients
M = 25 									# just a number
N = 2*M + 1 							# number of coefficients
r = 10**(-3) 							# sidelobe attenuation
x0 = cosh(1./2*M)*pow(cosh(1./r), -1) 	# stop band frequency
for n in range(-M, M+1):
	S = 0; theta_n = 2*pi*n*N
	for m in range(1, M+1):
		theta_m = 2*pi*m/N
		S += cheb(2*M, x0*cos(theta_m/2.))*cos(m*theta_n)
	wn.append((1./N)*(1+2*r*S))
PS ссылка на статью: http://www.maths.tcd.ie/~plynch/Publications/Dolph.pdf

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: окно Дольфа-Чебышёва

Сообщение Бахурин Сергей »

Думаю ошибка в формуле

Код: Выделить всё

x0 = cosh(1./2*M)*pow(cosh(1./r), -1)    # stop band frequency
поскольку в статье в формуле 4 обозначение означает не , а функцию обратную гиперболическому косинусу (аркосинус гиперболический acosh)

eiyawii
Сообщения: 60
Зарегистрирован: 20 окт 2011, 07:38

Re: окно Дольфа-Чебышёва

Сообщение eiyawii »

исправил

Код: Выделить всё

x0 = cosh((1./2*M)*acosh(1./r)) 
но теперь возникает ошибка в этой строке

Код: Выделить всё

S += abs(cos(2*M*arccos(x0*cos(theta_m/2.))))*cos(m*theta_n)
возможно из-за того, что х0 является большим числом: 1.27999744e+26 (это при r=0.001, а если меньше брать r, то х0 увеличивается крайне быстро)

Аватара пользователя
Бахурин Сергей
Администратор
Сообщения: 1114
Зарегистрирован: 05 окт 2010, 19:55
Контактная информация:

Re: окно Дольфа-Чебышёва

Сообщение Бахурин Сергей »

Код: Выделить всё

S += abs(cos(2*M*arccos(x0*cos(theta_m/2.))))*cos(m*theta_n)
не понял зачем точка после 2: theta_m/2. может так поможет

Код: Выделить всё

S = S+abs(cos(2*M*arccos(x0*cos(theta_m/2))))*cos(m*theta_n)
насколько помню матлаб не допускает операторов типа +=

eiyawii
Сообщения: 60
Зарегистрирован: 20 окт 2011, 07:38

Re: окно Дольфа-Чебышёва

Сообщение eiyawii »

Бахурин Сергей писал(а):

Код: Выделить всё

S += abs(cos(2*M*arccos(x0*cos(theta_m/2.))))*cos(m*theta_n)
не понял зачем точка после 2: theta_m/2. может так поможет

Код: Выделить всё

S = S+abs(cos(2*M*arccos(x0*cos(theta_m/2))))*cos(m*theta_n)
насколько помню матлаб не допускает операторов типа +=
1. просто 2 означает целый тип, 2. - означает тип с плавающей точкой

2. да, матлаб += не переваривает, но код написал на python :)

кажется что в х0 ошибка, хотя я как ни всматривался - ошибку не нашел:

Код: Выделить всё

x0 = cosh((1./2*M)*acosh(1./r))

Ответить