power spectral density PSD?

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?

Answer 1

[Youssef Khmou]

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)

[Mary Jon]

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]

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]

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]

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

Answer 2

[Wayne King]

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)))

[Mary Jon]

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]

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]

which are values =zero ? values of time you are meaning?

[Wayne King]

your data values are almost zero They are all on the order of 10^(-16)

[Youssef Khmou]

you can also try :

 T=abs(fft(c13,512));
T=T.*conj(T);
figure, plot(T)

[Mary Jon]

what is 512 in your code?

[Youssef Khmou]

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