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

As said, I compare this to the cpsd command, i.e.

[Pxy,F] = cpsd(wdatax,wdatay,window(L),ov_num,nfft,Fs);

The results "coincide" in the middle of the frequency interval but as mentioned earlier, the differ in the sides of the spectrum.

Can anyone give me any advice on what to implement or re-do in my code to obtain comple accordance?

Thanks in advance and kind regards, Martin.

Answer 1

[Martin]

I can see that I messed up with the code insertion. I hope you can "understand" it despite.