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.

http://www.mathworks.cn/support/solutions/en/data/1-AO3A1Q/index.html?product=SG&solution=1-AO3A1Q