조회 수: 42(최근 30일)

표시 이전 댓글

I have a signal measured @40Hz. I want to compute the average value of signal vs frequency using FFT and also the PSD of the signal.

When I ran the following code, the amplitude that I see in the plot is very small compared to the input amplitudes. When the number of data points is lesser (say 1000), it seems to be better (although I still think it may not be correct) [see Fig-1]. But when the data is large (say 20,000), it is clearly incorrect [see Fig-2]. I zoomed into the large sample data to verify if the average amplitude of signal was indeed so low, but it is not so (see Fig-4). The minimum signal amplitude is typically more than 1 mm, with the average probably around 3-4 mm visually. But the plot shows 0.25mm at best.

Where am I going wrong? (in choice of fs?) I am also not sure if the computed frequency is correct. How can I be sure?

I have attached the data file. The output seems quite sensitive to the chosen fs. Is there any guideline to choose fs? My signal is sampled at 40 Hz, and the frequency of interest is less than 5 Hz. All higher frequencies in the signal are noise and I am not intersted in them. What value of fs should I choose? And is the frequency (freq) that I am calculating or getting after PSD in Hz or in rad/s? SInce I have not used 2*pi in any the calculation, it should be rad/s, but then it doesn't match my expectation based on the time history plot of the signal. I have assumed it is in Hz (but I could be wrong).

load 'signal.mat' % this contains the signal x (sampled at 40 Hz) used for the current analysis

fs=100;

tm=0:1/fs:(length(x)-1)/fs;

figure;

plot(tm,x); % time history plot of inputs signal

xlabel('time (s)')

ylabel('Input Signal (mm)')

figure;

X=fft(x,length(x))*(2/length(x)); % fft of input signal. Output X is a complex number

mag_X = abs(X); % magnitude of complex number

freq=0:1/tm(end):fs/2-(1/tm(end)); % freq vector

plot(freq,mag_X(1:length(freq)));

xlabel('Frequency (Hz)')

ylabel('Amplitude (mm)')

figure;

psdest = psd(spectrum.periodogram,x,'fs',fs,'NFFT',length(x)); % PSD of signal

plot(psdest.Frequencies,psdest.Data);

xlabel('Frequency (Hz)')

ylabel('PSD (mm^2/Hz)')

dpb
2021년 10월 2일

편집: dpb
2021년 10월 2일

Several issues here -- the first being the RMS value masks the signal amplitude with the DC value unless you first detrend the time signal. And, it shows some nonstationarity on out in the time trace but I've made no compensations for anything but the mean. You can look at the doc for it and see other options available or "roll your own".

Attached are three figures of a single-sided PSD from the FFT -- they are first without removing the mean; the overall amplited scale is large so on a linear y axis there's nothing of note at all in the plot -- the attached is semilogy so one does see the peak around 1 Hz at at roughly 0.6 or so. It also shows the PSD estimate is quite noisy.

The code that did the above is

v=readmatrix('data.csv'); % read the data file

t=v(1:end-1,1);y=v(1:end-1,2); % create convenient variables; round to even number points

Fs=40; % sampling frequency, 40 Hz

L=numel(t); % length of time signal

P2=abs(fft(y)/L); % two-sided PSD from FFT normalized to signal

P1=P2(1:L/2+1); % one-side PSD elements

P1=2*P1,P1(1,end)=P1(1,end)/2; % normalize to total power (DC, Fmax are already)

f=Fs*(0:(L/2))/L; % compute the frequency vector for plotting

plot(f,P1)

set(gca,'yscale','log')

xlim([0 10])

Repeat the above except detrend the time series first...

...

y=detrend(y); % remove DC from time trace

P2=abs(fft(y)/L); % two-sided PSD from FFT normalized to signal

...

This produced the following figure (on linear y axis instead of log). Now you can see the spectrum; the DC component is still sizable because of the nonstationarity of the time trace but it's enough smaller that it doesn't swamp everything else. Of course, one can just set it to zero arbitrarily afterwards.

Lastly, because it is such a noisy-looking signal, let's try averaging and see if it really is just uncorrelated noise -- again we start with detrended y and add

...

NAvg=16; % use 16 averages

L=L/NAvg; % adjust spectrum length to match

Y=fft(reshape(y,numel(y)/NAvg,[])); % reshape y to NAvg columns

P2=mean(abs(Y/L),2); % take FFT and average the columns

P1=P2(1:L/2+1); % rest is the same from 2-sided to 1-sided

P1=2*P1;P1(1,end)=P1(1,end)/2;

f=Fs*(0:(L/2))/L; % have to recompute f vector for new L

plot(f,P1)

xlim([0 10])

Now we have quite a lot nicer result...

dpb
2021년 10월 3일

Experiment with the example above -- for starters change the frequency by a little bit in one of the sinusoids and observe what happens when the binning isn't an exact multiple of the frequency of the sinusoid.

For specific example, there

>> df=Fs/L; % the FFT frequency binning delta-f

>> 50/df % which bin is nearest to the 50 Hz signal?

ans =

75.00

>>

Turns out they chose numbers that makes it exactly hit a bin -- how convenient!

Now try

>> 51/df % 51 Hz instead is halfway between two bins

ans =

76.50

>>

So regenerate the time series in the example changing the 50 Hz signal to 51 Hz and follow through the rest of the example and see what the peak amplitudes are then...

The amplitude isn't exactly 0.7 even for the binning-matching case owing to the noise; as the example also shows, if it were noise free you would, in fact, retrieve identically 0.7 and 1.0 for the two peaks.

The key thing to note is that in the 51 Hz case the energy is spread between the two nearest frequency bins since there isn't a bin that matches precisely 51 Hz. Ergo, to get an estimate of the total energy in the peak you have to integrate the peak area; you can't just take the one peak value.

That's the simple case where you know precisely what the input signal is; you've got the case of sampled (very) noisy data with a system that clearly has the response frequency spread out around the expected center frequency; it isn't at all a pure sine wave going in at that frequency but a smear of energy roughly centered around it.

The time history contains all energies at a given instant of time in its amplitude; the frequency spectrum contains only the energy of the specific point frequency over all time; if all the energy isn't at one frequency, not all the energy is going to show up in just the one point value; you've got to add up all the bins that constitute the area of interest. With noise, even that isn't going to be perfect although averaging over longer time series will help IF the system is stationary -- and your time history shows some indications that that isn't a perfect assumption, either.

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

Start Hunting!