Spectral structure of sinusoidal model

Question

Let us consider the following code:

function [ x ] = generate1(N,m,A3)
f1 = 100;
f2 = 200;
T = 1./f1;
t = (0:(N*T/m):(N*T))'; %'
wn = randn(length(t),1); %zero mean variance 1
x = 20.*sin(2.*pi.*f1.*t) + 30.*cos(2.*pi.*f2.*t) + A3.*wn;
%[pks,locs] = findpeaks(x);
 %plot(x);
end

I know peaks in the Fourier domain are represented at these frequencies, which are present in the signal.

For example, let us take the plot of the Fourier transform of this signal:

y=generate1(3,500,1);

and plot

plot(abs(fft(y)))

but clearly it does not show peaks at the frequencies given in the signal. What is the problem?

Answer 1

Lets analyze the computation. Your function takes $N$ periods of the $100Hz$ sine wave and correspondingly $2N$ periods of the $200Hz$ sine wave and samples them in $m$ equidistant points.

Then you perform an FFT of size $m$. If we forget for the moment that there is a factor of 100Hz used in the production, we have sampled waves of $N$ and $2N$ full oscillations, and the FFT detects or handles sinusoids with up to $m/2$ full oscillations.

In the plot from of the FFT result scaled from 0 to m one would expect the first wave at position N, the second at position 2N, and the mirror frequencies at m-1-N and m-1-2N.

In your example with N=3 and m=500, the 100Hz component is at position 3 and the 200Hz component at position 6.

To get more striking results, either reduce the range in the display to $[0:10N]$ or similar, or choose the magnitudes of m and N closer together, m=8N to m=20N should work well.

---

Added matlab interna: the fft command produces the DC at position 1, the first array position. The FFT with sampling frequency $f_s$ and $m$ samples produces a sampling of the frequency spectrum ranging from $-f_s/2$ to $f_s/2$ in $m$ steps, assuming $m$ is even.

The fft command arranges these frequencies by wrapping around at $fs/2$, so that a component for $f>fs/2$ in reality is a component for $f-f_s<0$. For real signals, the negative spectrum is essentially the same as the positive, complex conjugation is the only difference. So one can just disregard the second half of the result of fft for visualization purposes.

The generate1 function should at least also return the sampling frequency fs=(f1*m)/N. Then

m=500;
y,fs=generate1(3,m,1);
A=abs(fft(y))/m;
f=(0:1:m-1)*fs;
plot(f(1:m/2),A(1:m/2))

produces the full spectrum with correct scale. Replace m/2 in the plot command with a smaller number to display only parts of the spectrum.