Problem with Fourier Transform fft

22 views (last 30 days)
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
  3 Comments
NASI
NASI on 10 Feb 2023
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

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 10 Feb 2023
Edited: Star Strider on 10 Feb 2023
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).
.
  4 Comments
NASI
NASI on 27 Feb 2023
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
Star Strider on 27 Feb 2023
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.

Sign in to comment.

More Answers (1)

Mathieu NOE
Mathieu NOE on 10 Feb 2023
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
  1 Comment
Mathieu NOE
Mathieu NOE on 10 Feb 2023
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

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!