fs = 120;

N = 500;

T = 0:1/fs:N/fs;

win_len = N;

noverlap = 0;

nfft = 500;

trials = 1000;

x = randn(trials,N); y = randn(trials,N); r1 = 0.9; r2 = 0.7;

for trial = 1:trials

for t = 4:N

x(trial,t) = 2*r1*cos(40*(1/fs)*2*pi)*x(trial,t-1) - r1^2*x(trial,t-2) + randn(1);

y(trial,t) = -0.356*x(trial,t-1) + 2*r2*cos(10*(1/fs)*2*pi)*y(trial,t-1) + 0.7136*x(trial,t-2) - r2^2*y(trial,t-2) -0.356*x(trial,t-3) + randn(1);

end

end

xA = reshape(x',1,size(x,1)*size(x,2));

yA = reshape(y',1,size(y,1)*size(y,2));

[pxx,fxx] = pwelch(xA,hanning(win_len),noverlap, nfft, fs);

[pyy,fyy] = pwelch(yA,hanning(win_len),noverlap, nfft, fs);

[pxy,fxy] = cpsd(xA,yA,hanning(win_len),noverlap,nfft, fs);

px = fft_manual(x,win_len,fs);

py = fft_manual(y,win_len,fs);

plot(fxx, 10*log10(pxx), fyy, 10*log10(pyy))

figure; plot(fxx,pxx,fyy,pyy)

figure; plot(fxx,px,fyy,py)