Why is signal amplitude so low after applying FFT?

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]

Attach the data file; we can't tell anything about the signal just from the code.

There's an example at doc fft of a one-sided PSD that is normalized to match the input amplitude of the time signal to illustrate the process...it also shows the frequency calculation (I didn't try to read the code above; it may be ok, too, but you can compare your result with that from the example; they should be the same).

NB: the RMS of the time signal is that of all energies contained in the signal at any instant in time -- the PSD shows the energy in a given frequency bin -- it's quite possible if it is a broadband signal that the energy is spread out over the full frequency range and so the total energy to match the time signal has to be the sum/integral over the frequency range.

Also, unless the sampling frequency is high enough and matches the actual frequencies in the signal, the PSD peaks will not be identically at the frequency of the FFT and so the total energy at a given frequency will have to also be the integral of the peak, not just the single bin point magnitude.

[M Prabhu]

Thank you MVP. An extract of the data file is attached. It contains 20000 records in two columns (the original datafile contains about 200,000 records). The first column is time and second column contains the signal of interest. When I use a small segment from the data file, the amplitudes are larger, but if I use a large segment or the entire data, it shows smaller amplitudes.But then neither of them seem to match the average amplitude seen visually when I plot a time history of the data. They seem much smaller.

I will go through the references you have indicated.

If I change fs (in code) there is a difference in output (amplitude and frequency changes). I do not know what value of fs should I select. Is there any guideline on this? Is it linked to the sampling rate (which is 40Hz in my case)? Is it linked to the expected largest frequency of the signal? That would be < 5 Hz in my case, all higher frequencies are noise. What would you recommend for my data?

[dpb]

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.

[M Prabhu]

Thank you! I've amended my code accordingly.

Answer 1

[G A]

https://www.mathworks.com/matlabcentral/answers/162846-amplitude-of-signal-after-fft-operation

https://www.sjsu.edu/people/burford.furman/docs/me120/FFT_tutorial_NI.pdf

Answer 2

[dpb]

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]

Actually, the above might be improved some more if one were to apply detrend to the result of the reshape() operation to detrend each segment individually. This would account for at least some of the nonstationarity...

...
NAvg=16;                                        % use 16 averages    
L=L/NAvg;                                       % adjust spectrum length to match
y=detrend(reshape(y,numel(y)/NAvg,[]));         % reshape y to NAvg columns and detrend
Y=fft(y);
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])    

[M Prabhu]

Thanks a ton for the detailed explantion and code. Was very helpful. There is clearly a significant improvement in the output after the above corrections. Looks a lot cleaner. The amplitudes have also gone up marginally now, but I still don't have what I felt I should be getting (my original question). I expected amplitudes of at least 2 mm by looking at the time history plot, whereas P1 is still showing 0.2 to 0.3 mm. By the way, my input signal is in mm. I understand that, after FFT, the output unit remains the same as input unit.

[M Prabhu]

By the way, when I generate a periodogram for this signal. It is also very noisy. I tried detrending and averaging the signal before running the above code, but it did not help. Further, when I detrend, many curves are superimposed over each other, as can be seen in the attached plots. This is not seen if I only average. Any way to remove noise from periodograms?

y=x;                        % x is the signal vector
L=numel(x_time);            % x_time is the time vector
Navg=12;                    % 12 sample averages
L=L/Navg;                   % adjust spectrum length to match
y=detrend(reshape(y,numel(y)/Navg,[])); % Reduce DC component & noise with detrend and averaging. But results in many overlapping PSD curves
% y=reshape(y,numel(y)/Navg,[]); % This does not result in multiple output PSD curves, if used instead of previous line 
figure;
periodogram(y,rectwin(length(y)),length(y),Fs)%
title('original signal');
axis([0 20 -80 40]);

[dpb]

The rms of the time history will only match the one frequency PSD amplitude when the time signal is a pure tone at the exact frequency of one of the bins of the FFT (or the integral of the single peak).

Look again at the example in the FFT doc --

>> Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;
>> S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
>> X = S + 2*randn(size(t));
>> rms(X)
ans =
          2.17
>>

but the peaks are 0.7 and 1.0 -- and, comparatively this is a much less noisy signal than yours and also had all the energy comprised in the two sinusoids whereas your signal is very broadband, not narrow at all.

And, again, you have to integrate the peak that you do get in order to have the total energy -- not just the peak bin amplitude.

Averaging only works to remove the random uncorrelated noise when done in the frequency domain, it does NOT help in the time domain.

[M Prabhu]

Thank you. But I don't think I understand you completely. Let me try reading up a bit and then get back.

[dpb]

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.