MATLAB: Power cross-spectral density

cross-spectral densityfftfrequency domainpower spectral density

Hi all,

I am curently working with frequency response functions on the basis of 1 input data sample and 1 output data sample. As I would like to plot the coherence function, I need to compute the cross-spectral densities. Since I am new to this field, I would like to ask some of you experts whether it is correct to compute a power cross-spectral density (pcsd) as

pcsd_xy = psd_xx.*conj(psd_yy)

where psd = power spectral density. Is this correct? I know that MATLAB has commands which can solve this for me, but I would really like to programme it myself.

Thanks in advance.

Kind regards, Martin Ulriksen.

Best Answer

That is not correct. The cross spectrum (the Fourier transform of the cross correlation) is not real-valued. The cross spectrum is (ignoring the scaling factors) the DFT of X times the complex conjugate of the DFT of Y.

The way you have it written there is no need to do conj(psd_yy) because the PSD of y is real-valued.

The cross-correlation of X and Y is not an even function of the lag, k, unlike the autocorrelation of X OR the autocorrelation of Y. That is why the Fourier transform of the cross-correlation (the cross spectrum) is complex-valued.

Related Solutions

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

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.

MATLAB: Normxcorr2 function – spatial / frequency domain selection

Hi Sam, normxcorr2 computes an estimate of the computation time for spatial domain correlation vs. frequency domain. Whichever is the smallest wins.

You could save normxcorr2.m as a different name and then modify the code so that it always did one or the other. For example, to always use spatial domain, you can modify the xcorr2_fast function inside of normxcorr2 as follows

function cross_corr = xcorr2_fast(T,A)
T_size = size(T);
A_size = size(A);
outsize = A_size + T_size - 1;
% figure out when to use spatial domain vs. freq domain

%conv_time = time_conv2(T_size,A_size); % 1 conv2

%fft_time = 3*time_fft2(outsize); % 2 fft2 + 1 ifft2

%if (conv_time < fft_time)

    cross_corr = conv2(rot90(T,2),A);
%else

%    cross_corr = freqxcorr(T,A,outsize);
%end

To only use the Fourier transform:

function cross_corr = xcorr2_fast(T,A)
T_size = size(T);
A_size = size(A);
outsize = A_size + T_size - 1;
% figure out when to use spatial domain vs. freq domain
%conv_time = time_conv2(T_size,A_size); % 1 conv2
%fft_time = 3*time_fft2(outsize); % 2 fft2 + 1 ifft2
%if (conv_time < fft_time)
%    cross_corr = conv2(rot90(T,2),A);
%else
    cross_corr = freqxcorr(T,A,outsize);
%end

Best Answer

Related Solutions

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

MATLAB: Normxcorr2 function – spatial / frequency domain selection

Related Question