22

I asked a question earlier but I didn't get any answer for it. So now I am simplifying it: what are Cross-Spectral Density (CSD) and Power-Spectral Sensity (PSD)? What is their application? How can I get them in MATLAB?

$$S_{kl}(\omega)=\lim_{T\to\infty}\frac{1}{T}E\{Y_k^*(\omega)Y_l(\omega)\} $$ $$S_{kk}(\omega)=\lim_{T\to\infty}\frac{1}{T}E\{Y_k^*(\omega)Y_k(\omega)\} $$

$S_{kl}(\omega)$ is the cross-spectral density (CSD) function between general signals $y_k(t)$ and $y_l(t)$, $S_{kk}(\omega)$ is the power-spectral density (PSD) of signal $y_k(t)$, $Y_k(\omega)$ is the finite Fourier transform of signal $y_k(t)$ at frequency $\omega$, $Y_k^*(\omega)$ is the complex conjugate of $Y_k(\omega)$, and $E\{\cdot\}$ is the expectation operator.

My earlier question was: What does 'wavelet power spectrum', 'Auto-power spectrum','cross-power spectrum' means in wavelet application? I was studying about mode shape identification with wavelet method and these terms confused me.

jojeck
  • 11,107
  • 6
  • 38
  • 74
SAH
  • 1,317
  • 4
  • 17
  • 39
  • Can you post a reference to some of the material that you've been studying? It's easier to help you if you do. – Phonon Aug 19 '13 at 18:35
  • @Phonon Hi phanon. I edit my question and post the link. Can you guys access to the paper or you want me to upload it somewhere? tnx – SAH Aug 19 '13 at 19:00

2 Answers2

16

Power-Spectral Density (PSD) is the distribution of power along the frequency axis. It is generally used for non-finite energy signals (mostly not limited in time signals), who aren't square-summable. The signal's PSD is the Fourier Transform of the signal's autocorrelation, as stated by the Wiener–Khinchin theorem. In Matlab:

N = length(S);
F = fft(S);
F = F(1:N/2+1);
PSD = (1/(2*pi*N)) * abs(F).^2;
PSD(2:end-1) = 2*PSD(2:end-1);
freq = 0:(2*pi)/N:pi;

See: https://de.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

Cross-Spectral Density is the same, but using cross-correlation, so you can find the power shared by a given frequency for the two signals using its squared module, and the phase shift between the two signals at that frequency using its argument.

Cross-Spectral Density can be used to identify the frequency response of a noisy LTI system : if the noise is not correlated to the input or output of the system, its frequency response can be found from the CSD of the input and output.

lennon310
  • 3,590
  • 19
  • 24
  • 27
Florian Castellane
  • 539
  • 1
  • 5
  • 14
  • ,Thanks for your answer, would you write matlab code for CSD too please? And would you write a example of CSD to identify frequency response of a noisy LTI system? – SAH Oct 29 '13 at 13:02
  • @Electricman The MATLAB Signal Processing Toolbox already has functions to do that. In particular, cpsd() does what you need. – Phonon Oct 29 '13 at 18:19
  • @Phonon, I think that uses FFT. how can I run a CSD by wavelet transform? Thanks Phonon – SAH Oct 29 '13 at 18:58
  • @Electricman You should ask that as a separate question. – Phonon Oct 29 '13 at 19:01
  • @Phonon, If someone write the FFT based CSD code in matlab. I can do wavelet based CSD myself.cpsd() function doesn't help me. Thanks loads – SAH Oct 29 '13 at 19:14
  • What I meant wa that you should create a separate question on the site, instead of asking it in comments to this question. – Phonon Oct 29 '13 at 19:24
  • @Phonon,The answer is fine and I accepted the question, But I would like to know detail matlab code of CSD and have a example of identifying frequency response of a noisy LTI system by CSD. so I may ask this question separately soon, But still I say a complete answer on this post is better. thanks loads. – SAH Oct 29 '13 at 22:03
  • I don't have an example I'm sure of (and correct me if I'm wrong) but if you take the FFT of the cross-correlation between the input and the output, you get the CSD of the input and output. Note that this is equal to the product of frequency response of the system and the PSD of the input. So you can divide the CSD by the PSD to get an estimate of the frequency response, the less the noise is coherent with the input the better the estimate. – Florian Castellane Oct 30 '13 at 07:37
  • @FlorianCastellane, Thanks for ur comment, I'm working on sth else now, so will work on CSD later and ask my questions here. best regard. – SAH Oct 30 '13 at 16:08
  • 3
    "The signal's PSD is the autocorrelation of the signal's Fourier Transform": the other way around? $PSD(x,\omega) = FT(acorr(x,t), \omega)$. – scrx2 Dec 06 '18 at 09:23
  • in the above code, the FFT calculation should be corrected as F = fft(S)/N; https://electronics.stackexchange.com/questions/25900/scaling-fft-output-by-number-of-points-in-fft – Kasun Jun 16 '19 at 10:46
7

To add to the above well stated explanation, in the case of wavelets, which are finite in time, it is more correct not to use the term 'power' but 'energy'. For Fourier who has as basis functions the sinusoid that extends infinite in time, power spectral density is the correct term. For wavelets, who has basis functions as finite in time deflections, we should use 'energy'.

forsker_for_dsp
  • 171
  • 2
  • 4