Old Matlab example of 1D FFT filter
13 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Valeriy
el 20 de Abr. de 2024
Editada: Star Strider
el 28 de Abr. de 2024
I remember that in one of the old Matlab version (2010 or even earlier), in its Help was shown example of the application FFT - IFFT filter to remove noise frequency components of signal, which was close to the sinusoidal. Example was short and useful, but now I need something similar, and can't find it. Is it possible to recover it?
2 comentarios
Paul
el 27 de Abr. de 2024
Hi Valeriy,
Do you recall if that example was showing how to remove components of a periodic signal? Or was it showing how to apply an LTI filter to a finite duration signal? Or maybe something else (though nothing comes to mind as what that could be)?
Respuesta aceptada
Star Strider
el 20 de Abr. de 2024
The only function that I am aware of that might do what you want is the fftfilt function (introduced before R2006a). It requires a previously-designed FIR filter denominator ‘b’ vector, however this is straightforward in MATLAB.
10 comentarios
Star Strider
el 25 de Abr. de 2024
Editada: Star Strider
el 28 de Abr. de 2024
As always, my pleasure!
I do not remember actually seeing that in any documentation. I have posted answers on that in the past, so that may be where you saw the code. (As I mentioned, I do not recommend that sort of approach to filtering, however it has come up from time to time.)
An example could be something like this —
Fs = 500; % Sampling Frequency (Hz)
Lt = 5; % Signal Length (s)
t = linspace(0, Fs*Lt, Fs*Lt+1)/Fs; % Time Vector
N = 50; % Number Of Frequencies In Signal
A = rand(1,N); % Signal Component Amplitudes
freqs = randi(250, 1, N); % Signal Component Frequencies
s = sum(A(:).*sin(2*pi*t.*freqs(:))); % Create Signal
figure
plot(t, s)
grid
xlabel('Time')
ylabel('Amplitude')
title('Original Signal')
Ls = numel(t);
NFFT = 2^nextpow2(Ls);
wf = hann(Ls); % Window Function
FTs = fft(s(:).*wf, NFFT)/sum(wf); % Fourier Transform
FTss = fftshift(FTs);
Fvs = Fs*(-(NFFT/2) : (NFFT/2))/NFFT; % Frequency Vector
Fvs(ceil(numel(Fvs/2))) = [];
figure
plot(Fvs,abs(FTs))
xt = xticks;
xtl = [xt(fix(numel(xt)/2)+1:end) flip(xt(fix(numel(xt)/2)+1:end-1))];
xticklabels(xtl)
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Fourier Transform Of Original Signal')
figure
plot(Fvs, abs(FTss))
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Shifted Fourier Transform Of Original Signal')
LPF = zeros(size(Fvs)); % Lowpass Filter
LPF((Fvs >= -100) & (Fvs <= 100)) = 1; % Lowpass Filter
figure
plot(Fvs, LPF)
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Lowpass Filter')
axis('padded')
FTss_filt = FTss .* LPF(:); % Filter Signal
figure
plot(Fvs, abs(FTss_filt))
grid
xlabel('Frequency')
ylabel('Magnitude')
title('Shifted Fourier Transform Of Filtered Signal')
s_filt = ifft(ifftshift(FTss_filt), 'symmetric'); % Inverse Fourier Transform
s_filt = s_filt(1:numel(t)); % Eliminate Zero-Padding (Added At End Of Signal)
figure
plot(t, s_filt)
grid
xlabel('Time')
ylabel('Amplitude')
title('Filtered Time-Domain Signal')
The filter eliminates frequencies above 100 Hz, as well as eliminating a large amount of the signal energy in the original signal.
EDIT — (28 Apr 2024 at 03:44)
The original signal has the fft as two essentially symmetric halves, the second half being the complex conjugate of the first half. The fftshift operation creates the vector as the second half being the flipped mirror image of the first half, with the two halves flipped as well. That makes it easier to do the symmetric filtering, since it is easier to write the filter code. Once the filtering is accomplished, the ifftshift function returns the vectors to their original orientation, so the ifft function can produce the filtered version of the original signal.
I added an extra plot of the original fft result (with an appropriate frequency axis) and a ‘sort of’ Bode plot of the filter transfer function (with padded axes) to illustrate this, and the filtering operation.
The code is otherwise unchanged.
.
Más respuestas (0)
Ver también
Categorías
Más información sobre Get Started with Signal Processing Toolbox en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!