x=maledesire(1:30,1);
y=femaledesire(1:30,1);
N = length(x);
nfft=N/2;
Fs=1;
window=rectwin(18);
lfft=2*nfft;
overlap=12; %12
winsize = length(window);
numwin=((N-winsize)/(winsize-overlap))+1; %3
stride = winsize - overlap; %6
iter = numwin; %3

istart = 1;
iend = istart + winsize - 1; %18


for i=1:iter
    %apply window, perform FFT, and Normalize
    d1m=(x(istart:iend,1).*window); %(18x1)
    d1f=(y(istart:iend,1).*window); %(18x1)
    xdft1 = fft(d1m,lfft);
    ydft1 = fft(d1f,lfft);
    xdftnorm= xdft1;
    ydftnorm= ydft1;

    % calculate CPSD
    Sxy1a=xdftnorm.*conj(ydftnorm);
    Sxy1b=Sxy1a./N^2;
    Sxy1c=Sxy1b(1:lfft/2+1);
    Sxy1c(2:end-1)=2*Sxy1c(2:end-1);
    Sxy1d=abs(Sxy1c)';
    Sxycum(1:lfft/2+1,i)=Sxy1d';

    %move to next window segment
istart = istart + stride;
iend = istart + (winsize-1);
end

Sxy5=mean(Sxycum,2);

Compare Sxy5 with Sxy7:

Sxy6=cpsd(x,y,window,overlap,lfft,Fs);
Sxy7=abs(Sxy6);