power spectral density PSD?

2013 11 28
2 个回答
回答已采纳
更新时间：2013 12 3
26 次查看（30 天）

0 个投票

If I am have signal with length(33),or 13 signals each with length(33), How finding PSD to each signal individually?and plot its individually?

请先登录，再回答此问题。

 采纳的回答

Youssef Khmou
Youssef Khmou 2013-11-29

2 个投票

You can try this way :
t=linspace(0,1,33);% 1 seconde
Fs=inv(t(3)-t(2));
f=Fs/10;
P=13; % number of signals
X=zeros(33,P);
for n=1:P
X(:,n)=sin(2*pi*t*(f+n));
end
N=512; % number of points for computing DFT
frequency=(0:N-1)*Fs/N; % frequency axis
frequency=frequency(1:floor(end/2)); % One sided spectrum
PSD=zeros(N,P);
for n=1:P
PSD(:,n)=fft(X(:,n),N);
PSD(:,n)=PSD(:,n).*conj(PSD(:,n));
end
PSD(floor(end/2)+1:N,:)=[]; % one sided spectrum
% Plotting them all in one figure
figure,plot(frequency,PSD)

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Mary Jon
Mary Jon 2013-11-30
Thank you Youssef for answer,
execuse me ,I will only change X=zeros(33,p) by my 13 signals ,
and ,what is floor in your code?
frequency=frequency(1:floor(end/2));
Youssef Khmou
Youssef Khmou 2013-12-1
suppose that the length of frequency vector is 201, end/2=100.5 which is false the length should be an integer, use floor or ceil or fix, clear?
Mary Jon
Mary Jon 2013-12-3
Hello Youssef I note in your code above ,there is sin(),my data or signals is not pure sin,why you are used sin()?
X(:,n)=sin(2*pi*t*(f+n));
If this is just example to show me how find PSD,how can I modified this code ,to finding PSD to my signals of length (33)?
when I am used this code
[Pxx(:,nn),Fxx] = periodogram(X(:,nn),[],64,1);
the final result not same result of your code?when plot two result of codes I get different plots?
Youssef Khmou
Youssef Khmou 2013-12-3
mary, they are the same, it is just a problem of scale, Fxx is the frequency axis and it is 33x1, then we should get the same number of points using fft:
plot(Fxx,Pxx);
hold on;
F1=fft(X(:,1),33*2); % 33*2 points
F1=abs(F1(1:33));
plot(Fxx,F1,'r');
I used sinus just as example, Good luck mary

请先登录，再进行评论。

更多回答（1 个）

Wayne King
Wayne King 2013-11-28

1 个投票

Hi Mary, If you have the Signal Processing Toolbox, the easiest thing is to use periodogram()
I'll create some simulated signals. I'll assume your sampling frequency is 1.
X = randn(33,13);
for nn = 1:13
[Pxx(:,nn),Fxx] = periodogram(X(:,nn),[],64,1);
end
You can plot each one individually by selecting the column.
plot(Fxx,10*log10(Pxx(:,1)))

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Mary Jon
Mary Jon 2013-11-28
c13=[4.0342e-016 4.0337e-016 4.0337e-016 4.0338e-016 4.034e-016 4.0342e-016 4.0344e-016 4.0347e-016 4.035e-016 4.0353e-016 4.0357e-016 4.036e-016 4.0363e-016 4.0365e-016 4.0367e-016 4.0368e-016 4.0369e-016 4.0368e-016 4.0367e-016 4.0365e-016 4.0363e-016 4.036e-016 4.0357e-016 4.0353e-016 4.035e-016 4.0347e-016 4.0344e-016 4.0342e-016 4.034e-016 4.0338e-016 4.0337e-016 4.0337e-016 4.0342e-016];
This is one of 13 signal ,how can specified sampling frequency? what mean 64,1 in your code
[Pxx(:,nn),Fxx] = periodogram(X(:,nn),[],64,1);
Wayne King
Wayne King 2013-11-28
64 is zero-padding because your signals are pretty short in length. To know the sampling frequency you have to know how the data were acquired. What is the time spacing between each value in c13 for example.
By the way, these values are almost zero 10^(-16)
Mary Jon
Mary Jon 2013-11-28
which are values =zero ? values of time you are meaning?
Wayne King
Wayne King 2013-11-28
your data values are almost zero They are all on the order of 10^(-16)
Youssef Khmou
Youssef Khmou 2013-11-29
you can also try :
T=abs(fft(c13,512));
T=T.*conj(T);
figure, plot(T)
Mary Jon
Mary Jon 2013-11-30
what is 512 in your code?
Youssef Khmou
Youssef Khmou 2013-12-1
512 is the number of points for calculating DFT, you can put any number, higher number gives good resolution, any number that is multiple of 2 (128,512,1024...) is faster DFT becomes FFT

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by