nicmatlabscripts/coherence.m at master · nicolassmith/nicmatlabscripts · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
function [A]= coherence(dataX, dataY, Fs, resolution, windowFunction, ...
   OverlapPercent, verbose)
%===============================================
% Function A = coherence(dataX, dataY, Fs, resolution, windowFunction, ...
%   OverlapPercent, verbose)
%
%  computes the coherence like mscohere
%  default rate = 16384
%  default_resolution = 1;
%  default_window     = @hanning; <- this is a function handle
%  default_overlap    = 50;
%  default_verbose    = 0;
%
% Based on asd by L. Cadonati.
% Nicolas Smith-Lefebvre 2011
%
% Ex 1:  A = coherence(data, 16384, 0.125)
%===============================================

narginchk(2, 7);
%====================================
default_rate = 16384;
default_resolution = 1;
default_window     = @hanning;
default_overlap    = 50;
default_verbose    = 0;
%====================================
if nargin == 2,
  Fs = default_rate;
  resolution = default_resolution;
  windowFunction = default_window;
  OverlapPercent = default_overlap;
  verbose = default_verbose;

elseif nargin == 3,
  resolution = default_resolution;
  windowFunction = default_window;
  OverlapPercent = default_overlap;
  verbose = default_verbose;
elseif nargin == 4,
  windowFunction = default_window;
  OverlapPercent = default_overlap;
  verbose = default_verbose;
elseif nargin == 5,
  OverlapPercent = default_overlap;
  verbose = default_verbose;
elseif nargin == 6,
  verbose = default_verbose;
end
%====================================


% force row vector
dataX = dataX(:).';
dataY = dataY(:).';

% segment lengths
SegmentDuration = 1/resolution;
SegmentLength   = Fs*SegmentDuration;

noverlap = OverlapPercent*SegmentLength/100;

% ensure integer power of two segment length
log2SegmentLength = log(SegmentLength) / log(2);
if log2SegmentLength ~= floor(log2SegmentLength),
  error('rate / resolution must be an integer power of two');
end

if (verbose),
  fprintf('duration       = %g\n',SegmentDuration);
  fprintf('segmentLength  = %g\n',SegmentLength);
end

%h = spectrum.welch(windowFunction,SegmentLength,OverlapPercent);
%Hpsd = psd(h,data,'Fs',Fs);
[Cxy,F] = mscohere(dataX,dataY,windowFunction(SegmentLength),noverlap,SegmentLength,Fs);


A = struct('Name','Coherence',...
           'Cxy',Cxy,'f',F);