[kostya-m 0]

Разбирая алгоритмы, предоставляемые Матлабом по теме 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

[kostya-m 0]

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

 

Первичная задача (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');

Изменено 3 апреля, 2011 пользователем Мусатов Константин

[chipovod 0]

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

[kostya-m 0]

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

[kostya-m 0]

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

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

[Mad_kvmg 0]

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

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

 

 

[qxov 0]

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

doc frpintf

[bahurin 0]

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

[bahurin 0]

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

[bahurin 0]

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

 

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

 

[kostya-m 0]

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

doc frpintf

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

[kostya-m 0]

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

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 0]

Пробовал воспользоваться приведенным здесь примером.

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

При подстановке фазы - увы нет.