Why do I get zero frequencies in a fft when the constant offset is already removed?

The goal is to find the frequencies (if any at all) in the data. When performing the fft, no significant peaks can be distinguished and repetitive "humps" are seen in the fft. Additionaly, there seem to be a lot of peaks near 0 Hz, which could be caused by a constant offset. However this is also removed with S= fft( data -mean(data)). I am looking for frequencies betweeh 0 and 500 Hz, as I sample with a frequency of around 1300 Hz (the total sample time is 300 seconds, and around 390000 samples are made.

What can I do to find valuable peaks in the fft?

Thanks in advance

Answer 1

[Mathieu NOE]

hello

the fft seems to me inappropriate in this case. If the goal was to find the time difference between the succesive pulses , see the %%% NEW CODE %% I added at the bottom of your original code

this new portion of code will define the time stamps for positive and negative slope (aka rising / falling) trigger.

results are given by :

mean_period_pos = mean(diff(t0_pos)) % (seconds) /  mean period computed on positive slope crossings 
mean_period_neg = mean(diff(t0_neg)) % (seconds) /  mean period computed on negative slope crossings 

results

mean_period_pos = 21.5847 (seconds)

mean_period_neg = 21.6408 (seconds)

As we can see the period is more or less the same if we use the rising or falling slopes of the signal

full code

FYI I changed for another smoothing function as movmean is not in my package. You can of course go back to your own version

just make sure that the signal is clean were you want the trigger to operate (here at y = 0.5) . if your data has zigzags all over the place due to noise , the crossing function will generate uneeded extra trigger points

clear variables; close all; clc;

%% Uploading the data from the concentration measurements
data = readtable('Test20_Position3_2kg_90d_5mm_InducedVibration20s_09_06_2022');
t1 = table2array( data(:,2) )./1000000 - table2array(data(1,2))./1000000;
I1 = table2array( data(:,1));

%% linearising the data
I1_I0_1= (10.^((I1-37)/(458*0.5))-1)/100;  
%% This section will remove the present noise, via a moving average

% Window = 21; %number of periods for the moving average
% Avrg_1 = movavg(I1_I0_1,'simple',Window);
Avrg_1 = smoothdata(I1_I0_1,'gaussian',100); % my code 

%% Plotting all data: Laser signal I over time

figure()
hold on
plot(t1,Avrg_1)
legend('Position 3')
title('Laser signal I over time:Induced Pulsations each 20 seconds')
xlabel('Time [s]'); 
ylabel('I [-]')
grid on

%%
% Discrete Fourier Transform (DFT)
S = fft(Avrg_1);      
% S = fft(detrend(Avrg_1));      
Fs1 = length(Avrg_1)/300;            % Sampling frequency    
T = 1/Fs1;                           % Sampling period 
L1 = length(Avrg_1);                 % Length of signal         
f1 = Fs1*(0:(L1/2))/L1;        


%% Converting magnitude to db and taking the left part of the mirrored Fourier transform
%Y = fft(X) computes the discrete Fourier transform (DFT) of X using a fast Fourier transform (FFT) algorithm.

%Left side of the spectrum
P2 = abs(S/L1);
P1 = P2(1:L1/2+1);                     
P1(2:end-1) = 2*P1(2:end-1);

%Removing lowest frequencies
for i = 1:length(P1)
    if P1(i) < 10^(-10)
        P1_angle(i) = 0;
    end
end 


%% FFT Plot:     Magnitude (Db) versus frequency (Hz)
figure()
% hold on
semilogy(f1,P1)                      %Frequency plotted by Fourier Transform

title('Single-Sided Amplitude Spectrum: Position 5, 90 degrees')
xlabel('f (Hz)')
ylabel('|P1(f)|')
legend('Data1')
% xlim([0 600])
% ylim([0 0.0003])

%%% NEW CODE %%%

threshold = 0.5; % your value here

[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(Avrg_1,t1,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

figure()
plot(t1,Avrg_1,'b',t0_pos,s0_pos,'*r',t0_neg,s0_neg,'*g','linewidth',2,'markersize',12);grid on
legend('signal','signal positive slope crossing points','signal negative slope crossing points');
% subplot(2,1,2),plot(t_period,period,'linewidth',2,'markersize',12);grid on

mean_period_pos = mean(diff(t0_pos)) % (seconds) /  mean period computed on positive slope crossings 
mean_period_neg = mean(diff(t0_neg)) % (seconds) /  mean period computed on negative slope crossings 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%


% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if ~isempty(ind)

    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)

                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end

    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);

    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end

[Mathieu NOE]

Just a few remarks

so the pulses are time distant of dt = 21.6 seconds , which is a frequency of 0.0462 Hz

as you have noticed, your frequency vector generated by the fft has resolution = mean(diff(f1)) = 0.0033 Hz

the good point is that the fft resolution is small enough to show the pulses rate and harmonics, but you really needed to zoom the plot down to 0 to 1 Hz to see it clearly.

if you now zoom back , you will have more and more difficulties to distinguish the higher harmonics , but at least you can see their amplitude goes down as with frequency

here for 0 to 10 Hz

and so forth...until you see only a "continuous" spectrum with no peaks , but there are present but very close , by definition their frequency distance equals the fundamental frequency = 0.0462 Hz

the fact that you have a strong DC output from the fft comes also that we have a baseline which is not strictly dc constant all the time, it's more a piecewise "almost" constant with some very low frequency rumble. That also create the strong fft output below 0.04 Hz

[Tyl van Gestel]

Dear Mathieu,

Thank you for the extended answer!

However, I'm mostly interested in finding the frequencies between the induced pulsations. As in your fft, there are no clear high frequency peaks, unless you zoom in far enough. The reason why I induced the pulses after dt = 21.6 seconds (thank you for the code defining the exact time stamps), is to see whether I could create distinguishable frequencies with higher magnitudes. Unfortunately, the magnitude of frequency of the induced pulsations is not clearly different than other frequency peaks.

I added another data set without the induced pulsations:

Which creates the following fft:

I'm wondering what the physical meaning is of the repetitive "humps", and if the top of the hump says something about occuring frequencies in the measurement. Also, when adjusting the window of the Moving Average filter (smoothdata( [ ], [ ], 100) these humps are removed when increasing or decreasing the window size, for example for 10 or 1000. Do you know more about the meaning of these repetitive humps?

[Mathieu NOE]

hello Tyl

so this new data sets is your "background" noise ? I have no idea what is the physical process behind it , but it generates some cyclic patterns with a spectral content that has itself a cyclic content. Does it mean your system has a physical processing system that is repetitive by design ?

[Tyl van Gestel]

Dear Mathieu,

The data sets are measured data and naturally contain some noise. However, the question is how catogorise the measured data as "true" data and noise. This is more a problem of my applied measuring method which I cannot ask your help for. In any case, thank you a lot for the support around the fast fourier transform.

Kinds regards,

Tyl

[Mathieu NOE]

hello Tyl

OK thank you for the explanations

good luck for the future !