How can I calculate the cross power spectral density from fft data?

Hello, I need to calculate the cross power spectrum of two signals

I have calculated the fourrier transform of them by using the following formula. Could you help me to write a code to calculate and plot cross power spectrum?

X=DATA1;

xdft = fft(X(:,2));

% sampling interval -- assuming equal sampling

DT = X(2,1)-X(1,1);

% sampling frequency

Fs = 1/DT;

DF = Fs/size(X,1);

freq = 0:DF:Fs/2;

xdft = xdft(1:round(length(X)/2));

%plot(freq,abs(xdft))

Y=DATA2;

xdft2 = fft(Y(:,2));

% sampling interval -- assuming equal sampling

DT2 = Y(2,1)-Y(1,1);

% sampling frequency

Fs2 = 1/DT2;

DF2 = Fs2/size(Y,1);

freq2 = 0:DF2:Fs2/2;

xdft2 = xdft2(1:round(length(X)/2));

plot(freq,abs(xdft2))