PSDs by code compared to cpsd command

Question

Martin 2012-11-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/52543-psds-by-code-compared-to-cpsd-command

Hi experts,

I have tried to make a code to generate, e.g., cross-spectral density of an input and output signal with the use of Welch's method and different windows (e.g. Hann and Hamming). However, when I compare this to the command cpsd my results differ in the first components (not just the first component) and the last components (not just the last component). The code I've made looks like this:

function [frf_mag,frf_phase,f] = spectral_analysis(t,x,y,window,seg_num,seg_ov);

if true
dt=t(2)-t(1);
Fs=1/dt;
fn=Fs/2; 
real_v=isfinite(x);
x=x(real_v); 
x=x(:);
y=y(real_v);
y=y(:);
t=t(real_v);
N=sum(real_v);
L=ceil(N/(seg_num-seg_num*seg_ov+seg_ov));
ov_num=floor(seg_ov*L);
w=eval([window,'(',num2str(L),')']);
R=w'*w;  
nfft=2^nextpow2(L);
f=(1:nfft/2)/(nfft/2)*fn; 
f=f(:);
el_seg=1:L;
h_par=2:nfft/2+1;
sx=zeros(length(f),1);
sy=zeros(length(f),1);
sxy=zeros(length(f),1);
syx=zeros(length(f),1);
x=datawrap(x,L*seg_num);
y=datawrap(y,L*seg_num);
end
if true
  for i=1:seg_num;    
  xx=w.*detrend(x(el_seg));
  yy=w.*detrend(y(el_seg));   
  cx=fft(xx,nfft); 
  cy=fft(yy,nfft);
  mxi=cx(h_par); 
  myi=cy(h_par); 
  sx=sx+mxi.*conj(mxi);
  sy=sy+myi.*conj(myi);
  sxy=sxy+mxi.*conj(myi);
  syx=syx+myi.*conj(mxi);
  el_seg=el_seg+L-ov_num;   
end
end
if true
sx=2*sx/seg_num/R;
sy=2*sy/seg_num/R;
sxy=2*sxy/seg_num/R;
syx=2*syx/seg_num/R;
psdxx=sx*dt;
psdyy=sy*dt;
psdxy=sxy*dt;
psdyx=syx*dt;
frf_mag=abs(psdyy./psdyx);
frf_phase=angle(psdyy./psdyx)*180/pi;
end