# Fast fourier tranformer for Time series data

207 views (last 30 days)

Show older comments

Dear All.

I am a biginner and for the 1st time in my life I am programming an FFT code.

The problem:

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.

My objective:

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.

##### 10 Comments

Bjorn Gustavsson
on 8 Aug 2022

The data I get out of that xlsx-file looks nothing like the plot you've presented. What is going on?

### Answers (3)

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...

tData=readtable('Data.xlsx');

tData.T=tData.timeInS*1000;

plot(tData.T,tData.Displacement)

ylim([-0.520 -0.519])

xlim([tData.T(1) tData.T(100)])

Then explored the bad data sections...

find(tData.Displacement<-0.6)

xlim([tData.T(290) tData.T(320)])

ylim([-0.525 -0.519])

hold on

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

N1=848;tData.Displacement(N1:N1+1)=tData.Displacement([N1:N1+1]-4);

tData.D=tData.Displacement-mean(tData.Displacement); % remove DC

tData.T=tData.T-tData.T(1); % conver T to 0 origin

figure

plot(tData.T,tData.D)

xlabel('Time (msec)')

ylabel('Displacement from mean')

figure

plot(tData.T(1:100),tData.D(1:100)) % another look at the detail after detrending

xlim([0 tData.T(100)])

ylim([-4 4]*1E-4)

yline(0)

% now do a PSD

L = height(tData);

Y=fft(tData.D);

P2=abs(Y)/L;

P1=P2(1:L/2+1);

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

figure,tiledlayout(2,1)

nexttile

plot(f,P1)

xlabel('Frequency (Hz)')

ylabel('Displacement Magnitude')

nexttile

plot(f,convert2db(P1))

ylim([-140 -70])

xlabel('Frequency (Hz)')

ylabel('Displacement Mag (dB)')

##### 1 Comment

Bjorn Gustavsson
on 9 Aug 2022

Bjorn Gustavsson
on 8 Aug 2022

It seems that you're a bit too close to the Nyquist frequency. If you try:

[NUM,TXT,RAW]=xlsread('data.xlsx');

dt = mean(diff(NUM(:,5)));

fs = 1/dt;

% Calculate stft/spectrogram

[S,F,T] = spectrogram(NUM(1:1024,1),128,32,[],fs);

subplot(3,1,1)

plot(NUM(:,5),NUM(:,1),'.-')

hold on

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

plot(NUM(ilM,5),NUM(ilM,1),'r.')

plot(NUM(ilm,5),NUM(ilm,1),'c.')

subplot(3,1,2)

pcolor(T,F,log10(abs(S))),shading flat

subplot(3,1,3)

% look at the variability of the period-time between local max and local

% min:

plot(diff(ilm),'c.-')

hold on

plot(diff(ilM),'r.-')

From there you could also zoom in in the top panel to look at the signal.

HTH

##### 4 Comments

dpb
on 8 Aug 2022

I think if you don't remove the bum data spikes somehow they completely corrupt everything else -- if, instead of making the substitutions noted above before compute the PSD, just remove the mean of the initial time series, the result then looks like

The magnitude of those three spikes has so overwhelmed everything else in the signal that all you have is the artifacts from those three peaks alone. To prove that I then just zeroed out everything in the input signal excepting for those six points as

>> ixSpike=find(tD.Displacement<-0.6).';

ixSpike =

314 315 757 758 848 849

>> tD.DD(ixSpike)=tD.Displacement(ixSpike);

and redid the PSD with it and plotted on top of the above -- almost indistinguishable results showing that the actual data you were interested in is totally swamped by those three anomolous readings. I'm guessing that's causing much of the difficulties in the other interpretation, too.

dpb
on 8 Aug 2022

Edited: dpb
on 10 Aug 2022

Both of @Bjorn Gustavsson's comments above are true as well, I believe -- the data are definitely undersampled and there's enough jitter in the time sampling to have some effect -- both of those problems in the data collection should be fixed rather than trying to make something from the present data set which also suffers from the bad data that's really tough to replace given the above issues -- although even just putting in the mean is probably not nearly as bad as the effects of leaving it in.

Look up "antialising" -- if you don't have any at hand, I've a pair of Waveteks that have corner from 1 Hz to 100 kHz, each two channels. I've retired entirely from the consulting gig and would let them go for cheap; they've just been sitting on the bench for years, now...

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

Y=fft(tData.D,L2);

P2=abs(Y)/L;

P1=P2(1:L2/2+1);

f = (0:L2/2)*(Fs/L2);

P1(2:end-1)=2*P1(2:end-1);

[~,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;

if P1(P1il)<P1(P1iu)

P1il = find(P1(P1il+1:P1i-1)>=P1(P1iu),1,"first")+P1il;

elseif P1(P1il)>P1(P1iu)

P1iu = find(P1(P1i+1:P1iu-1)>P1(P1il),1,"last")+P1i-1;

end

figure

plot(f(1:P1il),P1(1:P1il),'b-',f(P1il:P1iu),P1(P1il:P1iu),'r-',f(P1iu:end),P1(P1iu:end),'b-')

Ypl = unwrap(angle(Y(P1i:-1:P1il)));

Yp = [Ypl(end:-1:2);unwrap(angle(Y(P1i:P1iu)))]';

hold on

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

plot(f(P1il:P1iu),Yp*Yps+Ypo,'g-',Xpo,Ypo,'r*')

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;

figure

plot(tData.T/1000,tData.D,Tx4,Ypo*cos(Zpo+Tx4*Xpo*2*pi))

##### 12 Comments

dpb
on 15 Aug 2022

Yeah, you'd already mentioned this -- am so engrained with expericence of using a tach independently to measure shaft speed and then collecting the vibration data with conventional multi-channel analyzers, still responded to @Jeffrey Clark's comment as if that were the case here as well.

The data then match the expectation as you note; what that implies is that one should definitely use the procedures appropriate for non-uniform sampling; it's probably not the best idea to try to resample back to a uniform time base since the DUT itself is actually somewhat variable.

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!