Fast fourier tranformer for Time series data
207 views (last 30 days)
I am a biginner and for the 1st time in my life I am programming an FFT code.
well I aquired data ( 1052 values ) each values captured at a time period ( it is not perfectly uniform ), this time period equalt to nearly 0.01132 s (+/- 0.000005 s)
the vector of time is recorded also and it contain 1052 values as explained. below the plot of the time data signal.
Plot the FFT of this DATA to identify the spectrum.
What I did so far:
I followed this page to code https://www.mathworks.com/help/matlab/ref/fft.html, yet I did not have anything and this due that I am not understanding the process on how to do it.
Thanks in advance for all of your thoughts and suggestions.
dpb on 8 Aug 2022
No disagreements at all with what @Bjorn Gustavsson notes, just do a baseband PSD and look at the signal some...
First, look at a short section of the sampled waveform in enough detail to see it shows us --
which indicates as I suspected that the signal is undersampled for a best estimate of amplitude -- note the peaks are often flattened-of approximating square waves and most are just a spike; a truly adequate sample rate would show those as more rounded peaks; the actual magnitude of a the waveform at the peak will almost always be at a time instant that wasn't sampled.
Second, for the points that were saturated -- what caused that to happen, btw?, ought to try to fix when collect some additional data with better technique/higher sampling rate -- I looked at a region of that as well -- and just arbitrarily filled in an adjacent section -- trying to interpolate such an undersampled signal is pretty-much a pure guess, anyways...
Did same thing for other locations as well...
That done, I subtracted the mean from the resulting time signal to remove the DC component and did a one-sided PSD with the mean sampling frequency.
The top is on linear scale, the lower in dB; there is a strong peak at 28.4 Hz; not sure what that would be related to, but there is only the one dominant frequency in the signal; everything else is just noise it would appear.
For further work, you need to either increase the sample rate to something approaching 150 Hz or so (3-4X the highest frequency of interest is a common recommendation) or incorporate the use of analog lowpass filters before the A/D to limit the higher-frequency input energy.
I've no klew how/what use you are trying to make of this; I don't see how it relates to the problem description at all, but there are certainly clever techniques for solving problems I'm not aware of -- do you have a reference link to the method in application?
The above was done with code that (roughly) matches the following (it's not exact, I just poked around at the command line and then pulled out pertinent pieces to attach...)
% examine narrow range to ensure not aliased badly...
Then explored the bad data sections...
hL(1)=plot(tData.T(313:316),tData.Displacement([313 310 311 316]),'x-r'); % the plot presented...
tData.Displacement(314:315)=tData.Displacement(310:311); % subsitute for bad
N1=757;tData.Displacement(N1:N1+1)=tData.Displacement([N1:N1+1]-4); % ditto
tData.D=tData.Displacement-mean(tData.Displacement); % remove DC
tData.T=tData.T-tData.T(1); % conver T to 0 origin
ylabel('Displacement from mean')
plot(tData.T(1:100),tData.D(1:100)) % another look at the detail after detrending
% now do a PSD
L = height(tData);
Fs=1/mean(diff(tData.T/1000)); % mean sample rate
f = (0:L/2)*(Fs/L); % compute frequency vector to go with...
P1(2:end-1)=2*P1(2:end-1); % one-sided PSD is half total energy except DC and Fmaxc
ylabel('Displacement Mag (dB)')
Bjorn Gustavsson on 8 Aug 2022
It seems that you're a bit too close to the Nyquist frequency. If you try:
dt = mean(diff(NUM(:,5)));
fs = 1/dt;
% Calculate stft/spectrogram
[S,F,T] = spectrogram(NUM(1:1024,1),128,32,,fs);
% identify the local max and min
ilM = find(NUM(1:end-2,1)<NUM(2:end-1,1)&NUM(2:end-1,1)>NUM(3:end,1))+1;
ilm = find(NUM(1:end-2,1)>NUM(2:end-1,1)&NUM(2:end-1,1)<NUM(3:end,1))+1;
% look at the variability of the period-time between local max and local
From there you could also zoom in in the top panel to look at the signal.
Jeffrey Clark on 10 Aug 2022
@Mohammed Lamine Mekhalfia, I think the answer provided by @dpb is valid although @Bjorn Gustavsson has valid points as well. By better approximation of the peak actual value in amplitude, frequency and phase a cosine can be overlayed on the signal and has good correlation with the mean adjusted measurements in both where the "noise" seems strong and where mostly absent:
Their points about the missing data, slight time variances, and probable noise level do need to be addressed in future samples. The noise could be mitigated by a much longer sample period. For the record my calculations building upon dpb's code produced:
- amplitude: 2.0559e-04
- frequency: 28.3870 no chirp detectable
- phase: -2.6946 at time=0
%Compute cosine model of signal using dpb's code given in another answer
L2 = 2^(nextpow2(L)+1);
f = (0:L2/2)*(Fs/L2);
[~,P1i] = max(P1);
P1il = find(P1(1:P1i-1)>P1(2:P1i),1,"last")+1;
P1iu = find(P1(P1i+1:end)>P1(P1i:end-1),1,"first")+P1i-1;
P1il = find(P1(P1il+1:P1i-1)>=P1(P1iu),1,"first")+P1il;
P1iu = find(P1(P1i+1:P1iu-1)>P1(P1il),1,"last")+P1i-1;
Ypl = unwrap(angle(Y(P1i:-1:P1il)));
Yp = [Ypl(end:-1:2);unwrap(angle(Y(P1i:P1iu)))]';
Yps = max(abs(diff(P1(P1i-1:P1i+1))))/(max(Yp)-min(Yp));
csP1 = cumsum(P1(P1il:P1iu).^2);
pP1 = polyfit(f(P1i-1:P1i+1),P1(P1i-1:P1i+1),2);
pYp = polyfit(f(P1i-1:P1i+1),Yp([P1i-1:P1i+1]-P1il+1),2);
Xpo = -pP1(2)/(2*pP1(1));
Ypo = polyval(pP1,Xpo);
Zpo = polyval(pYp,Xpo);
axis([f([P1il-1 P1iu+1]) min(P1(P1il-1:P1iu+1)) max(0,max(Yp)*Yps)+Ypo*1.1])
Tx4 = [0:0.25:length(tData.T)]/Fs;