MATLAB Answers

Why do I get too high dominant frequencies in FFT

6 views (last 30 days)
Timo Mayer
Timo Mayer on 11 Nov 2020
Edited: David Goodmanson on 12 Nov 2020
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

  5 Comments

Show 2 older comments
David Goodmanson
David Goodmanson on 12 Nov 2020
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
Timo Mayer on 12 Nov 2020
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
David Goodmanson on 12 Nov 2020
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.

Sign in to comment.

Answers (2)

Mathieu NOE
Mathieu NOE on 12 Nov 2020
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, kirill@plume.mit.edu
% 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 ..........

  0 Comments

Sign in to comment.


David Goodmanson
David Goodmanson on 12 Nov 2020
Edited: David Goodmanson on 12 Nov 2020
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

  0 Comments

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!