How to make a "cross power spectral density" based on AR model (or BURG method)

the matlab function 'cpsd' can make a cpsd estimate, but the cpsd curve is not smooth enough and has bias error.It is a poor estimate because the 'cpsd' function is based an WELCH method --- a periodogram method.

[Pxx,f]=cpsd(x,y,window,noverlap,nfft,fs);[Pxx,f]=pwelch(xn,window,noverlap,nfft,fs,)

So I'm trying to make a cpsd estimate based on AR model, because the psd curve is more smooth and more accurate by using a PBURG function [Pxx,f] = pburg(x,p,nfft,fs).

The problem is: the psd estimate is available --- pburg function, how to make a cpsd estimate based on burg method??