Old Matlab example of 1D FFT filter
13 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Valeriy
am 20 Apr. 2024
Bearbeitet: Star Strider
am 28 Apr. 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 Kommentare
Paul
am 27 Apr. 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)?
Akzeptierte Antwort
Star Strider
am 20 Apr. 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 Kommentare
Star Strider
am 25 Apr. 2024
Bearbeitet: Star Strider
am 28 Apr. 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.
.
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Get Started with Signal Processing Toolbox finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!