Home > chronux_1_15 > spikesort > utility > datatools > mtm_cohere.m

mtm_cohere

PURPOSE ^

CMTM Coherence function estimate using the multitaper method.

SYNOPSIS ^

function [Cxy, F] = cmtm(X,Y,NW,Fs)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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;

Generated on Tue 15-Aug-2006 22:51:57 by m2html © 2003