Форум разработчиков электроники ELECTRONIX.ru

Мусатов Константин

Apr 2 2011, 12:49

Разбирая алгоритмы, предоставляемые Матлабом по теме FIR, не смог найти нужную функцию.
Есть некий фильтр, заданный по контрольным точкам с амплитудой и фазой.

Код

Fo=192000/128;
F=[0 ,10 ,15 , 20,30 ,40 ,50 ,60 ,80 ,100,150,200,250,300,350,Fo/2];
A=[-200,-100,-20, -3,0 ,0 ,0 , 0, 0 , 0 , -6,-100,-130,-160,-190,-200];
P=[0,10 ,15, 20, 40,80 ,130,160,200,240,310,330,340,350,355,360];

привел его к нормальной комплексной форме.

Код

RealA=times(power(10,A/20),cosd(P));
ImgA=times(power(10,A/20),sind(P));
H=complex(RealA,ImgA);
Fn=times(F,2/Fo);

Нужно построить оптимальный FIR фильтр. Нашел функцию fir2, она согласилась принять задачу

Код

B=fir2(1024,Fn,H);
fvtool(B,1);

Однако, в результате, фаза не такая, как я просил, а линейная в полосе пропускания. Есть ли способ построить фильтр и с заданной ФЧХ?
Думал уже пойти напрямую, через FFT, но не смог пока сгладить АФЧХ и отсемплировать ее, перед тем как провести FFT

Мусатов Константин

Apr 3 2011, 23:28

Если кому интересно, решение получено.

Первичная задача (F в Герцах, А в дБ, Р в градусах):

Код

Fo=192000/256;
N=1024;
F=[0   ,10  ,15 , 20,30 ,40 ,50 ,60 ,80 ,100,125, 150,200,250,300,350,370,Fo/2];
A=[-200,-100,-20, -3,0  ,0  ,0  ,  0, 0 , 0 , -1, -6,-100,-130,-160,-190,-200,-200];
P=[0   ,10  ,15 , 20, 40,80 ,130,160,200,240,270,310,330 ,340 ,350 ,355 ,359,360];

Строю оптимальный фильтр по амплитуде

Код

An=power(10,A/20);
Fn=times(F,2/Fo);
Bc=fir2(1024,Fn,An)

;
Строю фазовую характеристику Ro

Код

Fi=freqspace(N+1,'whole');
Fi=times(Fi-1,Fo);
F(1)=1e-12;
Pi=interp1(horzcat(flipdim(times(F,-1),1),F),horzcat(flipdim(times(P,-1),1),P),Fi,'cubic');
Ro=complex(cosd(Pi),sind(Pi));

Применяю ее к фильтру

Код

Hi=Bc*dftmtx(length(Bc));
Hi=times(Hi,Ro);
Bi=real(Hi*conj(dftmtx(length(Hi)))/length(Hi));

Смотрю результат и проверяю полученный фильтр на тестовом сигнале

Код

fvtool(Bi,1);
y=wavread('d:\test.wav','double');
y=filter(Bi,1,y);
wavwrite(y,48000,24,'d:\resp.wav');

Чиповод

Apr 13 2011, 08:41

Еще в matlab есть функции fdesign.arbmag и fdesign.arbmagnphase, я вам в вашей соседней теме посоветовал. Если вдруг придется использовать - отпишите плиз как они у вас работают. Именно в железе. У меня руки не дошили попробовать, уж больно все красиво получается.

Мусатов Константин

Apr 13 2011, 21:36

У меня функция свободного построения фильтра после попытки генерации уходит в себя навсегда... Не задалось.

Мусатов Константин

Jun 19 2011, 10:26

Возник еще один вопрос.
Если я создаю Fir фильтр в fdatool, то там же можно коэффициенты выгрузить в файл. А как выгрузить коэффициенты для фильтра, синтезированного не в этой утилите. Вот в примере выше был получен фильтр Bi. Его можно рассмотреть в fvtool, даже можно увидеть там коэффициенты, но вывести в файл он не дает. Наверно есть некая процедурка, но я не нашел.

Mad_max

Jun 24 2011, 07:51

Цитата(Мусатов Константин @ Jun 19 2011, 14:26)

Возник еще один вопрос.

Как-то удалось красиво решить вопрос?

qxov

Jun 24 2011, 08:53

Цитата(Mad_max @ Jun 24 2011, 11:51)

Как-то удалось красиво решить вопрос?

doc frpintf

bahurin

Jun 24 2011, 10:52

если задана АЧХ и даже задана ФЧХ, то расчет ких фильтра сводится к взятию ifft от вашего комплексного к-та передачи. Если задана просто ачх, то можно добавить линейную фчх и снова ifft. Это называется метод расчета ких фильтра на основе частотной выборки. Подробно рассмотрен здесь

bahurin

Jun 24 2011, 10:52

если задана АЧХ и даже задана ФЧХ, то расчет ких фильтра сводится к взятию ifft от вашего комплексного к-та передачи. Если задана просто ачх, то можно добавить линейную фчх и снова ifft. Это называется метод расчета ких фильтра на основе частотной выборки. Подробно рассмотрен здесь

bahurin

Jun 24 2011, 10:54

если задана АЧХ и даже задана ФЧХ, то расчет ких фильтра сводится к взятию ifft от вашего комплексного к-та передачи. Если задана просто ачх, то можно добавить линейную фчх и снова ifft. Это называется метод расчета ких фильтра на основе частотной выборки. Подробно рассмотрен здесь

Просьба удалить повторяющиеся сообщения. Сорри.

Мусатов Константин

Jul 2 2011, 18:48

Отвалился на некоторое время от темы, потому не реагировал.

Цитата

doc frpintf

Спасибо, помогло

Мусатов Константин

Jul 2 2011, 21:51

Вот итоговый вариант с сохранением коэффициентов

Код

Fs=192000/256; % Частота и дециматор

N=512; % Длина
% Частоты
F=[1e-12   ,10   , 20,30 ,40 ,50 ,60 ,80 ,100,125, 150,160,Fs/2-1];
% Амплитуды в дБ
A=[-200,-200, 0,0  ,0  ,0  ,  0, 0 , 0 , -1, -6,-200,-200];
% Фазы в градусах с приведением на ВЧ к 0 (нормировка задержки)
P=[0   ,5  , 90, 180,39 ,13,0,-10,-150,-27,-10,0, 0];
%P=[0,0,0,0,0,0,0,0,0,0,0,0,0];

P = P * pi / 180  - (F * 2 * pi ); % перевод в радианы и добавка линии
An=power(10,A/20); % перевод в абсолютные значения усиления

% Таблица частот
Fi=freqspace(N,'whole');
Fi=times(Fi,Fs/2);
% Удвоенные частоты для вещественной функции
Ft=horzcat(F, flipdim(-F + Fs ,2));
% Меняем знак фазы
Pt=horzcat(P, flipdim(-P,2));
% Амплитуду не трогаем
At=horzcat(An, flipdim(An,2));
% Интерполяция амплитуды и фазы
Pi=interp1(Ft,Pt,Fi,'cubic');
Ai=interp1(Ft,At,Fi,'cubic');
% Комплексный коэффициент передачи
Ro=complex(times(Ai,cos(Pi)),times(Ai,sin(Pi)));
% Обратное преобразование Фурье
Ti=Ro*conj(dftmtx(length(Ro)))/length(Ro);
% Берем реальную часть фильра
Bi=real(Ti);
% Оконная функция
%Ok = rot90(kaiser(length(Ro),0.5));
%Bi=times(Ok,Bi);
fvtool(Bi,1);
%y=wavread('d:\test.wav','double');
%y=filter(Bi,1,y);
%wavwrite(y,48000,24,'d:\resp.wav');
fid = fopen( 'D:\WORK\Matlab\Prj1\Sub.dat', 'w');
fprintf( fid, '%e,\r\n', Bi );
fclose(fid);

mml

Dec 16 2012, 11:08

Пробовал воспользоваться приведенным здесь примером.
При подстановке желаемых координат амплитуды - все работает.
При подстановке фазы - увы нет.