Why do I get too high dominant frequencies in FFT

Hi there,

I am measuring a signal with a vibrometer and am getting the following raw result from data very similar to that of the attachment (this was taken from a different measurement but is basically the same):

However, when I try to do the FFT using the example from here: https://se.mathworks.com/help/matlab/ref/fft.html , I get mostly large frequency spikes in the hundreds or thousands of Hz range. I also do not get the main dominant frequency from the signal above, i.e. 18 Hz when counting the number of peaks in one second. Why are these different and why am I not getting the main dominant frequency in the signal? The frequency spikes I get in the FFT can be seen below:

I know the scale is a bit bad but none of those signals are around 18 Hz, and the dominant one is about 2000 Hz. The code I am using can also be seen below:

Fs = s.Rate;          % Sampling frequency in Hz
T = 1/Fs;             % Sampling period       
L = (s.DurationInSeconds)*1000;          % Length of signal in ms
t = (0:L-1)*T;        % Time vector

tic
s.DurationInSeconds = 2;

[data, time] = s.startForeground;
toc
figure
plot(time,data)
xlabel('Time (s)')
ylabel('Voltage (V)')

Y = fft(data);

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 
plot(f,P1) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

Any help is greatly appreciated. Cheers

[dpb]

What is the sample rate?

W/o the data not much anybody here can do, but while there are the large spikes with the 18 Hz frequency (presuming the time scale is correct), that they are individually of large magnitude they appear to be of very short duration so they don't necessarily have very much of the total energy of the signal which appears to be in the very high frequency stuff in between those.

Couple things to try while we wait for more info -- if the high frequency stuff isn't real and of interest, use a lowpass filter to get rid of a lot of it -- if you could use an analog anti-aliasing filter on the data collection would be best.

It's possible what you're seeing is aliased stuff being folded back if the sampling rate isn't high enough.

Not having any idea what the signal is, no ideas on what might be things to try to clean it up otherwise.

[Timo Mayer]

Hi,

thanks for the reply. The sampling rate is 48000 Hz (sorry I thought it was in the code but just realised it wasn't). Is that high enough?

The data is from a laser vibrometer measuring vibrations from impact hammering on a piece of material. A 5 kHz low-pass filter in the controller is already being used. What is an example of a filter I can implement in the code for Matlab?

Let me know if you need more info.

Cheers.

[David Goodmanson]

Hi Timo,

something does not appear to be consistent here. The t array is not used. Could you supply the answers for size(time), and time(end)?

[Timo Mayer]

Hi David,

the time is 2 seconds so we start at 0 s and end at 2 s. What t array are you referring to?

Cheers,

Timo

[David Goodmanson]

Hi Timo,

the data you are actually using for the fft comes from

[data, time] = s.startForeground;

so I was interested in the size of time. I see the range is shown in the first plot, but there is no way to know the number of points. Also I am speculating that L = 2000 but but it's not so clear because s.DurationInSeconds is defined as equal to 2 only after the calculation of L using the original s.DurationInSeconds is done.

Answer 1

[Mathieu NOE]

hello

so your data file has two channels of data and no time data (a contrario from what is written in the header)

below some example of code for fft averaging (pwelch) and sepcgram

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 4096*4;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting : dB(L))
% option_w = 1 : A weighted spectrum  : (dB (A) ) for acoustic measurements
option_w = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


[data,head] = readclm('Data.txt',2,1);
signal = data(:,1);

samples = length(signal);

Fs = 48000 ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(1),semilogx(freq,sensor_spectrum_dB);grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sg,fsg,tsg] = specgram(signal,NFFT,Fs,w,NOVERLAP);  

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end

% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');


function pondA_dB = pondA_function(f)
% dB (A) weighting curve
n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
pondA = n./r;
pondA_dB = 20*log10(pondA(:));
end

% %%%%%%%%%%% other infos %%%%%%%%%%%%%%%%%%%%%%%%
% 
% 
% %% notch filter section %%%%%%
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % notch :      H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fc = 900;       % notch freq
% 
% Q = 10;         % the higher the Q factor, the deeper the notch
% 
% %%%%%%
% w0 = 2*pi*fc;
% 
% alpha = sin(w0)/(2*Q);
% 
% 
%             b0 =   1;
%             b1 =  -2*cos(w0);
%             b2 =   1;
%             a0 =   1 + alpha;
%             a1 =  -2*cos(w0);
%             a2 =   1 - alpha;
% % 
% num1z=[b0 b1 b2];
% den1z=[a0 a1 a2];
% 
% % now let's filter the signal 
% signal_filtered = filter(num1z,den1z,signal);

in the attachement some bonus on digital filters you can use (in conjunction with filter command)

function  [outdata,head] = readclm(filename,nclm,skip,formt)

% READCLM Reads numerical data from a text file into a matrix.
%	Text file can begin with a header or comment block.
%	[DATA,HEAD] = READCLM(FILENAME,NCLM,SKIP,FORMAT)
%	Opens file FILENAME, skips first several lines specified
%	by SKIP number or beginning with comment '%'.
%	Then reads next several lines into a string matrix HEAD
%	until the first line with numerical data is encountered
%	(that is until first non-empty output of SSCANF).
%	Then reads the rest of the file into a numerical matrix
%	DATA in a format FORMAT with number of columns equal
%	to number of columns of the text file or specified by
%	number NCLM. If data does not match the size of the
%	matrix DATA, it is padded with NaN at the end.
%
%	READCLM(FILENAME) reads data from a text file FILENAME,
%	skipping only commented lines. It determines number of
%	columns by the length of the first data line and uses
%	the floating point format '%g';
%
%	READCLM uses FGETS to read the first lines and 	FSCANF
%	for reading data.

%  Kirill K. Pankratov,  [e-mail]
%  03/12/94, 01/10/95.

 % Defaults and  parameters ..............................
formt_dflt = '%g';  % Default format for fscanf
addn = nan;         % Number to fill the end if necessary

 % Handle input ..........................................
if nargin<1, error('  File name is undefined'); end
if nargin<4, formt = formt_dflt; end
if nargin<3, skip = 0; end
if nargin<2, nclm = 0; end
if isempty(nclm), nclm = 0; end
if isempty(skip), skip = 0; end

 % Open file ............................
[fid,msg] = fopen(filename);
if fid<0, disp(msg), return, end

 % Find header and first  data line ......................
is_head = 1;
jl = 0;
head = ' ';
while is_head  % Add lines to header.....
  s = fgets(fid);           % Get next line
  jl = jl+1;
  is_skip = jl<=skip;
  is_skip = jl<=skip | s(1)=='%';
  out1 = sscanf(s,formt);   % Try to read this line
   % If unreadable by SSCANF or skip, add to header
  is_head = isempty(out1) | is_skip;
  if is_head & ~is_skip
    head = str2mat(head,s(1:length(s)-1)); end
end
head = head(2:size(head,1),:);

 % Determine number of columns if not specified
out1 = out1(:)';
l1 = length(out1);
if ~nclm, nclm = l1; end

 % Read the rest of the file ..............................
if l1~=nclm  % First line format is different from ncolumns

  outdata = fscanf(fid,formt);
  lout = length(outdata)+l1;
  ncu = ceil(lout/nclm);
  lz = nclm*ncu-lout;
  outdata = [out1'; outdata(:); ones(lz,1)*addn];
  outdata = reshape(outdata,nclm,ncu)';

else              % Regular case

  outdata = fscanf(fid,formt,[nclm inf]);
  outdata = [out1; outdata'];  % Add the first line

end

fclose (fid);     % Close file ..........

Answer 2

[David Goodmanson]

Hi Timo,

From your attached data I can see that it's 2 sec at 48 kHz as you said. In your frequency plot, the data seems to be off by a factor of 10 for reasons I don't understand, but with the code below things come as they should, using your attached data.

There are N = 47999 data points. Hard to know what the value of L is, but it is really the length of the time record in milliseconds, that is not correct. In the context you are using it, L has to be the number of points in the fft which is N.

I took the second column of your attached data. Using the code below, the plot shows the 0 to 1kHz portion of the fft (absolute value). In the time domain, the data is basically a pulse train repeated at a frequency 18 Hz, plus noise. Since the fft of a pulse train is another pulse train, in the fft response we would expect to see a series of peaks at 18, 36, 54, 72, 90 etc. Hz. And that is exactly what we do see. The fact that the largest peak happens to be at 8 x 18 = 144 Hz is a less obvious result, having to do with details of the data. Eventually the higher frequency peaks get washed out by noise.

y = str2double(table2array(Data4));
y = y(:,2);
N = length(y);
fs = 48000;
delt = 1/fs;
t = delt*(0:N-1)';
fig(1)
plot(t,y)
Y = fft(y);
P2 = abs(Y/N);
P1 = P2(1:floor(N/2)+1);        % modified since N is odd
P1(2:end-1) = 2*P1(2:end-1);    % see below
f = fs*(0:floor(N/2))'/N;       % modified since N is odd
fig(2)
plot(f,P1)
xlim([0 1e3])

% this line has not yet been modified to work for both odd and even,
% but it's a small detail that does not affect the bottom line result at all