使用 fft 的 Matlab 低通滤波器

Matlab Low Pass filter using fft

我使用前向和后向 fft 在 matlab 中实现了一个简单的低通滤波器。 原理上是可以的,但是最小值和最大值和原来的不一样

signal = data;
%% fourier spectrum
% number of elements in fft
NFFT = 1024;
% fft of data
Y = fft(signal,NFFT)/L;
% plot(freq_spectrum)

%% apply filter
fullw = zeros(1, numel(Y));
fullw( 1 : 20 ) = 1;
filteredData = Y.*fullw;

%% invers fft
iY = ifft(filteredData,NFFT);
% amplitude is in abs part
fY = abs(iY);
% use only the length of the original data
fY = fY(1:numel(signal));
filteredSignal = fY * NFFT; % correct maximum

clf; hold on;
plot(signal, 'g-')
plot(filteredSignal ,'b-')
hold off;

生成的图像看起来像这样

我做错了什么?如果我对这两个数据进行归一化,过滤后的信号看起来是正确的。

只是为了提醒我们自己 MATLAB 如何存储 Y = fft(y,N):

的频率内容
  • Y(1)是常量偏移量
  • Y(2:N/2 + 1)是正频率的集合
  • Y(N/2 + 2:end) 是一组负频率...(通常我们会绘制这个 left 的垂直轴)

为了制作真正的低通滤波器,我们必须保留低正频率频率。

这是一个在频域中使用乘法矩形滤波器执行此操作的示例,就像您所做的那样:

% make our noisy function
t = linspace(1,10,1024);
x = -(t-5).^2  + 2;
y = awgn(x,0.5); 
Y = fft(y,1024);

r = 20; % range of frequencies we want to preserve

rectangle = zeros(size(Y));
rectangle(1:r+1) = 1;               % preserve low +ve frequencies
y_half = ifft(Y.*rectangle,1024);   % +ve low-pass filtered signal
rectangle(end-r+1:end) = 1;         % preserve low -ve frequencies
y_rect = ifft(Y.*rectangle,1024);   % full low-pass filtered signal

hold on;
plot(t,y,'g--'); plot(t,x,'k','LineWidth',2); plot(t,y_half,'b','LineWidth',2); plot(t,y_rect,'r','LineWidth',2);
legend('noisy signal','true signal','+ve low-pass','full low-pass','Location','southwest')

完整的低通滤波器做得更好,但您会注意到重建有点 "wavy"。这是因为在频域中与矩形函数相乘与 convolution with a sinc function in the time domain 相同。与 sinc 函数的卷积用其邻居的非常不均匀的加权平均值替换每个点,因此 "wave" 效果。

高斯滤波器具有更好的低通滤波器特性,因为 the fourier transform of a gaussian is a gaussian。高斯很好地衰减到零,因此它在卷积期间的加权平均值中不包括遥远的邻居。这是一个保留正负频率的高斯滤波器的例子:

gauss = zeros(size(Y));
sigma = 8;                           % just a guess for a range of ~20
gauss(1:r+1) = exp(-(1:r+1).^ 2 / (2 * sigma ^ 2));  % +ve frequencies
gauss(end-r+1:end) = fliplr(gauss(2:r+1));           % -ve frequencies
y_gauss = ifft(Y.*gauss,1024);

hold on;
plot(t,x,'k','LineWidth',2); plot(t,y_rect,'r','LineWidth',2); plot(t,y_gauss,'c','LineWidth',2);
legend('true signal','full low-pass','gaussian','Location','southwest')

如您所见,这样重建效果更好。