使用fft的Matlab低通滤波器 [英] Matlab Low Pass filter using fft
问题描述
我在Matlab中使用前向和后向fft实现了一个简单的低通滤波器. 原则上可以使用,但是最小值和最大值与原始值不同.
I implemented a simple low pass filter in matlab using a forward and backward fft. It works in principle, but the minimum and maximum values differ from the original.
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;
生成的图像如下所示
我做错了什么?如果我将两个数据都归一化,则滤波后的信号看起来正确.
What am I doing wrong? If I normalize both data the filtered signal looks correct.
推荐答案
只是提醒我们自己MATLAB如何存储Y = fft(y,N)
的频率内容:
Just to remind ourselves of how MATLAB stores frequency content for Y = fft(y,N)
:
-
Y(1)
是常量偏移量 -
Y(2:N/2 + 1)
是一组正频率 -
Y(N/2 + 2:end)
是一组负频率...(通常我们会绘制垂直轴的这个 left )
Y(1)
is the constant offsetY(2:N/2 + 1)
is the set of positive frequenciesY(N/2 + 2:end)
is the set of negative frequencies... (normally we would plot this left of the vertical axis)
为了制作一个真正的低通滤波器,我们必须同时保留低的正频率和低的负频率.
In order to make a true low pass filter, we must preserve both the low positive frequencies and the low negative frequencies.
下面是一个使用频域中的乘法矩形滤波器执行此操作的示例,就像您所做的那样:
Here's an example of doing this with a multiplicative rectangle filter in the frequency domain, as you've done:
% 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')
完整的低通拟合器做得更好,但是您会注意到重建有点波浪".这是因为在频域中与矩形函数的乘法与在当时与Sinc函数的卷积相同域.具有正弦函数的卷积会用邻居的加权平均值非常不均匀的方式替换每个点,从而产生波动"效应.
The full low-pass fitler does a better job but you'll notice that the reconstruction is a bit "wavy". This is because multiplication with a rectangle function in the frequency domain is the same as a convolution with a sinc function in the time domain. Convolution with a sinc fucntion replaces every point with a very uneven weighted average of its neighbours, hence the "wave" effect.
高斯滤波器具有更好的低通滤波器属性,因为高斯的傅立叶变换是高斯.高斯很好地衰减到零,因此在卷积期间加权平均值中不包括遥远的邻居.这是一个使用高斯滤波器保留正负频率的示例:
A gaussian filter has nicer low-pass filter properties because the fourier transform of a gaussian is a gaussian. A gaussian decays to zero nicely so it doesn't include far-off neighbours in the weighted average during convolution. Here is an example with a gaussian filter preserving the positive and negative frequencies:
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')
如您所见,这种方式的重建要好得多.
As you can see, the reconstruction is much better this way.
这篇关于使用fft的Matlab低通滤波器的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!