Filter a signal to obtain a smooth fft

Hi,

I have a sequence of samples obtained in an accelerometer installed in a structure. If I do the fft of the signal, the natural frequency of the structure appears as a high peak in the fft graph.

I need to filter the accelerometer signal in order to obtain a smooth fft curve. If I use a bandstop filter to eliminate the natural frequency, all the frequency band dissapears and the curve has a deep valley.I want to eliminate only the peak to obtain a continuous distribution of frequencies.

Do you have any idea of how can I obtain this? Any advice will be great!

Lots of thanks for reading!

Answer 1

[Wayne King]

Maria, how about using a notch filter? I don't know what version of MATLAB you are running, but do you have the DSP System Toolbox, or the Filter Design Toolbox?

Fs = 1e3;
t = 0:1/Fs:1-(1/Fs);
x = cos(2*pi*100*t)+0.5*randn(size(t));
d  = fdesign.notch('N,F0,Q,Ap',6,100,10,1,1000);
Hd = design(d);
% View notch filter magnitude response
fvtool(Hd)
% filter signal
y = filter(Hd,x);
plot(psd(spectrum.periodogram,y,'Fs',Fs,'NFFT',length(y)))
figure;
plot(psd(spectrum.periodogram,x,'Fs',Fs,'NFFT',length(y)))

Wayne

[MariaLA]

The notch filter is like the bandstop filter, all a band of frequencies dissapear.

Imagine the deep valley in the Figure 2 not so deep, with the minimum in -35dB/Hz instead of in -62.I want to obtain something like that.

Maybe it is imposible with a filter, I don't know...

[Wayne King]

Hi Maria, yes, you can adjust the depth of the notch you obtain with the notch filter. Here is how you can get 30 dB of attenuation in the notch at 100 Hz.

d = fdesign.notch('N,F0,Q,Ast',6,100,30,30,1000);

Hd = design(d);

fvtool(Hd)

[MariaLA]

You are right.I've been trying the filter with different factors and the result seems valid.

Lots of thanks, Wayne.

Answer 2

[Wayne King]

Hi, You don't need to filter to produce a smooth fft(). Using a fft() directly on a time series and plotting the magnitude squared of the DFT coefficients is proportional to the periodogram. The periodogram is a low-degree-of-freedom estimate of the true power spectral density. The periodogram has great resolution, but suffers from high variance and inconsistency.

You can smooth the periodogram in a number of ways. You will lose some resolution but you may find the trade off acceptable.

You can obtain a smoother plot by using a Welch estimate or multitaper estimate

     Fs = 1e3;
     t = 0:.001:1-0.001;
     x = cos(2*pi*100*t)+randn(size(t));
     % compare
     plot(psd(spectrum.periodogram,x,'NFFT',length(x),'Fs',1e3));

     %with
     figure;
     hwelch = spectrum.welch;
     hwelch.SegmentLength = 100;
     plot(psd(hwelch,x,'NFFT',length(x),'Fs',1e3));

     % or
     figure;
     plot(psd(spectrum.mtm,x,'NFFT',length(x),'Fs',1e3));

Note that you do lose frequency resolution with both the Welch and multitaper estimators --- there is no free lunch.

[MariaLA]

Hi Wayne, thanks for the comment.

With the code of your post I can obtain a smoother spectrum of frequencies.

But what I really need is modifying my accelerometer signal to obtain another one with a fft() without high peaks in natural frequencies.

Usually if you need to eliminate a band of frequencies a "bandstop" filter is used. An example is shown below:

Fs = 1e3;

t = 0:.001:1-0.001;

x = cos(2*pi*100*t)+randn(size(t));

Fn=Fs/2;

%normalized band of frequencies to stop

Wn(1)=900/Fn;

Wn(2)=105/Fn;

n=7;

[b,a]=butter(n,Wn,'stop');

Y=filter(b,a,x);

% compare periodograms

subplot(2,1,1)

plot(psd(spectrum.periodogram,x,'NFFT',length(x),'Fs',1e3));

title('Original signal')

subplot(2,1,2)

plot(psd(spectrum.periodogram,Y,'NFFT',length(Y),'Fs',1e3));

title('Filtered signal')

I want to eliminate from x the peak at 100Hz using a 'bandstop' filter. The fft() of the result (y) is shown in 'Filtered signal'. In the 95-105Hz band appears a deep valley.

Neither the signal x with the peak nor the signal Y with the deep valley are valid. I want to obtain a signal with an almost continuous frequency distribution.

Do you know what I mean?

Thanks in advance.

[Wayne King]

Yes, I do.

Do you want to just eliminate a peak at 100 Hz, or really a band of frequencies around 100 Hz?

[MariaLA]

I want to eliminate only the peak of resonance frequency, but to configure the bandstop filter I had to define a band of frequencies around the peak.

A comment: In the code above I made a mistake

Wn(1)=90/Fn or Wn(1)=95/Fn for example, not 900