%% SISO and MIMO Cross Power Spectral Densities
% Generate two two-channel sinusoids sampled at 1 kHz for 1 second. The channels of the first signal have frequencies of 200 Hz and 300 Hz. The channels of the second signal have frequencies of 300 Hz and 400 Hz. Both signals are embedded in unit-variance white Gaussian noise.

fs = 1e3;
t = (0:1/fs:1-1/fs)';

q = 2*sin(2*pi*[200 300].*t);
q = q+randn(size(q));

r = 2*sin(2*pi*[300 400].*t);
r = r+randn(size(r));

%
% Compute the cross power spectral density of the two signals. Use a 256-sample Bartlett window to divide the signals into segments and window the segments. Specify 128 samples of overlap between adjoining segments and 2048 DFT points. Use the built-in functionality of cpsd to plot the result.

cpsd(q,r,bartlett(256),128,2048,fs)

% By default, cpsd works column-by-column on matrix inputs of the same size. Each channel peaks at the frequencies of the original sinusoids.
%
% Repeat the calculation, but now append 'mimo' to the list of arguments.

cpsd(q,r,bartlett(256),128,2048,fs,'mimo')

% When called with the 'mimo' option, cpsd returns a three-dimensional array containing cross power spectral density estimates for all combinations of input columns. The estimate of the second channel of q and the first channel of r shows an enhanced peak at the common frequency of 300 Hz.
% Copyright 2012 The MathWorks, Inc.