# Compute the power spectrum using FFT method

106 views (last 30 days)
An tran on 31 Jan 2015
Edited: Walter Roberson on 9 Apr 2020
I have a project as follows: there are 2 sinusoids in the white noise background. 32 received samples are u(n)=exp(j2pif1n)+exp(j2pif2n+phase)+w(n), n=0,1,2..31 where phase is a random phase and w(n) is the white noise. f1=0.115 and f2=0.135, signal to noise ration is 20dB. Compute the power spectrum by FFT method.
I am just the beginner to Matlab programming. Could you please tell me what I went wrong in my code below so that I can learn from my mistakes. Any recommendations to enhance my codes? I would greatly appreciate it
clear all; format long; f1=0.115; f2=0.135; n=0:31;
phase=2*pi*rand();
w=0.1*randn();
u=exp(j*2*pi*f1*n)+exp(j*2*pi*f2*n+phase)+w;
%Using FFT
U=fft(u);
fvals=(0:length(U)-1)/length(U);
mag=abs(U);
subplot(1,2,1);plot(fvals,mag.^2);
xlabel('Frenquency'); ylabel('P(f)'); title('Power spectrum by FFT method');
subplot(1,2,2);plot(fvals,10*log10(mag.^2));
xlabel('Frequency'); ylabel('P(f)'); title('Power spectrum by FFT method in dB');

Youssef Khmou on 1 Feb 2015
Edited: Youssef Khmou on 1 Feb 2015
The essay is correct, Some remarks are :
1. Try to give variables names different than those of already built in functions , like phase.
2.The random variables is scalar, instead use randn(size(n)), to generate a vector.
3.the function abs(U) computes the magnitude if U is complex, no need to square the result.
4.the defined SNR is correct only for real values signal, which is not in this case, for this example use real(sinusoid) instead.
5. You have to scale the fft result in order to obtain the correct power ( std(u)).
6. You can increase the number of points N for FFT computation as follows fft(u,N).
7. With low number of time samples, you can not detect both frequencies f1 and f2, you need larger number of samples, this is related to Heisenberg uncertainty principle in signal processing, as example for n=310, you can detect the frequencies.
clear all;close all;
f1=0.115; f2=0.135; n=0:310-1;
%PHASE=2*pi*rand();
w=0.1*randn(size(n));
u=real(exp(j*2*pi*f1*n))+real(exp(j*2*pi*f2*n))+w;
%Using FFT
N=600;
U=2*fft(u,N)/N;
fvals=(0:N-1)/N;
mag=abs(U);
subplot(1,2,1);plot(fvals,mag);
xlabel('Frenquency'); ylabel('P(f)'); title('Power spectrum by FFT method');
subplot(1,2,2);plot(fvals,10*log10(mag.^2));
xlabel('Frequency'); ylabel('P(f)'); title('Power spectrum by FFT method in dB');
##### 2 CommentsShowHide 1 older comment
Youssef Khmou on 1 Feb 2015
You can scale by any chosen factor in order to get amplitude or power (std(u)), the fft result contain negative and positive frequencies if you use the frequency axes as :
fvals=(-N/2:N/2-1)*Fs/N; % Fs is the sampling rate which is 1 in this case.
mag=fftshift(abs(U));
In the example above you obtain the correct frequencies f1 and f2, plus a replicate at around f'1=0.86 and 0.88.