Fast fourier tranformer for Time series data

84 次查看（过去 30 天）
Mohammed Lamine Mekhalfia
评论： dpb 2022-8-15
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 个评论
显示 8更早的评论
Bjorn Gustavsson
Bjorn Gustavsson 2022-8-8
The data I get out of that xlsx-file looks nothing like the plot you've presented. What is going on?
Mohammed Lamine Mekhalfia
Mohammed Lamine Mekhalfia 2022-8-8
@Bjorn Gustavsson
attached the signal image, there is a 3 points that are divergening from the values, so what I did is to interpolate those values only.
and then I continue with the written code

请先登录，再进行评论。

请先登录，再回答此问题。

回答（3 个）

dpb
dpb 2022-8-8
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 个评论
Bjorn Gustavsson
Bjorn Gustavsson 2022-8-9
To add to @dpb's answer you might get slightly better results if you use nufft since that function sould take the jitter in the sampling into account. But none of our suggestions will overcome th fundamental problem that you're a bit too close to the Nyquist-frequency.

请先登录，再进行评论。

Bjorn Gustavsson
Bjorn Gustavsson 2022-8-8
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 个评论
显示 2更早的评论
dpb
dpb 2022-8-8
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
dpb 2022-8-8
编辑：dpb 2022-8-10
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
Jeffrey Clark 2022-8-10
@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 个评论
显示 10更早的评论
Mohammed Lamine Mekhalfia
Mohammed Lamine Mekhalfia 2022-8-15
@Jeffrey Clark @Bjorn Gustavsson @Bjorn Gustavsson
to answer your questions, the time is capture from the omce per revolution sensor named in the figure above as the speed sensor, it is used to measure the time of one rotation and at each revolution we capture also the time the blade cross the Bj sensors in the front
so time is from speed sensors ( it has some difference because we used a very small time scale were the speed flectuation is obvious, we were measureing at a constant speed yet there must be always an error gap due to speed controllers ...etc so the value reflect the real measured value).
for the displacement it is the difference between the time of blade cross the sensor - time of 1 rotation.
dpb
dpb 2022-8-15
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Whos 的更多信息

标签

产品


版本

R2018a

Community Treasure Hunt

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

Start Hunting!

Translated by