How to get the phase shift, the amplitude ratio and cross power spectral density Pxy between two signals?
amanda
Hi, I am new to matlab and signal processing.
I have a problem:
I have 2 signals, one is regarded to be the input (x) and the other is the output (y). They have the same frequency. The relation between input x and output y is given by the frequency response function H, which is the Fourier transform of the unit impulse response function of the system. And the frequency response function can be calculated optimally from the input power spectral density Pxx and the cross power spectral density Pxy as
H(f)=Pxy(f)/Pxx(f）
Generally speaking, H(f) is a complex function and it can be given in polar notation as
H(f)=|H（f)|e^(-i*a(f))
where a(f) is the phase shift and |H（f)| is the amplitude ratio.
The modulus of H at frequency f is the amplitude ratio, could be estimated from the auto power spectral densities by
|H（f)|=(Pyy(f)/Pxx(f))^0.5
And the phase a(f) is the phase shift of the output signal compared to the input signal...
So how could I get the phase shift a(f), amplitude ratio |H(f)| and H?
I tried this codes, I used two different methods to make H and phase shift, but I dont know if it is available..
% PSD using Welch's method
clear;
clc;
Fs=1000;
s=xlsread('123.xlsx');
x1=s(:,2);
x2=x1-mean(x1);
nfft=1024;
window=hanning(1024);
noverlap=512;
[Pxx,f]=pwelch(x2,window,noverlap,nfft,Fs);
plot(f,Pxx)
x3=s(:,4);
x4=x3-mean(x3);
nfft=1024;
window=hanning(1024);
noverlap=512;
[Pyy,f]=pwelch(x4,window,noverlap,nfft,Fs);
plot(f,Pyy)
[Pxy,f] = cpsd(x2,x4,window,noverlap,nfft,Fs);
plot(f,Pxy);
[Txy,f] = tfestimate(x2,x4,window,noverlap,nfft,Fs);
H=norm(Txy);
H2=norm(Pxy);
Phaseshift=angle(log(Txy/H)/(-i))
Phaseshift2=angle(log(Pxy/H2)/(-i))
thanks for help!
Regards!
amanda
In the codes, the Txy is equal to the frequency response H.
and H and H2 are both equal to its modulus |H|.
Phaseshift and phaseshift2 are both equal to the phase shift a(f).
These three value are needed.
Thanks again.
Alexander Sotnikov
to apply 'tfestimate'. It gives you an estimate of the transfer
function. The amplitude response of the system is the absolute value of
Txy in your code, and the phase response is the angle between real and
imaginary components of Txy. So
Amplitude = 20*log10(abs(fftshift(Txy)))
Phase = angle(fftshift(Txy))
And if you need to estimate phase shift and amplitude ratio at a
particular frequency, then simply pick the corresponding elements from
Amplitude and Phase arrays.
HTH
--
Alexander
amanda
Hi, Alexander, thanks a lot.
I am not sure 'tfestimate' means the frequency response H (corresponding to the equation 1, H(f)=Pxy(f)/Pxx(f）). and with another two equations (H(f)=|H（f)|e^(-i*a(f)) , |H（f)|=(Pyy(f)/Pxx(f))^0.5), I wanna get the a(f).
1.By using the 'tfestimate', I got the H(f)(which is equal to Txy).
2.And is the amplitude ratio's code not like this ? |H(f)|=abs(H(f)). What you wrote looks like signal-to-noise ratio (HNR).
3.To get the phase shift (a(f)), I solved the equation (H(f)=|H（f)|e^(-i*a(f)), and then got a(f). is 'fftshift ' equal to this equation?
I am really sorry I have so many questions... Thanks for your help
Alexander Sotnikov
> I am not sure 'tfestimate' means the frequency response H (corresponding
> to the equation 1, H(f)=Pxy(f)/Pxx(f）).
quotient of Cross Power Spectral Density (CPSD) of X and Y, Pxy and the
Power Spectral Density (PSD) of X, Pxx. It can be represented in polar
coordinates as H(f) = abs(H(f))*exp(j*angle(H(f))).
>
> 2.And is the amplitude ratio's code not like this ? |H(f)|=abs(H(f)).
Yes, it is
> What you wrote looks like signal-to-noise ratio (HNR).
I have just converted it to decibel, which is more convenient, but not
absolutely necessary
3.To get the
> phase shift (a(f)), I solved the equation
> (H(f)=|H（f)|e^(-i*a(f)), and then got a(f). is 'fftshift ' equal
> to this equation?
No. I assumed complex-valued x and y, in which case 'tfestimate' returns
transfer function estimate for the frequencies between 0 and Fs. The
'fftshift' is used to swap the two halves of the spectrum when working
with functions like 'fft' and so on.
However, for real-valued data 'tfestimate' returns one-sided frequency
response between 0 and Fs/2, so if your data sequences are real-valued
you shouldn't use 'fftshift', because it will give you erroneous
results. Simply apply 'angle' function to the TF array. Be careful with
interpreting the results because the output of 'angle' function is
defined between -pi and pi, so you may need to unwrap it.
--
Alexander
amanda
Thanks a lot. I got it.