What is the difference between CSD and CPSD in Signal Processing Toolbox 6.11 (R2009a)?

I was using CSD from the Signal Processing Toolbox to compute power spectral densities but since Release R14 I am receiving a warning that CSD has been replaced by CPSD and CSD still works but may be removed in future.

I was trying to adapt my code to use CPSD however I was not able to produce the same results.

What is the difference between CSD and CPSD and how can produce the same results using CPSD?

Answer 1

[MathWorks Support Team]

CPSD returns the power spectral density, while CSD returns neither power nor power spectral density. Depending on the application, this difference may be critical.

In addition, when calling CSD(x,y), it returns the FFT of autocorrelation E(yx*) while CPSD(x,y) returns the FFT of autocorrelation E(xy*). This explains the phase inverse observed in the result. CPSD's approach is more natural if one wants to find out the cross power spectral density between x and y.

In summary, to get CSD result with CPSD, please do the following:

1. Use Hanning window instead of Hamming (default).

2. Set 'noverlap' to zero.

3. Scale the result with Fs.

4. Pass the input arguments in opposite order.

Here is some additional information regarding step 3 - how to scale the result. In CSD, if the signal is real, the function returns the positive side of the

spectrum, period. However, in CPSD, it does a little more to come up with

this one-sided power spectrum. Specifically, the power from the negative

frequencies are added to the positive frequencies. However, if the number of

FFT is even, there are two unique frequency points, i.e., the DC and the

Nyquist point (Fs/2). Therefore, magnitude in all points are doubled except

these two points.

Here is an example on how to produce the same result as CSD using CPSD:

t=0:0.00025:3; %3s @ 4kHz sample rate

y=chirp(t,100,1,1000);

semilogx(t,y)

%CSD

[Pcsd F1]=csd(t,y, 2^12,4000);

%CPSD

window=hanning(4096);

noverlap=0;

[Pcpsd,F2]=cpsd(y,t,window,noverlap,2^12, 4000);

Pcpsd = Pcpsd*4000;  % scale with Fs

Pcpsd(2:end-1) = Pcpsd(2:end-1)/2; 

% except the two end points, everything is 

% already doubled, which is not present in CSD

sum(Pcpsd-Pcsd)

Note the difference is almost zero - the results can not be exactly equal due

to different mathematical code paths and consequent floating-point round

off differences.