spectral structure of sinusoidal model

Question

let us consider 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

as i know peaks in Fourier domain represent at this frequency,which are present in signal,for example let us take plot of Fourier transform of this signal let us run this signal

y=generate1(3,500,1);

and plot

plot(abs(fft(y)))

but clearly it does not shows me peaks at frequency given in signal,what is problem?please help me

Answer 1

Well, first off, you haven't given plot any specific values for the x-axis, so it just uses the array indices – why should these correspond to frequencies? The frequencies corresponding to Fourier coefficients can be computed like this:

f = (0 : numel(y) - 1) / numel(y) * fs;

where fs is the sampling frequency of your signal y. You don't specify the sampling frequency, but from looking at your code it is clear it is the inverse of N*T/m. Extending your function to return this information:

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

Putting this together,

[y, fs] = generate1(3,500,1);
f = (0 : numel(y) - 1) / numel(y) * fs;
plot(f, abs(fft(y)))

you get

If you zoom in, you will see that you have two peaks at the correct frequencies, 100 and 200. But you also have two peaks at 16467 and 16567.

That's because the Fourier transform of a real-valued signal is always symmetric. Usually this symmetry is stated to be around 0, which you get with a different interpretation of the frequencies attached to the Fourier coefficients:

f = (0 : numel(y) - 1) / numel(y) * fs;
f(f > fs/2) = f(f > fs/2) - fs;

These two interpretations are equivalent with respect to the reconstructed signal. The result:

Since the frequency vector f is no longer increasing, it is better to use fftshift for plotting:

plot(fftshift(f), fftshift(abs(fft(y))))

Zoomed in view on the result:

Please note that a simple Fourier transform is in general not a good way to assess the spectral content of a signal, for two reasons:

1. The discrete Fourier transform treats the signal effectively as periodic, i.e. it assumes the signal continues after the end by repeating the beginning. This assumption will in many cases introduce a discontinuity, which distorts the spectral content. This problem can be solved e.g. by using a taper or window function.

2. Especially for noisy signals, the absolute square of Fourier coefficients is a very imprecise estimate of the signal's power. It is better to cut the signal into shorter windows, and average the different estimates from the different windows. See Welch's method, or in Matlab, pwelch.