How Do You Compute Power Cross-Spectral Density for Coherence Functions?

[MartinDU]

Homework Statement

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

Homework Equations

pcsd_xy = psd_xx.*conj(psd_yy)

where psd = power spectral density.

The Attempt at a Solution

Is this a correct approach to obtain the cross-spectral density?

Kind regards,
Martin.

 

[marcusl]

No, that's not going to do it. Look at the definition of the cross-power PSD from the Wiener-Khinchin theorem
S_{xy}(\omega)=\int_{-\infty}^{\infty}{R_{xy}(\tau)e^{-i\omega \tau} d\tau},where R is the cross-correlation function
R_{xy}(\tau)=\int_{-\infty}^{\infty}{x(t) y^*(t+\tau)dt}.You can see that you can't get the above from what you wrote (note that your PSD S_xx is found by setting x=y).

 

[MartinDU]

Thanks for the answer. Then I have probably made an error on an earlier state since I obtain compliance between the frequency response function when I use

S_xx(w)=|H(w)|^2 S_ff(w)

and

S_xx(w)=H(jw)S_xf(w).

I have posted a small MATLAB script in which I compute these quantities. I know it is a lot to ask for but can you spot something wrong with my computation of the PSDs?

% Time increment %
dt=t(2)-t(1);

% Number of samples %
N=2^nextpow2(length(t));

% Frequency band %
f=1/dt/N*(0:(N-2)/2);

% --- Fast Fourier transformation --- %

% Input data %
Ft_in=fft(wdata_in,N);

% Output data %
xt_out=fft(wdata_out,N);

% --- Power spectral densities --- %

% Power spectral density of input data %
PSD_in=Ft_in.*conj(Ft_in)/N;

% Power spectral density of output data %
PSD_out=xt_out.*conj(xt_out)/N;

% Power cross-spectral densities of input and output data %
PCSD_in_out=Ft_in.*conj(xt_out)/N;
PCSD_out_in=xt_out.*conj(Ft_in)/N;

% Frequency response function %
H=PSD_out(1:N/2)./PCSD_out_in(1:N/2);

Your answer would be highly appreciated.

Thanks in advance and kind regards,
//Martin.

 

[marcusl]

Your input and output PSD's are correct, your cross power calculation is not. You are getting a pleasing result because cross-power would be the wrong thing to use, were you to calculate it correctly. The transfer function is H(\omega)=\frac{Y(\omega)}{X(\omega)}, where Y and X are the output and input spectra of your system. Look at your code--you in fact calculate H*, albeit through a very convoluted process.

 

[MartinDU]

Again, thank you very much for the answer.

The reason why I want to use the cross-power is that I want informations about the phase angle as well as the magnitude.
So what you are saying is that I through

H=PSD_out(1:N/2)./PCSD_out_in(1:N/2)

actually get the complex conjugate of H?

 

[marcusl]

Again, thank you very much for the answer.

The reason why I want to use the cross-power is that I want informations about the phase angle as well as the magnitude.
So what you are saying is that I through

H=PSD_out(1:N/2)./PCSD_out_in(1:N/2)

actually get the complex conjugate of H?

Take a look at your own code:

PSD_out=xt_out.*conj(xt_out)/N;

PCSD_out_in=xt_out.*conj(Ft_in)/N;

% Frequency response function %
H=PSD_out(1:N/2)./PCSD_out_in(1:N/2);

What does it look like to you?