FFT of Unsteady Temperature Data Resulting a Peak at 0 Hz

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));

[David Goodmanson]

Hi Oguzhan,

Your 50 Hz comment and the Fs = 100 in the code do not agree. I went with 100, so 100 Hz sampling for 60 sec but the basic conclusion is the same whichever one is used. Either way the result looks ok. I did the transform slightly differently, not using nexpow2, and saw very similar results. A lot of frequency content in the region of 0.1 Hz, which makes sense because that's 6 oscilllations in the time window and you can see a lot of that kind of thing, and slower, going on.

I'm not a fan of nextpow2 and I think it's not useful the vast majority of the time, but it looks to have done all right in this situation.

[Oguzhan M]

Hi David,

Thanks for your reply. Is Fs not supposed to be twice as much the sampling rate?

Can you also please furher explain what you meant vy " A lot of frequency content in the region of 0.1 Hz, which makes sense because that's 6 oscilllations in the time window and you can see a lot of that kind of thing, and slower, going on."

Many thanks.

[Daniel M]

No, sampling rate is the same as Fs. You are thinking of Nyquist rate, which is half of Fs, and is typically used for filtering (as you did with cut_off = 2 / (Fs/2) ).

For the second comment, he is saying that he can visually identify an underlying 0.1 Hz oscillation in the time-domain data. 0.1 Hz means one oscillation every 10 seconds.

Answer 1

[Daniel M]

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.

[Oguzhan M]

Hi Daniel,

Thanks for your reply. I'll definetely compare my results with the methods you suggest, which hopefully will be done in couple of hours.

Can you please share the codes you used for these plots? If not, no worries.

Many thanks.

[Daniel M]

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

[Oguzhan M]

Hi again Daniel,

I've been looking into details of the code you shared here for some time now, and I have few questions regarding results generated by semilogy and gradient plots.First, please see the attached plots (and data if you like). Also, please bear with me as I try to make sense out of all these.

1- Looking at two gradient plots, the dominant frequency for 250 RPM is 0.75, and for 600 RPM is 0.36. These results make perfect sense considering nature of the physics here. ( At specific measurement point for both measurements, 250 RPM is expected to show more oscillations than 600 RPM, so all good here.)

I wonder how to interpret the number of dominant frequency contents . As you can see, there are few peaks in 250 RPM whereas there are 2-3 distinct peaks in 600 RPM.

2- Why did you call gradient method as tough one?

3- Looking at two semilogy plots, what is the interpretation for the frequency domain plot? I can't see any distinguishable differences between 250 and 600 RPM as both frequency peaks are at 0 Hz.

4- Is there anything I need to pay attention to when using gradient code?

Many thanks.