Why is signal amplitude so low after applying FFT?

55 次查看（过去 30 天）
M Prabhu
M Prabhu 2021-10-2
编辑： dpb 2021-10-3
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)')
  4 个评论
显示 2更早的评论
dpb
dpb 2021-10-2
编辑：dpb 2021-10-2
Fs is (and must be) identically the sampling rate to compute the correct frequency vector.
If your data are sampled at 40 Hz, then that's what it is for these data and the Nyquist frequency will be 20 Hz.
At 40 Hz, a 5 Hz signal is quite oversampled; if the signal is noisy you'll need to average.

请先登录，再进行评论。

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

回答（2 个）

G A
G A 2021-10-2
dpb
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...
  6 个评论
显示 4更早的评论
M Prabhu
M Prabhu 2021-10-3
Thank you. But I don't think I understand you completely. Let me try reading up a bit and then get back.
dpb
dpb 2021-10-3
编辑：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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Signal Generation and Preprocessing 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by