%% Coherence of MIMO System
% Generate two multichannel signals, each sampled at 1 kHz for 2 seconds. The first signal, the input, consists of three sinusoids with frequencies of 120 Hz, 360 Hz, and 480 Hz. The second signal, the output, is composed of two sinusoids with frequencies of 120 Hz and 360 Hz. One of the sinusoids lags the first signal by π/2. The other sinusoid has a lag of π/4. Both signals are embedded in white Gaussian noise.

fs = 1000;
f = 120;
t = (0:1/fs:2-1/fs)';

inpt = sin(2*pi*f*[1 3 4].*t);
inpt = inpt+randn(size(inpt));
oupt = sin(2*pi*f*[1 3].*t-[pi/2 pi/4]);
oupt = oupt+randn(size(oupt));

%
% Estimate the degree of correlation between all the input signals and each of the output channels. Use a Hamming window of length 100 to window the data. mscohere returns one coherence function for each output channel. The coherence functions reach maxima at the frequencies shared by the input and the output.

[Cxy,f] = mscohere(inpt,oupt,hamming(100),[],[],fs,'mimo');

for k = 1:size(oupt,2)
    subplot(size(oupt,2),1,k)
    plot(f,Cxy(:,k))
    title(['Output ' int2str(k) ', All Inputs'])
end

%
% Switch the input and output signals and compute the multiple coherence function. Use the same Hamming window. There is no correlation between input and output at 480 Hz. Thus there are no peaks in the third correlation function.

[Cxy,f] = mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo');

for k = 1:size(inpt,2)
    subplot(size(inpt,2),1,k)
    plot(f,Cxy(:,k))
    title(['Input ' int2str(k) ', All Outputs'])
end

%
% Repeat the computation, using the plotting functionality of mscohere.

clf
mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo')

%
% Compute the ordinary coherence function of the second signal and the first two channels of the first signal. The off-peak values differ from the multiple coherence function.

[Cxy,f] = mscohere(oupt,inpt(:,[1 2]),hamming(100),[],[],fs);
plot(f,Cxy)

%
% Find the phase differences by computing the angle of the cross-spectrum at the points of maximum coherence.

Pxy = cpsd(oupt,inpt(:,[1 2]),hamming(100),[],[],fs);
[~,mxx] = max(Cxy);
for k = 1:2
    fprintf('Phase lag %d = %5.2f*pi\n',k,angle(Pxy(mxx(k),k))/pi)
end

% Copyright 2012 The MathWorks, Inc.