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)