|
Hi there, I need to find the power spectral densities of the following signals: t = [0:1/1800:1]; x=sin(2*pi*300*t)+sin(2*pi*330*t); y=sin(2*pi*700*t)+sin(2*pi*900*t); X = fft(x,8192); Y = fft(y,8192); f=[0:8191]/8192; Pyy=Y.*conj(Y)/abs(Y); Pxx=X.*conj(X)/abs(X); when i plot (f,Pyy ) , the resulting graphich seems to be wrong.. Could you please tell me what am i doing wrong?.. Thank you very much Berra |
|
|
Power Spectral Density
Started by ●October 16, 2002
Reply by ●October 16, 20022002-10-16
|
Hi, You probably wanted to do Pyy=Y.*conj(Y)./abs(Y); Pxx=X.*conj(X)./abs(X); Note the 'dot' before / It is element by element division. But Pxx = abs(X); Pyy = abs(Y); will do the same faster. Navan --- berra tosun <> wrote: > Hi there, > > I need to find the power spectral densities of the > following signals: > > t = [0:1/1800:1]; > x=sin(2*pi*300*t)+sin(2*pi*330*t); > y=sin(2*pi*700*t)+sin(2*pi*900*t); > X = fft(x,8192); > Y = fft(y,8192); > > f=[0:8191]/8192; > > Pyy=Y.*conj(Y)/abs(Y); > Pxx=X.*conj(X)/abs(X); > > when i plot (f,Pyy ) , the resulting graphich seems > to be wrong.. Could you please tell me what am i > doing wrong?.. > > Thank you very much > > Berra > ------------------------ Yahoo! Groups Sponsor > > _____________________________________ > Note: If you do a simple "reply" with your email > client, only the author of this message will receive > your answer. You need to do a "reply all" if you > want your answer to be distributed to the entire > group. > > _____________________________________ > About this discussion group: > > To Join: > > To Post: > > To Leave: > > Archives: http://www.yahoogroups.com/group/matlab > > More DSP-Related Groups: > http://www.dsprelated.com/groups.php3 > > ">http://docs.yahoo.com/info/terms/ __________________________________________________ |
Reply by ●October 17, 20022002-10-17
|
Hi Berra! First of all, 1/1800 is too low sampling frequency for the second y component sin(2*pi*900*t), so change > t = [0:1/1800:1]; to > t = [0:1/8191:1]; Then, try with abs instead of > Pyy=Y.*conj(Y)/abs(Y); > Xa = abs(X); > plot(f,Xa); > Ya = abs(Y); > plot(f,Ya); Pyy as previously calculated was a constant, not an array to be plotted. Regards Predrag ----- Original Message ----- From: berra tosun <> To: <> Sent: 2002. listopad 16 20:33 Subject: [matlab] Power Spectral Density > Hi there, > > I need to find the power spectral densities of the following signals: > > t = [0:1/1800:1]; > x=sin(2*pi*300*t)+sin(2*pi*330*t); > y=sin(2*pi*700*t)+sin(2*pi*900*t); > X = fft(x,8192); > Y = fft(y,8192); > > f=[0:8191]/8192; > > Pyy=Y.*conj(Y)/abs(Y); > Pxx=X.*conj(X)/abs(X); > > when i plot (f,Pyy ) , the resulting graphich seems to be wrong.. Could you please tell me what am i doing wrong?.. > > Thank you very much > > Berra > > _____________________________________ > Note: If you do a simple "reply" with your email client, only the author of this message will receive your answer. You need to do a "reply all" if you want your answer to be distributed to the entire group. > > _____________________________________ > About this discussion group: > > To Join: > > To Post: > > To Leave: > > Archives: http://www.yahoogroups.com/group/matlab > > More DSP-Related Groups: http://www.dsprelated.com/groups.php3 > > ">http://docs.yahoo.com/info/terms/ > |
Reply by ●October 17, 20022002-10-17
|
Hi Berra, There are numerous ways to estimate power spectral density - and I recommend you look into which suits your needs best. But for now to answer your question. There are a couple of errors in your code. Most importantly, your sampling frequency is a little too small. As a rule of thumb allow a bit more than the Nyquist rate (2.5 times the maximum frequency should be enough). The PSD is a measurement at various frequencies, so: Pyy = Y.*conj(Y)/8192; % estimated over the FFT size Now to graph, you only need to plot the first 4096 points, as the other 4096 points are symmetric. Use: f=Fs*(0:4095)/8192; plot(f,Pyy(1:4096)) Put this all together gives: Fs'00; t = [0:1/Fs:1]; x=sin(2*pi*300*t)+sin(2*pi*330*t); y=sin(2*pi*700*t)+sin(2*pi*900*t); X = fft(x,8192); Y = fft(y,8192); f=Fs*(0:4095)/8192; Pyy=Y.*conj(Y)/8192; Pxx=X.*conj(X)/8192; plot(f,Pyy(1:4096)) hold plot(f,Pxx(1:4096),'g') You should see your spectral peaks at 300Hz, 330Hz (in green) and 700Hz, 900Hz in blue. I would also suggest you use a window to reduce spectral leakage, something like: W=blakman(length(y)); Yt((y.*W'),8192); You can then compare the windowed and non-windowed responses. Matlab has several PSD estimation algorithms, type > help psd at the matlab prompt for details. Hope this helps, Jeff -----Original Message----- From: berra tosun [mailto:] Sent: Wednesday, October 16, 2002 7:33 PM To: Subject: [matlab] Power Spectral Density Hi there, I need to find the power spectral densities of the following signals: t = [0:1/1800:1]; x=sin(2*pi*300*t)+sin(2*pi*330*t); y=sin(2*pi*700*t)+sin(2*pi*900*t); X = fft(x,8192); Y = fft(y,8192); f=[0:8191]/8192; Pyy=Y.*conj(Y)/abs(Y); Pxx=X.*conj(X)/abs(X); when i plot (f,Pyy ) , the resulting graphich seems to be wrong.. Could you please tell me what am i doing wrong?.. Thank you very much Berra _____________________________________ Note: If you do a simple "reply" with your email client, only the author of this message will receive your answer. You need to do a "reply all" if you want your answer to be distributed to the entire group. _____________________________________ About this discussion group: To Join: To Post: To Leave: Archives: http://www.yahoogroups.com/group/matlab More DSP-Related Groups: http://www.dsprelated.com/groups.php3 ">http://docs.yahoo.com/info/terms/ |
