2014-12-09 1 views
0

Я пытаюсь найти обратное преобразование Фурье простого фильтра в Matlab. В первом случае (фильтр sinc/«кирпичная стена») я использую функцию , чтобы найти функцию во временной области, которая является центром с центром в t = 0.Обратный ФТ фильтра в Matlab

Теперь я хочу найти время -домена для простого чебышевского фильтра. Однако по какой-то причине нижеприведенный код дает мне импульсный отклик, где временная ось неверна. Должен ли я ожидать, что подобная нетерминальная функция с центром в t = 0?

fo = 10; %frequency of the sine wave 
Fs = 100; %sampling rate 
Ts = 1/Fs; %sampling time interval 
t = -1+Ts:Ts:1-Ts; %sampling period 
freq = -Fs/2:(Fs/length(t)):Fs/2; 

%% Sinc with bandwidth = fo. This works! 
y = 0.5*sinc(2*fo*t); 
YfreqDomain1 = fft(y); 
figure('Name','Brick wall sinc filter (freq)'); 
plot(freq(1:length(y)),2/length(y)*fftshift(abs(YfreqDomain1))) 
y_ret1=ifft(YfreqDomain1,'nonsymmetric'); 
figure('Name','Brick wall sinc filter (time)'); 
plot(t,y_ret1); 

%% Chebyshev with bandwidth fo. This gives me a strange result. 
[b,a] = cheby1(6,0.1,2*fo/Fs); % 6th order, 0.1dB ripple 
[YfreqDomain2 w] = freqz(b,a,length(t),'whole'); 
figure('Name','Chebyshev Filter (freq)'); 
plot(freq(1:length(YfreqDomain2)), 2/length(y)*fftshift(abs(YfreqDomain2))); 
figure('Name','Chebyshev Filter (time)'); 
y_ret2=ifft(YfreqDomain2,'nonsymmetric'); 
plot(t,y_ret2); 
+0

Я думаю, что импульсная характеристика является желаемым результатом. Отфильтрованный сигнал 'y (t)' является сверткой между входным сигналом 'x (t)' и импульсной характеристикой 'h (t)', поэтому 'h (t) = IDFT [H (w)]' является импульсивный ответ. Или я там не прав? – hbaderts

+1

Проблема, по-видимому, связана с различиями в фазе двух сигналов. 'plot (t, fftshift (ifft (abs (YfreqDomain2))) выглядит корректно и сопоставимо с той же командой с YfreqDomain1. –

+1

Не причина вашей проблемы, но также стоит отметить, что YfreqDomain1 равен 199x1, тогда как YfreqDomain2 равен 1x199 –

ответ

2

Есть две проблемы с вашим вычислением:

Во-первых, вы оценить временные коэффициенты фильтра домена на очень узкое окно времени, которое обрезает фильтр и изменяет свои характеристики.

Во-вторых, вы не правильно отслеживаете, какие индексы в векторах временной области и частотной области соответствуют временным точкам 0 и частоте 0 соответственно. Это можно сделать по-другому, я выбрал всегда t(1) = 0 и f(1) = 0, и примените fftshift только для печати.

Вот моя исправленная версия кода:

fo = 10;  % frequency of the sine wave 
Fs = 100;  % sampling rate 
Ts = 1/Fs; % sampling time interval 
n = 1000;  % number of samples 

% prepare sampling time vector such that t(1) = 0 
t = (0 : n - 1)'; 
t = t - n * (t >= 0.5 * n); 
t = t/Fs; 

% prepare frequency vector such that f(1) = 0; 
f = (0 : n - 1)'/n; 
f = f - (f >= 0.5); 
f = f * Fs; 


%% sinc filter with bandwidth fo 

% define filter in time domain 
s = 2*fo/Fs * sinc(2*fo*t); 

% transform into frequency domain 
Hs = fft(s); 
% since the filter is symmetric in time, it is purely real in the frequency 
% domain. remove numeric deviations from that: 
Hs = real(Hs); 

subplot(2, 4, 1) 
plot(fftshift(t), fftshift(s)) 
ylim([-0.1 0.25]) 
title('sinc, time domain') 

subplot(2, 4, 2) 
plot(fftshift(f), real(fftshift(Hs)), ... 
    fftshift(f), imag(fftshift(Hs))) 
ylim([-1.1 1.1]) 
title({'sinc, frequency domain', 'real/imaginary'}) 

subplot(2, 4, 3) 
plot(fftshift(f), abs(fftshift(Hs))) 
ylim([-0.1 1.1]) 
title({'sinc, frequency domain', 'modulus'}) 


%% Chebyshev filter with bandwidth fo 

% define filter in frequency domain 
[b,a] = cheby1(6, 0.1, 2*fo/Fs); % 6th order, 0.1 dB ripple 
Hc = freqz(b, a, n, 'whole', Fs); 

% transform into time domain 
c = ifft(Hc); 

% determine phase such that phase(1) is as close to 0 as possible 
phase = fftshift(unwrap(angle(fftshift(Hc)))); 
phase = phase - 2*pi * round(phase(1) /2/pi); 

subplot(2, 4, 5) 
plot(fftshift(t), fftshift(c)) 
title('Chebyshev, time domain') 
ylim([-0.1 0.25]) 

subplot(2, 4, 6) 
plot(fftshift(f), real(fftshift(Hc)), ... 
    fftshift(f), imag(fftshift(Hc))) 
ylim([-1.1 1.1]) 
title({'Chebyshev, frequency domain', 'real/imaginary'}) 

subplot(2, 4, 7) 
plot(fftshift(f), abs(fftshift(Hc))) 
ylim([-0.1 1.1]) 
title({'Chebyshev, frequency domain', 'modulus'}) 

subplot(2, 4, 8) 
plot(fftshift(f), fftshift(phase)) 
title({'Chebyshev, frequency domain', 'phase'}) 

И вот результат:

Как вы можете видеть, синк и Чебышева фильтры сходны по модулю частотного отклика, но очень различны в отношении фазы. Причина в том, что фильтр Чебышева равен causal filter, что означает, что коэффициенты во временной области ограничены 0 для t < 0, естественным свойством для фильтра, который реализован в реальной физической системе.

+1

Блестящий ответ - спасибо за помощь. – Graham

 Смежные вопросы

  • Нет связанных вопросов^_^