# Compute the power spectrum using FFT method

37 次查看（过去 30 天）
An tran2015-1-31
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 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 个评论显示 1更早的评论隐藏 1更早的评论
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.

### 类别

Find more on Digital Filtering in Help Center and File Exchange

### Community Treasure Hunt

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

Start Hunting!