6 views (last 30 days)

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

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 ..........

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

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

Start Hunting!
## 5 Comments

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1127103

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1127103

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1127673

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1127673

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1127898

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1127898

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1128018

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1128018

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1128028

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/644828-why-do-i-get-too-high-dominant-frequencies-in-fft#comment_1128028

Sign in to comment.