-1
$\begingroup$

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); %******
Improve this question
$\endgroup$
9
  • $\begingroup$ Hint: because the cross power spectral density is NOT equal to $\sqrt{S_{xx} \cdot S_{yy}}$ ... $\endgroup$
    – user883521
    Dec 4, 2017 at 17:39
  • $\begingroup$ 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. $\endgroup$
    – Matt Vowels
    Dec 4, 2017 at 17:41
  • $\begingroup$ 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. $\endgroup$
    – user883521
    Dec 4, 2017 at 17:47
  • $\begingroup$ 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. $\endgroup$
    – Matt Vowels
    Dec 4, 2017 at 17:48
  • 1
    $\begingroup$ 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. $\endgroup$
    – user883521
    Dec 6, 2017 at 19:37

2 Answers 2

Reset to default
0
$\begingroup$

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);
Improve this answer
$\endgroup$
0
$\begingroup$

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;
samplesread = nfft-1;
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 
endpoint=startingpoint+samplesread; %read 255 data points every time
x1=resultx(startingpoint:endpoint)'.*window;
y1=resulty(startingpoint:endpoint)'.*window;
xdft1 = fft(x1,nfft);  
ydft1 = fft(y1,nfft);
startingpoint=startingpoint+(samplesread+1)/2;
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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.