How to obtain correct frequency and amplitude via FFT

Hi

I have tried to test a Fast Fourier Transform using a simple example signal. The signal has a length of 2 seconds, a frequency of f = 2 and an amplitude of a = 3. However, when I apply the FFT and determine the maxima via the command findpeaks, I get f = 1.8311 and a = 2.5763. This means that the results differ quite a lot from my input.

I have already looked around here in the forum, but so far I have not been able to find a real solution to my problem. I hope that someone here has an idea what I am doing wrong.

I'd appricate the help! Thanks

clear all; close all; clc;

% create sinusoidal signal

f_sign = 2          ; % signal frequency
a_sign = 3          ; % signal amplitude
t_sign = 2          ; % length of signal [s]
dt_sign = 0.0001    ; % sample rate [s]

nt_sign = t_sign/dt_sign        ; % discrete number of datapoints

t = 0:dt_sign:t_sign            ; % Time vector

y = a_sign*sin(2*pi*f_sign*t)   ; % create signal


% Fourier transformation

y_fft = y - mean(y);

nfft = 2^nextpow2(nt_sign);
FFT = fft(y_fft,nfft)/nt_sign;
FFT = 2*abs(FFT(1:nfft/2+1));
fs=1/dt_sign;
f = fs/2*linspace(0,1,nfft/2+1);

% find maximum amplitude and frequency
[pks,locs] = findpeaks(FFT,f);
[pks_max,idx_pks_max] = max(pks);
locs_max = locs(idx_pks_max);

% results
f_fft = locs_max

f_fft = 1.8311

a_fft = pks_max

a_fft = 2.5763

[Paul]

Are you trying to construct an example where you get back the expected result knowing that the amplitude is 3 and the frequency is 2 Hz?

Or are you trying to work an example for a signal as if you do not know the amplitude and frequency and playing around with fft parameters to see how they change the results of the analysis?

[Nicolas]

I have a signal whose amplitude and frequency I do not know. Since I am using the FFT for the first time in Matlab, I first wanted to test the result on a signal whose values I know to see if it works. Hence the very simple example above.

[Paul]

OK. Then I think you have what you need. Just one more thing to be aware of. Consider these lines of code:

nfft = 2^nextpow2(nt_sign);
FFT = fft(y_fft,nfft)/nt_sign;
FFT = 2*abs(FFT(1:nfft/2+1));
fs=1/dt_sign;
f = fs/2*linspace(0,1,nfft/2+1);

Here, f(1) = 0 and f(end) = fs/2 (nfft/2 is even, the computation of f would be different if nfft is odd). These frequencies only show up once in the full fft, so it doesn't make sense to multiply the corresponding values in FFT by 2 when only looking at the first half of the spectrum. The code subtracts the mean from y before taking the fft, so FFT(1) should be zero anyway, but that's not the case for FFT(nfft/2+1).

[Nicolas]

But the multiplication by 2 seems to be important for the correct amplitude, doesn't it? Without it, I would get an amplitude of 1.5 in the example, which is only half of the actual amplitude.

[Star Strider]

The multiplication-by-2 is important, and is in the documentation. The reason is that the signal energy is divided (in the two-sided fft) between the ‘positive’ and ‘negative’ frequencies. The multiplication-by-2 (in either the one-sided to two-sided fft corrects) for this, giving an appropriate amplitude.

[Paul]

Focusing on amplitude, because that's what @Nicolas seems to be interest is ...

I thought I was clear that my comment was only relevant to the FFT(1) and FFT(N/2+1), which are the elements of FFT that correspond to f(1) and f(end). So, I'm only talking about those two elements.

Suppose we have "sine wave" at zero frequency with amplitude = 3. For simplicity, I'll assume N = 8.

N = 8;
n = (0:(N-1));
x = 3*cos(0*n)

x = 1×8

3 3 3 3 3 3 3 3

X = fft(x)/N

X = 1×8

3 0 0 0 0 0 0 0

As we can see, the first point in X is 3, which is the amplitude of the input. So why multiply X(1) by 2?

Now, let the frequency be pi, i.e., the Nyquist frequency

x = 3*cos(pi*n)

x = 1×8

3 -3 3 -3 3 -3 3 -3

X = fft(x)/N

X = 1×8

0 0 0 0 3 0 0 0

The element X(N/2+1) that corresponds to the Nyquist frequency is 3, which is the amplitude of the input. Why multiply

X(N/2+1) by 2?

Keep in mind that if N is odd, no element of X corresponds to the Nyquist frequency. But we still have the same issue with X(1) regardless.

Check out this doc page link and you'll see these lines

P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

i.e., the first element of P1 is not multiplied by 2, and neither is the last element for L even.

Answer 1

[Star Strider]

Increase the frequency resolution (I doubled it here, although you can increase it further) with

nfft = 2^nextpow2(nt_sign)*2

and the results significantly improve.

With that:

f_fft = 1.9836

a_fft = 3.0070

and you can likely get improved accuracy with finer ffrequency resolution.

% clear all; close all; clc;

% create sinusoidal signal

f_sign = 2          ; % signal frequency
a_sign = 3          ; % signal amplitude
t_sign = 2          ; % length of signal [s]
dt_sign = 0.0001    ; % sample rate [s]

nt_sign = t_sign/dt_sign        ; % discrete number of datapoints

t = 0:dt_sign:t_sign            ; % Time vector

y = a_sign*sin(2*pi*f_sign*t)   ; % create signal


% Fourier transformation

y_fft = y - mean(y);

nfft = 2^nextpow2(nt_sign)*2

nfft = 65536

FFT = fft(y_fft,nfft)/nt_sign;
FFT = 2*abs(FFT(1:nfft/2+1));
fs=1/dt_sign;
f = fs/2*linspace(0,1,nfft/2+1);

% find maximum amplitude and frequency
[pks,locs] = findpeaks(FFT,f);
[pks_max,idx_pks_max] = max(pks);
locs_max = locs(idx_pks_max);

% results
f_fft = locs_max

f_fft = 1.9836

a_fft = pks_max

a_fft = 3.0070

.

[Nicolas]

I really thought I had tried that, but obviously I haven't. Thank you very much!

[Star Strider]

As always, my pleasure!