%% Cross Power Spectral Density of Exponential Sequences
% Generate two N-sample exponential sequences, x_a=a^n and x_b=b^n, with n\ge0. Specify a=0.8, b=0.9, and a small N to see finite-size effects.

N = 10;
n = 0:N-1;

a = 0.8;
b = 0.9;

xa = a.^n;
xb = b.^n;

%
% Compute and plot the cross power spectral density of the sequences over the complete interval of normalized frequencies, [-\pi,\pi]. Specify a rectangular window of length N and no overlap between segments.

w = -pi:1/1000:pi;
wind = rectwin(N);
nove = 0;

[pxx,f] = cpsd(xa,xb,wind,nove,w);

%
% The cross power spectrum of the two sequences has an analytic expression for large N:
% R(\omega)={1\over{1-ae^{-j\omega}}}\,{1\over{1-be^{j\omega}}}.
% Convert this expression to a cross power spectral density by dividing it by 2\pi N. Compare the results. The ripple in the cpsd result is a consequence of windowing.

nfac = 2*pi*N;

X = 1./(1-a*exp(-1j*w));
Y = 1./(1-b*exp( 1j*w));
R = X.*Y/nfac;

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

%
% Repeat the calculation with N=25. The curves agree to six figures for N as small as 100.

N = 25;
n = 0:N-1;
xa = a.^n;
xb = b.^n;

wind = rectwin(N);

[pxx,f] = cpsd(xa,xb,wind,nove,w);
R = X.*Y/(2*pi*N);

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

% Copyright 2012 The MathWorks, Inc.