使用 fft 的 Matlab 低通滤波器

Question

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

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;

生成的图像如下所示

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

Answer 1

提醒自己注意 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')

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

高斯滤波器具有更好的低通滤波器特性，因为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')

如您所见，这种方式的重建效果要好得多。