Compute the power spectrum using FFT method

51 次查看(过去 30 天)
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
Youssef Khmou 2015-2-1
编辑:Youssef Khmou 2015-2-1
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 个评论
An tran
An tran 2015-2-1
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
Youssef Khmou 2015-2-1
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.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by