Old Matlab example of 1D FFT filter

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?

[Paul]

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)?

[Valeriy]

Hi Paul, thanks for your message. As I remember, it was function of fft filtering of sinusoidal signal with some noise. Noise was removed by multiplication of the fft spectrum by set of ones and zeros, which worked as kind of filter, and filtered signal was obtained by ifft.

What I don't understand from example of Star Strider, is application, action of fftshift anad ifftshift.

Answer 1

[Star Strider]

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.

[Valeriy]

Thank you, Star Strider for your answer. No, it was simple fft transformation and to result of it was applied some function, which I forget, how it was expressed. This function was not obvious, at least for me, and it was key point in filtering of the spectrum, which then was processed by ifft and gives desired result. It was very simple and elegant example, I don't know why it was exculded from later releases.

[Star Strider]

My pleasure!

That is actually straightforward. It requires that the fftshift function be used on the original fft result, and then the filter needs to be designed to be similarly symmetric. (This is easier with a lowpass filter than with a bandpass or highpass filter.) That is element-wise multipliied by the shifted fft result, then that is back-shifted using ifftshift, and that result inverted using the ifft funciton. The ‘filters’ here are simple blocks of ones and zeros vectors that together equal the length of the shifted fft, although they can be more sophisticated. It likely disappeared from the documentation because that sort of filtering has problems, specifically because of the abrupt transitions in the filter passband edges.

For the record, I do not recommend this sort of approach to filtering (that being the reason I only provide a description and not an example), although others may have different opinions.

[Valeriy]

Thanks again, it seems now we are talking about the same function. I remember that I used several times this example for some noise-filtering properties of quasi-sinusoidal signals with intensive main harmonic, and it works correctly, at least for my cases. In fact, now I have something similar task to correct shape of such signal and suppress its high frequency components.

[Star Strider]

My pleasure!

There are much better filtering optionos available if you have the Signal Ptocessing Toolbox.

[Valeriy]

No, I don't have Signal Processing Toolbox and I'm not sure that it worth to buy it only for one small task, my manager will not understand me... What do you think, this function which you described, its application will produce the phase shift of the filtered result relatively original one?

And returning to my initial request, is it the way to recover the code of this function? With all of its drawbacks, it was present many years in corresponding Help section... I'd like to test it and to carefully analyze consequencies of its applications.

[Star Strider]

You will have to request that from MathWorks, using the Contact Support page. (Include the URL of this thread so you will not need to repeat everything here.) I am reasonably certain that it is not available in the Other Releases documentation page (that only goes back to R2019a).

[Valeriy]

Excellent idea, thanks a lot! I will do immediately...

> It likely disappeared from the documentation because that sort of filtering has problems, specifically because of the abrupt transitions in the filter passband edges.

If in the border area ones and zeros will be exchanged by something like 0.8, 0.6, 0.4, 0.2? Some kind of apodization?

[Star Strider]

As always, my pleasure!

That could work. I suggest that you consult the documentation for fftfilt (that I linked to) for an efficient way to do this. It may be necessary to transform a FIR filter to the frequency domain, or design appropriate transition regions in your filter to make this work correctly and to avoid the Gibbs phenomenon (also see Ringing (signal processing). (This is also the reason to avoid JPEG image compression. Always use PNG instead.)

[Valeriy]

As you recommended, I have requested help from Contact Support page. I have received numerous letters with many copies of Help pages, related to fft and ifft commands, but, unfortunately I didn't found between them code of the function we have discussed. Do you know, what is Matlab version, where it was shown and how this function/filter was named? This information will be very useful because my correspondent from Contact Support has no idea, how to answer these questions.

Thank you very much for our discussion and help

[Star Strider]

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.

.