Plotting Real Amplitude of Sine Wave in FFT

Question

Consider a sine wave having an amplitude of 9 and a frequency of 0.3 Hz. I wish to perform an fft and plot the frequency on the x-axis and the real amplitude (i.e. 9) on the y-axis. I am using the following MATLAB code:

Fs = 10;
data = xlsread('Signal.xlsx');
Y(1)=0;
Y = fft(data);
N = length(Y);
freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = abs(Y(1:floor(N/2)))/floor(N/2);
plot(freq,amplitude);
xlabel('Frequency [Hz]');
ylabel('Amplitude');

When running the code, an amplitude close to 9 is obtained at a frequency of 0.3 Hz. In reality, I have a sinusoidal signal of a rotating motor shaft. I am measuring the angle using an encoder. The data can be found at:

http://www.dropbox.com/s/pe3g44b48ugauqp/Signal.xlsx?dl=0&m=

I know that the amplitude of the sine wave should be about 0.6 m but when conducting an fft, it is indicating that the amplitude is around 0.4 m. I am sampling at 10Hz and have about 185 000 data points. When trying to reduce the data, the amplitude of the sine wave is approximating the desired one i.e. 0.6 m. As far as I know, the more data points you have, the better will the resolution be when conducting an fft. So why am I getting the opposite? Can someone please help me? Thanks.

Answer 1

Judging from your post, you have over 5 hours of data, your are using to calculate the spectrum. This is a lot of data points to process with FFT. Let's focus on your code and modify it a little bit to generate the $0.3\; \mathtt{Hz}$ sinusoid with an amplitude of $9$:

Fs = 10;
t = 0:0.1:1000;
A = 9;
f0 = 0.3;
data = A*sin(2*pi*f0*t);
N = length(data);    

Y(1)=0;
Y = fft(data);

freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = abs(Y(1:floor(N/2)))/floor(N/2);

subplot(211)
plot(t, data)
xlabel('Time [s]')
ylabel('Amplitude')
title('Time domain signal')
grid on

subplot(212)
plot(freq, amplitude);
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Amplitude spectrum')
grid on

This gives us what you expect to obtain.

Now let's assume that your $0.3\; \mathtt{Hz}$ signal is not always present (sometimes there is only noise present). I suspect this is how your data looks like, especially over so many hours. I will approximate this situation by adding the random noise after half of a sinusoid duration: data(5000:10000)=randn(1,5001); (consider it as someone turned off the machine halfway through the experiment). Obviously you can make it even more intelligent by having signal for 1 second, followed by 1 second of noise, and so on.

How it will affect your amplitude spectrum? Obviously you won't get $9$ as an amplitude because there is only half of the energy present before at $0.3\; \mathtt{Hz}$. Thus you will observe something like:

If you want to detect a true amplitude, then I suggest you to use the Time-Frequency analysis tools. Easiest one is the Short-Time Fourier Transform. You will be able to observe the amplitude of your component at a given time, and extract it later on.

---

After analysing your signal indeed you can expect the peak with an amplitude of $0.62$ at $0.2\; \mathtt{Hz}$ (not as you claiming $0.3\; \mathtt{Hz}$):

You problem relates to the leakage as you get lot's of frequency bins and likelihood of your frequency falling exactly between get's lower. Indeed, to many is also not a good idea. Below you can observe two spectra for two synthesized sinusoids with same amount of samples as your dataset and amplitude $0.62$. Frequencies are $f=0.2\; \mathtt{Hz}$, and $f=0.2003\; \mathtt{Hz}$. You can observe that for a perfect match, amplitude is estimated correctly, although for $0.2003\; \mathtt{Hz}$ you get $0.41$.

One workaround to deal with that problem in your case, would be to apply windowing to your signal. Although it will decrease amount of leakage, but still it won't remove it absolutely. Below you can observe example for the Hamming window.

In my opinion better solution would be to use less samples. For example if you expect your frequency to be around $0.2\; \mathtt{Hz}$, then you can choose $100$ samples (I know - very crude), which will give you resolution of $0.1\; \mathtt{Hz}$ and component will fall exactly inside. Here is the example:

One last remark, regarding your code. I do not understand what your declaration Y(1)=0 is suppose to do? If it's initialization, then it is unnecessary. If you are trying to remove the DC component, then it's not a good place nor way.

---

In case of using other window than rectangular one, i.e. Hamming you must normalize it properly. In case of DFT you must divide your output samples by sum of all window samples. In case of 'no windowing' which is in fact rectangular window you get sum of all samples being equal to N. Putting it into code:

win = hamming(N); % rectwin(N), hann(N), ...
data = data .* win;
Y = fft(data)/sum(win);

For more theory please refer to other posts on DSP SE or take a look into this great paper.

Answer 2

When doing an FFT of a length on the order of 100k samples, the frequency resolution will be on the order of 0.001%, which means even tiny variations in the signal frequency (from being perfectly integer periodic in the FFT length) or the sampling frequency can spread the energy of a sinusoid out into more than one FFT result bin, thus reducing the magnitude of the peak bin.