Compute the power spectrum using FFT method

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');

Answer 1

[Youssef Khmou]

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');

[An tran]

1.Why do you multiply by 2 and divide by N in " U=2*fft(u,N)/N;"?

2. In the plot, why are there 2 peaks at the frequency 0.8<f<0.9. Isn't it supposed to be only two peaks at the f=0.115 and f=0.135?

Thanks

[Youssef Khmou]

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.