【问题标题】:Inverse FT of a filter in MatlabMatlab中滤波器的逆FT
【发布时间】:2014-12-09 14:00:46
【问题描述】:

我试图在 Matlab 中找到一个简单滤波器的傅里叶逆变换。在第一种情况下(sinc 滤波器/“砖墙”),我使用ifft 函数找到时域函数,它是一个 sinc,以 t = 0 为中心。

我现在想找到一个简单 Chebyshev 滤波器的时域函数。但是,由于某种原因,下面的代码似乎给了我脉冲响应,其中时间轴不正确。我不应该期待以 t = 0 为中心的类似外观的 sinc 函数吗?

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);

【问题讨论】:

  • 我认为脉冲响应是想要的结果。滤波后的信号y(t)是输入信号x(t)和脉冲响应h(t)的卷积,所以h(t) = IDFT[H(w)]就是脉冲响应。还是我错了?
  • 问题似乎与两个信号的相位差异有关。 plot(t,fftshift(ifft(abs(YfreqDomain2))) 看起来是正确的,并且与使用 YfreqDomain1 的同一命令相当。
  • 不是问题的原因,但也值得注意的是 YfreqDomain1 是 199x1 而 YfreqDomain2 是 1x199
  • 感谢吉姆 - 成功了。但是,为什么绘图不像第一种情况那样是对称的类 sinc 函数?它看起来像一个脉冲响应,奇怪的是,周期是第一种情况的两倍。 hbaderts,如果 x(t) 是狄拉克函数,情况就是这样。但是,在这种情况下,我们只是在做一个变换?
  • 为了确保我们看到的是相同的东西......我上面给出的 plot 命令(假设你在最后添加了缺少的))应该产生一个对称的类似 sinc 的函数。至于为什么您的原始代码显示时移版本,目前对我来说是个谜

标签: matlab filter fft ifft


【解决方案1】:

您的计算存在两个问题:

首先,您在一个非常窄的时间窗口上评估时域滤波器系数,这会截断滤波器并改变其特性。

其次,您没有正确跟踪时域和频域向量中的哪些索引分别对应于时间点 0 和频率 0。这可以通过不同的方式完成,我在这里选择始终使用t(1) = 0f(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'})

结果如下:

如您所见,sinc 和 Chebyshev 滤波器在频率响应的模数方面相似,但在相位方面却大不相同。原因是切比雪夫滤波器是causal filter,这意味着时域中的系数在 t

【讨论】:

    猜你喜欢
    • 2014-12-06
    • 1970-01-01
    • 2011-10-02
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2011-02-15
    • 1970-01-01
    相关资源
    最近更新 更多