| |
|
  |
 КИХ и Матлаб, Знакомлюсь со средством |
|
|
|
|
 Apr 2 2011, 12:49
|
Частый гость
Группа: Участник
Сообщений: 188
Регистрация: 10-10-06
Пользователь №: 21 172
|
Разбирая алгоритмы, предоставляемые Матлабом по теме 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
|
Частый гость
Группа: Участник
Сообщений: 188
Регистрация: 10-10-06
Пользователь №: 21 172
|
Если кому интересно, решение получено. Первичная задача (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 3 2011, 23:30
|
|
|
|
|
|
|
|
Jun 24 2011, 07:51
|
Местный
Группа: Свой
Сообщений: 377
Регистрация: 23-12-06
Из: Зеленоград
Пользователь №: 23 811
|
Цитата(Мусатов Константин @ Jun 19 2011, 14:26)  Возник еще один вопрос. Как-то удалось красиво решить вопрос?
|
|
|
|
|
|
|
|
Jun 24 2011, 08:53
|
Частый гость
Группа: Свой
Сообщений: 86
Регистрация: 22-03-07
Из: Санкт-Петербург
Пользователь №: 26 406
|
Цитата(Mad_max @ Jun 24 2011, 11:51) Как-то удалось красиво решить вопрос? doc frpintf
|
|
|
|
|
|
|
|
Jul 2 2011, 18:48
|
Частый гость
Группа: Участник
Сообщений: 188
Регистрация: 10-10-06
Пользователь №: 21 172
|
Отвалился на некоторое время от темы, потому не реагировал. Цитата doc frpintf Спасибо, помогло 
|
|
|
|
|
|
|
|
Jul 2 2011, 21:51
|
Частый гость
Группа: Участник
Сообщений: 188
Регистрация: 10-10-06
Пользователь №: 21 172
|
Вот итоговый вариант с сохранением коэффициентов Код 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);
|
|
|
|
|
|
|
|
Dec 16 2012, 11:08
|
Участник
Группа: Участник
Сообщений: 43
Регистрация: 10-11-10
Из: Екатеринбург
Пользователь №: 60 777
|
Пробовал воспользоваться приведенным здесь примером.
При подстановке желаемых координат амплитуды - все работает.
При подстановке фазы - увы нет.
|
|
|
|
|
|
|
|
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|
