5 views (last 30 days)

Hi all,

I have a data that was sampled at 50 Hz for 120 seconds. When I apply fft, dominant frequency seems to be 0 Hz which I initially thought it could actually very well be. However, when I apply same fft to different set of data acquired from the experiment, I get the same peak at 0 Hz, which does not seem to be correct as I expect it to show relatively high frequency(ies) ( this is due to nature of experiment in which the unsteady temperature data is acquired, and the second set of data is expected to show some form of periodicity).

As you may see below in the code, DC offset is removed by "detrend" and a low pass filter is applied ( the result was the same even before low pass filtering).

Can you please take a look at the code (data is attached as well) for any possible corrections or recommendations? If everything is okay with the code, I will appreciate some comments on the physical meaning of the frequency analysis result in this particular case.

Best.

Fs = 100; % Sampling frequency

T = 1/Fs; % Sampling period

dt = 0 :T:120-T; % Time vector

nfft= length(S2TR3_0); %Length of FFT

nfft2 = 2.^nextpow2(nfft) ;

S2TR3_0 = detrend(S2TR3_0);

figure

plot(S2TR3_0)

fy = fft(S2TR3_0,nfft2);

fy = fy(1:nfft2/2);

xfft = Fs.*(0:nfft2/2-1)/nfft2;

plot(xfft,abs(fy/max(fy)));

%low pass filter

cut_off = 2/Fs/2;

order = 256;

h = fir1(order,cut_off);

con = conv(S2TR3_0,h);

figure

plot(con);

fh = fft(h,nfft2);

fh = fh(1:nfft2/2);

fh = fh';

mul = fh.*fy;

figure

plot(abs(mul));

Daniel M
on 25 Apr 2020

Edited: Daniel M
on 25 Apr 2020

Data looks fine. You shoud plot using semilogy to better visualize the higher frequencies. Your data shows a typical 1/f noise pattern, seen in white and pink noise dominated data. (Fig 1).

You can also try removing this 1/f curve by taking the gradient of the data, this is a rough method. (Fig 2)

Compare these semilogy plots with the periodic data and see if you can spot the differences now.

Daniel M
on 25 Apr 2020

I've attached three files for you. Put them all on your path and you can call them from anywhere.

Create the plots using this line of code

plotFFT(S2TR3_0,Fs,[],[],0,'semilogy',1);

plotFFT(gradient(S2TR3_0),Fs,[],[],0,'plot',1);

I'm also going to paste the contents of the plotFFT function, for the casual browser.

function [H,F,Y] = plotFFT(x,fs,begtime,endtime,demean,plotType,side,varargin)

% [H,F,Y] = plotFFT(x,fs,begtime,endtime,demean,plotType,side,varargin)

% This function takes a signal, x, and the sample rate, fs, and

% plots the FFT between begtime and endtime (defaults to

% beginning and end of signal). Times are in seconds, not

% indices.

%

% Inputs:

% x - must be a signal from one channel

% fs - the sampling frequency

% begtime/endtime - beginning and end in seconds. default is full amount

% demean - 1 or 0 to detrend data. 0 is default

% plotType - 'plot' (default), 'semilogy', 'semilogx', 'loglog'

% side - input to figurebig()

% varargin - plotting options (color, linestyle, etc.)

%

% Outputs

% H - figure handle

% F - one-sided frequency vector

% Y - one-sided Fourier transform

T = 1/fs; % sampling period

flipBack = false;

%%% Do some input error checking

if size(x,1) ~= 1

if size(x,2) == 1

% make it a row vector

x = x.';

flipBack = true; % return the outputs as column vectors

else

error('x must be a vector')

end

end

if nargin < 3 || isempty(begtime) || begtime < 0

begtime = 0;

end

if nargin < 4 || isempty(endtime) || endtime > (length(x)-1)*T

endtime = (length(x)-1)*T;

end

if begtime >= endtime

error('begtime must be before endtime')

end

if nargin < 5 || isempty(demean)

demean = 0;

elseif ~(demean==1 || demean==0)

error('demean must be true or false')

end

if nargin < 6 || isempty(plotType)

plotType = 'plot';

elseif ~any(strcmp(plotType,{'plot','semilogy','semilogx','loglog'}))

error('Incorrect plotType')

end

if nargin < 7 ||isempty(side)

side = 0;

end

% check if user passed 'Visible' into varargin

% [~,visloc] = ismember('Visible',varargin);

% if visloc == 0; visloc = []; end % avoid error in figurebig

visloc = find(strcmp('Visible',varargin)); % ismember errors if varargin is not cellstr

% Get indices from beg/end times

begind = floor(begtime*fs+1);

endind = floor(endtime*fs+1);

X = x(begind:endind); % new signal

if demean

X = detrend(X);

end

L = length(X); % length of signal

N = L; % n-point fft. This could be an input.

Y2 = fft(X,N); % fourier transform, 2-sided

t = (0:L-1)*T; % time vector

% Create one-sided frequency vector including Nyq freq

% This is the same for both even and odd numbered vectors

Ny = fix(N/2)+1; % Nyq index

F = linspace(0,1,Ny)*fs/2; % frequency vector, 1-sided

Y2 = abs(Y2/N); % 2-sided, scaled by N

Y = [Y2(1) 2.*Y2(2:Ny)]; % get positive freqs, scale power by 2 for non-dc components

%%% Make the plot

H = figurebig(side,varargin{[visloc,visloc+1]});

subplot(3,1,1)

plot(t,X,varargin{:})

grid minor

xlabel('Time [s]')

ylabel('Signal')

axis tight

subplot(3,1,[2 3])

switch plotType

case 'plot'

plot(F,Y,varargin{:})

case 'semilogy'

semilogy(F,Y,varargin{:})

case 'semilogx'

semilogx(F,Y,varargin{:})

case 'loglog'

loglog(F,Y,varargin{:})

end

xlabel('Frequency [Hz]')

ylabel('|P(f)|')

axis tight

grid on

H = tightfig(H);

if flipBack

% x was input as a column vector. Return F and Y as column

% vectors

F = F.';

Y = Y.';

end

end

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

Start Hunting!
## 3 Comments

## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/520568-fft-of-unsteady-temperature-data-resulting-a-peak-at-0-hz#comment_834807

⋮## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/520568-fft-of-unsteady-temperature-data-resulting-a-peak-at-0-hz#comment_834807

## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/520568-fft-of-unsteady-temperature-data-resulting-a-peak-at-0-hz#comment_834907

⋮## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/520568-fft-of-unsteady-temperature-data-resulting-a-peak-at-0-hz#comment_834907

## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/520568-fft-of-unsteady-temperature-data-resulting-a-peak-at-0-hz#comment_835023

⋮## Direct link to this comment

https://es.mathworks.com/matlabcentral/answers/520568-fft-of-unsteady-temperature-data-resulting-a-peak-at-0-hz#comment_835023

Sign in to comment.