FFT command on discrete signals and processing the results
in [Matlab]
|
From: Jos on 27 Jul 2010 05:58 Hi all, Currently I'm working on a project that involves an analysis in a FE package which calculates the loads on a part of our structure due to dynamic responses. This is stored as a load-time history file. I want to investigate the FFT and the PSD of this signal to eventually see if the signal contains any non-predicted frequencies and to check if the signal is narrow- or widebanded. (any input on the last subject is highly appreciated) The following parameter can be provided. Sample time is 0.5 sec. Total time history is 3hours Period of the inputsignal (sea state) is 5.5 sec. *Probably varied in the FEA program* I store the data in the vector 'sum_hor'. After investigating the files in the help and online I came to the following m-file. The question to you guys is if I'm making any mistakes (which I assume I do).. ------------------ H = sum_hor; nfft = 2^nextpow2(length(H)); Pxx = abs(fft(H,nfft)).^2 / length(H); %FFT figure(1) plot(Pxx(1:length(Pxx)/2) Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2)); %PSD figure(2) plot(Hpsd) ------------------ In addition, the results generated by this m-file are a bit noisy, any suggestions how to cope with this problem? Thanks in advance for all the help, Jos
From: Wayne King on 27 Jul 2010 06:34 "Jos " <j.g.h.haarman(a)student.utwente.nl> wrote in message <i2majc$462$1(a)fred.mathworks.com>... > Hi all, > > Currently I'm working on a project that involves an analysis in a FE package which calculates the loads on a part of our structure due to dynamic responses. This is stored as a load-time history file. I want to investigate the FFT and the PSD of this signal to eventually see if the signal contains any non-predicted frequencies and to check if the signal is narrow- or widebanded. (any input on the last subject is highly appreciated) > > The following parameter can be provided. > Sample time is 0.5 sec. > Total time history is 3hours > Period of the inputsignal (sea state) is 5.5 sec. *Probably varied in the FEA program* > > I store the data in the vector 'sum_hor'. > > After investigating the files in the help and online I came to the following m-file. The question to you guys is if I'm making any mistakes (which I assume I do).. > > ------------------ > > H = sum_hor; > nfft = 2^nextpow2(length(H)); > Pxx = abs(fft(H,nfft)).^2 / length(H); > > %FFT > figure(1) > plot(Pxx(1:length(Pxx)/2) > > Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2)); > > %PSD > figure(2) > plot(Hpsd) > > ------------------ > > In addition, the results generated by this m-file are a bit noisy, any suggestions how to cope with this problem? > > Thanks in advance for all the help, > Jos Hi Jos, The easiest way to arrive at a PSD estimate is using spectral analysis objects. Assumming x is your signal. The following will compute a Welch spectral estimate and periodogram (the periodogram is what you've done above). You will note that the Welch estimate exhibits markedly less variance than the periodogram. I'm assuming that you use the entire 3 hour recording. Of course, it may be beneficial to do a time-frequency analysis on that data if the statistical properties change over time. For this example, I'll assume the data is wide-sense stationary and use the entire data record. Fs = 2; % your sampling interval is 0.5 t = 0:1/Fs:(3600*3); % signal x = cos(2*pi*(1/5.5)*t)+randn(size(t)); % compute a Welch estimate % Hamming window, segment length of 500 samples hspec = spectrum.welch('Hamming',500); PsdEst = psd(hspec,x,'Fs',2); plot(PsdEst); % compute a periodogram hp = spectrum.periodogram; Psdper = psd(hp,x,'Fs',2); figure; plot(Psdper); I've used a segment length of 500. If you reduce this length, the variance of the PSD estimate will reduce further, but the width of the peak will broaden. Note that the Welch estimate is much less "noisy" than the periodogram, but this noise reduction comes with a loss of resolution. Another way would be to reduce the noise in the estimate is to first low pass filter your input. If your frequency of interest is 1/5.5 hertz, you can remove everything above a certain frequency, let's say 0.25. hfilt = fdesign.lowpass('Fp,Fst,Ap,Ast',0.25,0.3,0.1,50,2); d = design(hfilt); % filter the input data y = filter(d,x); % if you need zero-phased filtering y1 = filtfilt(d.Numerator,1,x); Hope that helps, Wayne
From: Jos on 27 Jul 2010 10:00 "Wayne King" <wmkingty(a)gmail.com> wrote in message <i2mcmr$k20$1(a)fred.mathworks.com>... > "Jos " <j.g.h.haarman(a)student.utwente.nl> wrote in message <i2majc$462$1(a)fred.mathworks.com>... > > Hi all, > > > > Currently I'm working on a project that involves an analysis in a FE package which calculates the loads on a part of our structure due to dynamic responses. This is stored as a load-time history file. I want to investigate the FFT and the PSD of this signal to eventually see if the signal contains any non-predicted frequencies and to check if the signal is narrow- or widebanded. (any input on the last subject is highly appreciated) > > > > The following parameter can be provided. > > Sample time is 0.5 sec. > > Total time history is 3hours > > Period of the inputsignal (sea state) is 5.5 sec. *Probably varied in the FEA program* > > > > I store the data in the vector 'sum_hor'. > > > > After investigating the files in the help and online I came to the following m-file. The question to you guys is if I'm making any mistakes (which I assume I do).. > > > > ------------------ > > > > H = sum_hor; > > nfft = 2^nextpow2(length(H)); > > Pxx = abs(fft(H,nfft)).^2 / length(H); > > > > %FFT > > figure(1) > > plot(Pxx(1:length(Pxx)/2) > > > > Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2)); > > > > %PSD > > figure(2) > > plot(Hpsd) > > > > ------------------ > > > > In addition, the results generated by this m-file are a bit noisy, any suggestions how to cope with this problem? > > > > Thanks in advance for all the help, > > Jos > > Hi Jos, The easiest way to arrive at a PSD estimate is using spectral analysis objects. > > Assumming x is your signal. The following will compute a Welch spectral estimate and periodogram (the periodogram is what you've done above). You will note that the Welch estimate exhibits markedly less variance than the periodogram. I'm assuming that you use the entire 3 hour recording. Of course, it may be beneficial to do a time-frequency analysis on that data if the statistical properties change over time. For this example, I'll assume the data is wide-sense stationary and use the entire data record. > > Fs = 2; % your sampling interval is 0.5 > t = 0:1/Fs:(3600*3); > % signal > x = cos(2*pi*(1/5.5)*t)+randn(size(t)); > % compute a Welch estimate > % Hamming window, segment length of 500 samples Would it be wise to use a Window in this situation, since the data is not directly measurement data, but more data extracted from a FEA. Windowing with Hamming would lead to neglection of the side lobs of the sample. > hspec = spectrum.welch('Hamming',500); > PsdEst = psd(hspec,x,'Fs',2); > plot(PsdEst); > % compute a periodogram > hp = spectrum.periodogram; > Psdper = psd(hp,x,'Fs',2); > figure; > plot(Psdper); > > > I've used a segment length of 500. If you reduce this length, the variance of the PSD estimate will reduce further, but the width of the peak will broaden. Note that the Welch estimate is much less "noisy" than the periodogram, but this noise reduction comes with a loss of resolution. You're right in this one where the Welch estimate contains less noise than the PSD and increasing the number of sample will create a more wispy (don't know if it is the correct word for this, meaning: more peaks) result. > > Another way would be to reduce the noise in the estimate is to first low pass filter your input. If your frequency of interest is 1/5.5 hertz, you can remove everything above a certain frequency, let's say 0.25. > > hfilt = fdesign.lowpass('Fp,Fst,Ap,Ast',0.25,0.3,0.1,50,2); > d = design(hfilt); > % filter the input data > y = filter(d,x); > % if you need zero-phased filtering > y1 = filtfilt(d.Numerator,1,x); > > Hope that helps, > Wayne
From: Wayne King on 27 Jul 2010 10:24 "Jos " <j.g.h.haarman(a)student.utwente.nl> wrote in message <i2moqm$24j$1(a)fred.mathworks.com>... > "Wayne King" <wmkingty(a)gmail.com> wrote in message <i2mcmr$k20$1(a)fred.mathworks.com>... > > "Jos " <j.g.h.haarman(a)student.utwente.nl> wrote in message <i2majc$462$1(a)fred.mathworks.com>... > > > Hi all, > > > > > > Currently I'm working on a project that involves an analysis in a FE package which calculates the loads on a part of our structure due to dynamic responses. This is stored as a load-time history file. I want to investigate the FFT and the PSD of this signal to eventually see if the signal contains any non-predicted frequencies and to check if the signal is narrow- or widebanded. (any input on the last subject is highly appreciated) > > > > > > The following parameter can be provided. > > > Sample time is 0.5 sec. > > > Total time history is 3hours > > > Period of the inputsignal (sea state) is 5.5 sec. *Probably varied in the FEA program* > > > > > > I store the data in the vector 'sum_hor'. > > > > > > After investigating the files in the help and online I came to the following m-file. The question to you guys is if I'm making any mistakes (which I assume I do).. > > > > > > ------------------ > > > > > > H = sum_hor; > > > nfft = 2^nextpow2(length(H)); > > > Pxx = abs(fft(H,nfft)).^2 / length(H); > > > > > > %FFT > > > figure(1) > > > plot(Pxx(1:length(Pxx)/2) > > > > > > Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2)); > > > > > > %PSD > > > figure(2) > > > plot(Hpsd) > > > > > > ------------------ > > > > > > In addition, the results generated by this m-file are a bit noisy, any suggestions how to cope with this problem? > > > > > > Thanks in advance for all the help, > > > Jos > > > > Hi Jos, The easiest way to arrive at a PSD estimate is using spectral analysis objects. > > > > Assumming x is your signal. The following will compute a Welch spectral estimate and periodogram (the periodogram is what you've done above). You will note that the Welch estimate exhibits markedly less variance than the periodogram. I'm assuming that you use the entire 3 hour recording. Of course, it may be beneficial to do a time-frequency analysis on that data if the statistical properties change over time. For this example, I'll assume the data is wide-sense stationary and use the entire data record. > > > > Fs = 2; % your sampling interval is 0.5 > > t = 0:1/Fs:(3600*3); > > % signal > > x = cos(2*pi*(1/5.5)*t)+randn(size(t)); > > % compute a Welch estimate > > % Hamming window, segment length of 500 samples > > Would it be wise to use a Window in this situation, since the data is not directly measurement data, but more data extracted from a FEA. Windowing with Hamming would lead to neglection of the side lobs of the sample. > > > hspec = spectrum.welch('Hamming',500); > > PsdEst = psd(hspec,x,'Fs',2); > > plot(PsdEst); > > % compute a periodogram > > hp = spectrum.periodogram; > > Psdper = psd(hp,x,'Fs',2); > > figure; > > plot(Psdper); > > > > > > I've used a segment length of 500. If you reduce this length, the variance of the PSD estimate will reduce further, but the width of the peak will broaden. Note that the Welch estimate is much less "noisy" than the periodogram, but this noise reduction comes with a loss of resolution. > > You're right in this one where the Welch estimate contains less noise than the PSD and increasing the number of sample will create a more wispy (don't know if it is the correct word for this, meaning: more peaks) result. > > > > > Another way would be to reduce the noise in the estimate is to first low pass filter your input. If your frequency of interest is 1/5.5 hertz, you can remove everything above a certain frequency, let's say 0.25. > > > > hfilt = fdesign.lowpass('Fp,Fst,Ap,Ast',0.25,0.3,0.1,50,2); > > d = design(hfilt); > > % filter the input data > > y = filter(d,x); > > % if you need zero-phased filtering > > y1 = filtfilt(d.Numerator,1,x); > > > > Hope that helps, > > Wayne Hi Jos, the Welch estimate used above multiplies each segment of length 500 by a Hamming window and overlaps the segments by 50%. You can modify the OverlapPercent property of the hspec if you wish this overlap percentage to be different. If you wish to just window the entire data vector, you can do hper = spectrum.periodogram('Hamming'); % Assuming x is your data vector PsdEst = psd(hper,x,'Fs',2); But I think you'll find that the Welch estimate is much smoother because you are averaging over segments. Wayne
From: Jos on 27 Jul 2010 10:36
"Wayne King" <wmkingty(a)gmail.com> wrote in message <i2mq67$372$1(a)fred.mathworks.com>... > "Jos " <j.g.h.haarman(a)student.utwente.nl> wrote in message <i2moqm$24j$1(a)fred.mathworks.com>... > > "Wayne King" <wmkingty(a)gmail.com> wrote in message <i2mcmr$k20$1(a)fred.mathworks.com>... > > > "Jos " <j.g.h.haarman(a)student.utwente.nl> wrote in message <i2majc$462$1(a)fred.mathworks.com>... > > > > Hi all, > > > > > > > > Currently I'm working on a project that involves an analysis in a FE package which calculates the loads on a part of our structure due to dynamic responses. This is stored as a load-time history file. I want to investigate the FFT and the PSD of this signal to eventually see if the signal contains any non-predicted frequencies and to check if the signal is narrow- or widebanded. (any input on the last subject is highly appreciated) > > > > > > > > The following parameter can be provided. > > > > Sample time is 0.5 sec. > > > > Total time history is 3hours > > > > Period of the inputsignal (sea state) is 5.5 sec. *Probably varied in the FEA program* > > > > > > > > I store the data in the vector 'sum_hor'. > > > > > > > > After investigating the files in the help and online I came to the following m-file. The question to you guys is if I'm making any mistakes (which I assume I do).. > > > > > > > > ------------------ > > > > > > > > H = sum_hor; > > > > nfft = 2^nextpow2(length(H)); > > > > Pxx = abs(fft(H,nfft)).^2 / length(H); > > > > > > > > %FFT > > > > figure(1) > > > > plot(Pxx(1:length(Pxx)/2) > > > > > > > > Hpsd = dspdata.psd(Pxx(1:length(Pxx)/2)); > > > > > > > > %PSD > > > > figure(2) > > > > plot(Hpsd) > > > > > > > > ------------------ > > > > > > > > In addition, the results generated by this m-file are a bit noisy, any suggestions how to cope with this problem? > > > > > > > > Thanks in advance for all the help, > > > > Jos > > > > > > Hi Jos, The easiest way to arrive at a PSD estimate is using spectral analysis objects. > > > > > > Assumming x is your signal. The following will compute a Welch spectral estimate and periodogram (the periodogram is what you've done above). You will note that the Welch estimate exhibits markedly less variance than the periodogram. I'm assuming that you use the entire 3 hour recording. Of course, it may be beneficial to do a time-frequency analysis on that data if the statistical properties change over time. For this example, I'll assume the data is wide-sense stationary and use the entire data record. > > > > > > Fs = 2; % your sampling interval is 0.5 > > > t = 0:1/Fs:(3600*3); > > > % signal > > > x = cos(2*pi*(1/5.5)*t)+randn(size(t)); > > > % compute a Welch estimate > > > % Hamming window, segment length of 500 samples > > > > Would it be wise to use a Window in this situation, since the data is not directly measurement data, but more data extracted from a FEA. Windowing with Hamming would lead to neglection of the side lobs of the sample. > > > > > hspec = spectrum.welch('Hamming',500); > > > PsdEst = psd(hspec,x,'Fs',2); > > > plot(PsdEst); > > > % compute a periodogram > > > hp = spectrum.periodogram; > > > Psdper = psd(hp,x,'Fs',2); > > > figure; > > > plot(Psdper); > > > > > > > > > I've used a segment length of 500. If you reduce this length, the variance of the PSD estimate will reduce further, but the width of the peak will broaden. Note that the Welch estimate is much less "noisy" than the periodogram, but this noise reduction comes with a loss of resolution. > > > > You're right in this one where the Welch estimate contains less noise than the PSD and increasing the number of sample will create a more wispy (don't know if it is the correct word for this, meaning: more peaks) result. > > > > > > > > Another way would be to reduce the noise in the estimate is to first low pass filter your input. If your frequency of interest is 1/5.5 hertz, you can remove everything above a certain frequency, let's say 0.25. > > > > > > hfilt = fdesign.lowpass('Fp,Fst,Ap,Ast',0.25,0.3,0.1,50,2); > > > d = design(hfilt); > > > % filter the input data > > > y = filter(d,x); > > > % if you need zero-phased filtering > > > y1 = filtfilt(d.Numerator,1,x); > > > > > > Hope that helps, > > > Wayne > > Hi Jos, the Welch estimate used above multiplies each segment of length 500 by a Hamming window and overlaps the segments by 50%. You can modify the OverlapPercent property of the hspec if you wish this overlap percentage to be different. > > If you wish to just window the entire data vector, you can do > > hper = spectrum.periodogram('Hamming'); > % Assuming x is your data vector > PsdEst = psd(hper,x,'Fs',2); > > But I think you'll find that the Welch estimate is much smoother because you are averaging over segments. > > Wayne Definitly true what you stated on the Welch estimate, the windowing over the whole set of data gives a much more noisy outcome, which is not very surprisingly, when you think about averaging over a whole window or 500 windows and trying to reduce the noice. Thank for your fast replies, I think I'm saved for the moment! |