Remove specific frequencies from FFT signal and reconstruct the signal after filtering those frequencies

Hi,

I have a signal that shows a very distinctive peaks in the FFT.

Those high amplitudes are the 'noise' of the signal. I would like to remove that values from the original signal and to plot the filtered signal.

Fs = 4500; % Sampling frequency (fps)

T = 1/Fs; % Sampling period (s)

L = 900; % Length of signal (how many frames)

tt = (0:L-1)*T; % Time vector

thickness = detrend(thickness);

Y = fft(thickness);

P2 = abs(Y/L);

P1 = P2(1:L/2+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

figure(1111)

h1=plot(f(1:end),P1(1:end)) ;

title('Amplitude Spectrum')

xlabel('f [Hz]')

ylabel('Power [mm]')

ylim auto

[B,IX] = sort(P1); %order the amplitudes

A1=B(end); %amplitude of the first peak

A2=B(end-1); %amplitude of second peak

f1=f(IX(end)); %frequency of first peak

f2=f(IX(end-1)); %frequency of second peak

AmpTab=[A1 A2];

FreTab=[f1 f2];

Answer 1

[Bjorn Gustavsson]

Those high amplitudes are a noise in the signal. Not the noise in the signal. To remove such interference-components you will be better off doing it this way (remember the Fourier-transform of a real signal has symmetric real components and an anti-symmetric imaginary components, and also that the DC-component is real.):

fD = fft(data_t_tichkness(:,2)); % Discrete Fourier-transform of your data
subplot(2,1,1)
plot(abs(fD)) % Plot of its absolute values
hold on
[safD,idx] = sort(abs(fD),'descend'); % Sort in descending order, this makes indexing simpler
plot(idx(2:5),abs(fD(idx(2:5))),'r.') % DC-component will be the first, then 
                                      % the positive and negative components will 
                                      % have equal magnitudes and appear consecutively in idx
fD(idx(2:5)) = 0;  % Set the Fourier-components of the first and second spike to zero.
plot(abs(fD))      % Yup, they're gone.
subplot(2,1,2)
ifD = ifft(fD);    % inverse-Fourier-transform
plot(data_t_tichkness(:,2))          
hold on
plot(ifD)

HTH

[paloma paleo]

That's brilliant and it does exactly what I was trying to do.

However, I was expecting that the 'clean' signal would give me a different mean value than the orignal signal.

mean(data_t_tickness(:,2))=0.7051

mean(ifD)=0.7051

I guess this is different issue than matlab issue?

Another concern is that when calculating single side FFT, the high peaks were at 135 and 270. I could relate those values to 1xrpm and 2xrpm (i.e. 8100rpm/60=135 etc).

When calculating double side FFT, I am not sure how to relate those peaks to the rpm?

Thanks in advance for your support. I really appreciate it.

[Bjorn Gustavsson]

The two Fourier-components that is removed will be of the form: which are periodic over your intervall - and with the average value of zero.

The frequencies will be the same, except you will also have the negative frequencies.

[paloma paleo]

Thanks indeed

[paloma paleo]

Dear Bjorn,

I would like to do exactly the opposite. I want to leave the low frequencies with the high peaks and to remove the high frequencies. I estimated that the high frequencies are any frequency higher than 216Hz. Therefore, in the fft graph, it will be between 86 and 814. I wrote the following code, but when I reconstruct the signal it doesnt make sense.

fD = fft(thickness); % Discrete Fourier-transform of your data

subplot(2,1,1)

absolute=abs(fD);

plot(abs(fD)) % Plot of its absolute values

hold on

%[safD,idx] = sort(abs(fD),'descend'); % Sort in descending order, this makes indexing simpler

%plot(idx(2:5),abs(fD(idx(2:5))),'r.') % DC-component will be the first, then

% the positive and negative components will

% have equal magnitudes and appear consecutively in idx

%fD(idx(2:5)) = 0; % Set the Fourier-components of the first and second spike to zero.

fD(86:814) = 0;

plot(abs(fD)) % Yup, they're gone.

subplot(2,1,2)

ifD = ifft(fD); % inverse-Fourier-transform

plot(thickness,'r')

hold on

plot(ifD,'b')

I would appreciate it if you could explain me what is going on.

I wonder if I could use in instead the following command:

fs=4500;

fpass=216;

yl = lowpass(thickness,fpass,fs);

Thanks

[Bjorn Gustavsson]

When it comes to "semi-manually" setting the fft-coefficient to zero one has to keep in mind that the DC-component is the first in the array, then the first non-zero frequency-fourier-coefficient is the second and the last component in the array. When I do these tasks I typically get the indexing wrong in that I dont index matching pairs that have the same frequency, but miss with one or the other. This I embarrassingly still have to count up/down by listing: 2 and end, 3 and end-1, 4 and end-2, etc. Doing that you should be able to figure out if you are zeroing the conjugate pairs or if you're off by one. That typically shows up as strange results where one of the +f_i or -f_i fourier-coefficient pair are set to zero and the other not - this makes for confusing signals...

I come from image-processing, and when I need that type of low-pass filtering I typically use a simple call to conv2:

I_filtered = conv2(I_raw,conv2(ones(3)/3^3,conv2(ones(2)/2^2,ones(2)/2^2,'full'),'full'),'same');

But you can use all the more clever filter-window-utilities of matlab.

HTH

[paloma paleo]

I can see the issue. I am not too queen to use a 'semi-manually' method since I have a huge number of tests and it would be tricky.

I am not sure what your low-pass filtering does. The type of filter that I need, it is a very simple. I just want to remove the values of the signal from certain frequency. In my case, the frequency is factor of the rotational speed, so I can use e.g., 5x rpm to define the threshold of the frequency. All values with higher frequency than the threshold will be consider a noise. So I would like a filter that only takes in considaration the frequency value.

I have been reading about the type of filters and how to define them and I believe I need something really simple. The type of filters that I found are too fancy for what I need. I work also with image procesing/anaylisis.

So I thought that something like

yl = lowpass(I_raw,fpass,fs);

But I am not 100% sure if this is the correct filter.

Thanks again for you kind responses.

[Bjorn Gustavsson]

Well the lowpass filter should do what you want. If you're unsure about the filter characteristics you can reasonably easy simply check the documentation (it contains a lot of juicy information!), test what frequency response you get from a single delta-spike as input etc.