Problem with Fourier Transform fft

Hello,

I have field data (Fig 1) and I want to apply a Fourier transform of the siganl signal. Here is the code I use:

load('dataset')
time_duration=(Time-Time(1))* 24 * 60 * 60;
Time_Continuous=0:60:(Time(end)-Time(1))* 24 * 60 * 60;
Data_Signal_interp=interp1(time_duration,Data_Signal,Time_Continuous); %Because the data set is not continuous in time
Fs = 1/60;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = time_duration(end)/T+1;             % Length of signal
t = (0:L-1)*T;        % Time vector
X=Data_Signal_interp;
Y = fft(X);
S_oneside=Y(1:L/2);
f=Fs*(0:L/2-1)/L;
S_meg=abs(S_oneside)/(L/2);
plot(f,S_meg)

The result is not what I expected (Fig 2). Could you please help me about that ?

Thanks a lot!

NS

[Mathieu NOE]

hello

what do you expected ?

[NASI]

hello Mathieu,

The raw signal clealry shows variations at different frequencies (around 1 day for instance). Here the Fourier transform does not highlight these variations

[Matt J]

open Fig1.fig

open Fig2.fig

Answer 1

[Star Strider]

The ‘Time’ vector is not uniformly sampled. Most signal processing techniques require regularly-sampled vecttors, so I used the resample function to achieve that here.

Try something like this —

LD = load(websave('dataset','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1291920/dataset.mat'))

LD = struct with fields:
    Data_signal: [1.0048e+03 1.0051e+03 1.0053e+03 1.0054e+03 1.0057e+03 1.0058e+03 1.0058e+03 1.0058e+03 1.0058e+03 1.0059e+03 1.0062e+03 1.0061e+03 1.0059e+03 1006 1006 1.0061e+03 1.0062e+03 1.0061e+03 1.0062e+03 1.0063e+03 1.0062e+03 1.0061e+03 … ]
           Time: [7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 7.3886e+05 … ]

Data_signal = LD.Data_signal;
Time = LD.Time;
CheckTime = [mean(diff(Time)) std(diff(Time))]

CheckTime = 1×2

0.0024 0.0014

Fs = ceil(1/mean(diff(Time)))

Fs = 420

Fn = Fs/2;
[Data_signalr, Timer] = resample(Data_signal, Time, Fs);

figure
plot(Timer, Data_signalr)
grid

L = numel(Timer)

L = 21727

NFFT = 2^nextpow2(L)

NFFT = 32768

FTs = fft((Data_signalr(:)-mean(Data_signalr)).*hann(L), NFFT);
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

[pks,locs] = findpeaks(abs(FTs(Iv))*2, 'MinPeakProminence',1E+4)

pks = 6×1

1.0e+04 * 6.8170 9.3346 3.1670 1.8018 1.7125 1.0741

locs = 6×1

5 10 23 34 44 57

figure
plot(Fv, abs(FTs(Iv))*2)
hold on
plot(Fv(locs), pks, '^r')
hold off
grid
xlim([0 5])
ylim(ylim*1.5)
text(Fv(locs), pks, compose('\\leftarrow Mag = %.2f, Per = %0.3f', [pks 1./Fv(locs).']), 'Rotation',45)

EDIT — (10 Feb 2023 at 17:54)

Changed text call to calculate and display period rather than frequency (original code).

.

[NASI]

Hello,

Thank you about that !

Just one question. The frequency is built from Matlab Time (date numbers), right ?

In this case, what is the x-axis ? Is it possible to get it in frequency (Hz) ? I tried to convert 'Time' in sec but in this case the result is not the same at all.

I'd be grateful if you could help me with this.

NS

[Star Strider]

My pleasure!

‘The frequency is built from Matlab Time (date numbers), right ?’

That depends on how you extracted it. It is calculated initially from the ‘Time’ vector provided in ‘dataset.mat’ and then, since the standard deviation of the time intervals (0.0014) is quite high with respect to the mean of the time intervals (0.0024), meaning that there is a significant variation in the time intervals, I resampled them (using the resample function I linked to earlier) to a sampling frequency of 420, so if the ‘Time’ vector is in seconds, that sampling frequency (and all the frequencies in the fft plot) will be in Hz. (Without knowing the units of ‘Time’ the frequency will be in . Here, that is Hz if the time units are in seconds.)

‘In this case, what is the x-axis ?’

The x-axis of the Fourier transforrm plot is in Hz if the time units are in seconds. Otherwise it is in

.

‘Is it possible to get it in frequency (Hz) ?’

If the ‘Time’ units are seconds, it already is. If they are not, then some sort of transformation to seconds will be needed first. In that instance, the transformation will have to be done first, and then the rest of the code will need to be run with the transformed time vector.

.

[NASI]

Thank you very much !

All is clear.

Last question (I promise), the unit of amplitude (yaxis) is given by the unit of the data ? or its square ? (since you plot plot(Fv, abs(FTs(Iv))*2))

[Star Strider]

My pleasure!

In this instance, the magnitude of the peaks are approximately equal to the amplitude of the time domain waveform relative to the mean of the signal. It is not the square because neither the signal itself or the Fourier transform of it are squared. (The square would be the power rather than the magnitude of the signal at each frequency.) The reason that it is not exactly equal has to do with the frequency resolution and to ‘spectral leakage’ of energy into neighbouring frequencies near the centre frequency due to the finite nature of the signal and of the fft calculations. Windowing the signal helps minimise it, although it does not completely prevent it.

I subtracted the mean of the signal here before calculating the Fourier transform because otherwise the mean magnitude at 0 frequency would make the other peaks very difficult to see in the plot. Subtracting the mean does not otherwise affect the other peak magnitudes.

Answer 2

[Mathieu NOE]

hello

when you do the fft with a buffer length equal your signal length (as you did) you get the maximal frequency resolution (nice) but periodic content in your signal may be difficult to see (very little emergence) as you don't average the fft process (here the signal will be cut in multiple chuncks (buffers) with some overlap to reduce the effect of noise)

here you can see the plot that overlay your fft (no averaging) and mine (with lot of averaging). what you lose in terms of frequency resolution is balanced with a better amplitude emergence, but that depends also if the signal is more or less stationnary. You have to try and compare the results.

beside that , if you look at the data you generated after interp1, it's a very long vector for a fairly low frequency oscillating signal. So if you are interested mostly in those low frequency waves , there is no need to resample the data so fast with T = 60

I probably get a better result if a choose a lower "sampling" rate , but not too low otherwise I will loose some valuable information (the interpolated signal will look too coarse)

here I changed from T = 60 to T = 6000 (yes !) so you can see that the fft frequency vector has also been reduced to match better the low frequency domain (max value of freq vector is also reduced by factor 100

and the original and "low FS" resampled time data do overlay quite well even with T = 6000

so , again, no need to resample super fast withT = 60.

Now , this is the cherry on the cake , I let you play with the averaged fft vs "one shot" fft to see the effect of NFFT and OVERLAP

load('dataset')
Time = Time-Time(1); % start Time at zero
time_duration =Time* 24 * 60 * 60;

T = 6000; % now I can change the sampling rate (60 => 6000)
Time_Continuous=0:T:(Time(end)-Time(1))* 24 * 60 * 60;
x=interp1(time_duration,Data_signal,Time_Continuous); %Because the data set is not continuous in time

figure(1),plot(time_duration,Data_signal,Time_Continuous,x);
legend('original','resampled');

Fs = 1/T;            % Sampling frequency                    
L = time_duration(end)/T+1;             % Length of signal
t = (0:L-1)*T;        % Time vector

Y = fft(x);
S_oneside=Y(1:L/2);
f=Fs*(0:L/2-1)/L;
S_meg=abs(S_oneside)/(L/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : "one shot" vs averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 256;
OVERLAP = round(0.95*NFFT);
[freq, spectrum_raw] = myfft_peak(x,Fs,NFFT,OVERLAP);

figure(2),semilogy(f,S_meg,'b',freq,spectrum_raw,'r');grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude');
legend('fft','averaged fft');

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

if channels> 10*samples  % data must be transposed
    signal = signal';
end
[samples,channels] = size(signal);


% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

[Mathieu NOE]

A "one per day" oscillation would mean a peak at f = 1.1574e-05 Hz

On the fft plot, I can see one emerging peak at 2.3e-5 Hz so twice as fast => 2 oscillations per day