# MATLAB's $\tt cpsd$ and $\tt pwelch$ - different results for cross spectral power density

I am wondering why the following code does not yield the same results for Sxy1 and Sxy2, where Sxy1 is calculated using the MATLAB's pwelch function, and Sxy2 is calculated using cpsd function.

fs = 1000;
t = 0:1/fs:5-1/fs;
win=hann(500);
overlap=300;
nfft=500;
fs=1000;

x = cos(2*pi*100*t)+randn(size(t));
y=sin(2*pi*160*t)+cos(2*pi*150*t)+randn(size(t));

[Sxx,f] = pwelch(x,win,overlap,nfft,fs);
[Syy,f] = pwelch(y,win,overlap,nfft,fs);

Sxy1=sqrt(Sxx.*Syy); %******

[Sxy2,f] = cpsd(x,y,win,overlap,nfft,fs);

Sxy2=abs(Sxy2); %******

• Hint: because the cross power spectral density is NOT equal to $\sqrt{S_{xx} \cdot S_{yy}}$ ... – user883521 Dec 4 '17 at 17:39
• Thanks! This is what I suspected. Do you know the equation for cross power spectral density in terms of power spectra? I believe, for amplitude spectra: Sxy = (1/L) .\times (X(ejw) .\times conj(Y(ejw)) is correct? Where X(ejw) and Y(ejw) are the complex spectra of the x and y signals. – Matt Vowels Dec 4 '17 at 17:41
• Did you compare the Sxy above to the output of cpsd? The Sxy you specify in the above comment is an estimate of the cross power spectral density (depending on how the DFT is scaled an additional scaling factor might be required). To get the same result as cpsd make sure you use the same window, overlap, number of averages etc. – user883521 Dec 4 '17 at 17:47
• Unfortunately, there is no function that I am aware of that outputs complex amplitude spectra using Welch's method (only power spectra). I am currently in the process of coding one that does, but it would be easier if I had an equation for CPSD in terms of (real valued) power spectra. – Matt Vowels Dec 4 '17 at 17:48
• Welch’s method uses overlapped, averaged FFT’s and their complex conjugates to estimate the PSD. You cannot use this method to directly calculate an amplitude spectrum as, in general, you cannot average the complex FFT (converges to zero if signals are noisy). You can calculate an amplitude spectrum from a PSD estimated using Welch’s method (AS = sqrt(2*PSD*df)) but, as the PSD, this will not contain any phase information and hence cannot be used to calculate the cpsd. – user883521 Dec 6 '17 at 19:37

Whilst I have not worked out how to calculate CPSD from power spectra, I have confirmed the following two methods are equivalent (using amplitude spectra, rather than power spectra or cross-correlation):

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));
y=sin(2*pi*160*t)+cos(2*pi*150*t)+randn(size(t));
N = length(x);
nfft=500;
lfft=2*nfft;
overlap=0;
win=rectwin(N);

xdft1 = fft(x,2*nfft);  <outputs two-sided complex amplitude spectra
ydft1 = fft(y,2*nfft);
Sxy1a=xdft1.*conj(ydft1);
Sxy1b=Sxy1a./N^2;
Sxy1c=Sxy1b(1:lfft/2+1);
Sxy1c(2:end-1)=2*Sxy1c(2:end-1);
Sxy1=abs(Sxy1c)';

%CPSD Method 2
Sxy2a=cpsd(x,y,win,overlap,2*nfft,Fs);   < uses the cpsd function with Welch’s overlapped method
Sxy2=abs(Sxy2a);


I reproduced the last piece of code and it works. I'll just add my piece of code to show how it works with any Fs and fft points but only 50% overlap.

resultx = EEG.data(5,1:16384);
resulty = EEG.data(3,1:16384);
Fs = 250;
Navg = 127;
nfft = 256;
window = hamming(nfft);
xmulconjyfft = zeros(nfft/2+1,Navg);
for i=1:Navg % number of iterations, that depends on fft points and the overlap, in my case with fftpoints=256, 50% overlap and inputlength=2^14, the number of iterations is 127
x1=resultx(startingpoint:endpoint)'.*window;
y1=resulty(startingpoint:endpoint)'.*window;
xdft1 = fft(x1,nfft);
ydft1 = fft(y1,nfft);
Sxy1a=ydft1.*conj(xdft1);
Sxy1b=Sxy1a./N^2;
Sxy1c=Sxy1b(1:nfft/2+1);
Sxy1c(2:end-1)=2*Sxy1c(2:end-1);
Sxy1c = Sxy1c';
xmulconjyfft(:,i)=Sxy1c;
%Sxy1=abs(Sxy1c);
end

avgxyfft = mean(xmulconjyfft,2); %average the spectrums
avgxyfft= avgxyfft*((length(resultx).^2)/(nfft*Fs)); % normalization factor

Sxy2a=cpsd(resultx,resulty,window,[],nfft,Fs);


From code above Sxy2a and avgxyfft result in two identical vectors with complex values.