Compute the power spectrum using FFT method

48 visualizaciones (últimos 30 días)
An tran
An tran el 31 de En. de 2015
Editada: Walter Roberson el 9 de Abr. de 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');

Respuesta aceptada

Youssef  Khmou
Youssef Khmou el 1 de Feb. de 2015
Editada: Youssef Khmou el 1 de Feb. de 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 comentarios
An tran
An tran el 1 de Feb. de 2015
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 el 1 de Feb. de 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.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Fourier Analysis and Filtering en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by