


CMTM Coherence function estimate using the multitaper method.
Cxy = CMTM(X,Y) estimates the coherence between two equal length
vectors vectors X and Y using Thomson's multitaper method. The
coherence is a complex function of frequency, with magnitude between 0
and 1, that estimates
E[X*(f) Y(f)] / sqrt(E[X*(f) X(f)] E[Y*(f) Y(f)])
where * denotes complex conjugation. **Note that this is different
from the matlab function COHERE, which returns the magnitude squared
of this value.**
Cxy = CMTM(X,Y,NW) specifies the time-bandwidth product for the
discrete prolate spheroidal sequences (DPSS) is specified by NW; this
value also determines the number of tapers as (2*NW-1). If not
given, the default NW is 4.
Cxy = CMTM(X,Y,NW,Fs) specifies a sampling frequency Fs.
[Cxy,F] = CMTM(...) also returns the vector of frequencies at which
the coherence is estimated. If Fs is specified, this vector ranges
from [0,Fs/2] and the units are the same as Fs. If Fs is not
specified, F ranges from [0,1/2] and the units are relative to the
sampling frequency.
CMTM(...) without output arguments plots the magnitude-squared
and phase of the coherence (in two subplots) in the current figure.

0001 function [Cxy, F] = cmtm(X,Y,NW,Fs) 0002 %CMTM Coherence function estimate using the multitaper method. 0003 % Cxy = CMTM(X,Y) estimates the coherence between two equal length 0004 % vectors vectors X and Y using Thomson's multitaper method. The 0005 % coherence is a complex function of frequency, with magnitude between 0 0006 % and 1, that estimates 0007 % E[X*(f) Y(f)] / sqrt(E[X*(f) X(f)] E[Y*(f) Y(f)]) 0008 % where * denotes complex conjugation. **Note that this is different 0009 % from the matlab function COHERE, which returns the magnitude squared 0010 % of this value.** 0011 % 0012 % Cxy = CMTM(X,Y,NW) specifies the time-bandwidth product for the 0013 % discrete prolate spheroidal sequences (DPSS) is specified by NW; this 0014 % value also determines the number of tapers as (2*NW-1). If not 0015 % given, the default NW is 4. 0016 % 0017 % Cxy = CMTM(X,Y,NW,Fs) specifies a sampling frequency Fs. 0018 % 0019 % [Cxy,F] = CMTM(...) also returns the vector of frequencies at which 0020 % the coherence is estimated. If Fs is specified, this vector ranges 0021 % from [0,Fs/2] and the units are the same as Fs. If Fs is not 0022 % specified, F ranges from [0,1/2] and the units are relative to the 0023 % sampling frequency. 0024 % 0025 % CMTM(...) without output arguments plots the magnitude-squared 0026 % and phase of the coherence (in two subplots) in the current figure. 0027 0028 %%%%% Debugging test code: 0029 % Fs = 100; 0030 % t = 0:1/Fs:(32-1/Fs); 0031 % X = sin(2*pi*5*t); Y = sin(2*pi*5*t) + 0.2*randn(size(X)); 0032 % cmtm(X,Y,4,Fs); 0033 0034 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Check Inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0035 if (nargin < 2), error('Three arguments are required.'); end 0036 if (nargin < 3), NW = 4; end; 0037 if (nargin < 4), Fs = 1; end; 0038 0039 if (size(X,1) > 1), X = X'; end; 0040 if (size(Y,1) > 1), Y = Y'; end; 0041 [Mx,Nx] = size(X); [My,Ny] = size(Y); 0042 if (Mx~=1 || My~=1), error('Matrix arguments are not supported.'); end; 0043 if (Nx ~= Ny), error('A Coherence estimate requires two vectors of equal length.'); end; 0044 if (length(NW) ~= 1), error('NW must be a scalar.'); end; 0045 if (length(Fs) ~= 1), error('Fs must be a scalar.'); end; 0046 0047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% Make Tapers %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0048 P = 2*NW-1; 0049 [E,V] = dpss(Nx,NW,P); 0050 E = E'; 0051 0052 %%%%%%%%%%%%%%%%%%%%%%%%%% Calculate Coherence %%%%%%%%%%%%%%%%%%%%%%% 0053 midway = floor(Nx/2)+1; 0054 Xf = fft(E .* repmat(X,P,1),[],2); Xf = Xf(:,1:midway); Pxx = conj(Xf).*Xf; 0055 Yf = fft(E .* repmat(Y,P,1),[],2); Yf = Yf(:,1:midway); Pyy = conj(Yf).*Yf; 0056 Cxy = (conj(Xf) .* Yf) ./ sqrt((Pxx.*Pyy)); 0057 Cxy = mean(Cxy,1); 0058 0059 0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%% Plot if Needed %%%%%%%%%%%%%%%%%%%%%%%%%% 0061 F = (0:(1/Nx):0.5)*Fs; 0062 0063 if (nargout == 0) 0064 subplot(2,1,1); plot(F,abs(Cxy)); ylabel('Magnitude Squared Coherence'); 0065 subplot(2,1,2); plot(F,angle(Cxy)); ylabel('Coherence Phase'); xlabel('Frequency'); 0066 0067 clear Cxy 0068 end 0069 0070 if (nargout < 2) 0071 clear F; 0072 end;