Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/mltfe.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % MLTFE compute log-frequency space TF
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: MLTFE compute log-frequency space TF
5 %
6 % CALL: Txy = mltfe(X,f,r,m,L,fs,win,order,olap,Lmin,method,variance)
7 %
8 % VERSION: $Id: mltfe.m,v 1.25 2009/11/27 15:18:45 miquel Exp $
9 %
10 % HISTORY: 19-05-2007 M Hewitson
11 % Creation
12 %
13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14
15 function varargout = mltfe(varargin)
16
17 import utils.const.*
18
19 % Get inputs
20 X = varargin{1};
21 f = varargin{2};
22 r = varargin{3};
23 m = varargin{4};
24 L = varargin{5};
25 K = varargin{6};
26 fs = varargin{7};
27 win = varargin{8};
28 order = varargin{9};
29 olap = varargin{10};
30 Lmin = varargin{11};
31 method = varargin{12};
32
33 % --- Prepare some variables
34 si = size(X);
35 nc = si(1);
36 nf = length(f);
37 Txy = zeros(nf,1);
38 dev = zeros(nf,1);
39 disp_each = round(nf/100)*10;
40
41 winType = win.type;
42 winPsll = win.psll;
43
44 % ----- Loop over Frequency
45 for fi=1:nf
46 [Txy(fi) dev(fi)]= computeTF(fs, L(fi), K(fi), m, winType, winPsll, X, olap, order, nc, f(fi), fi, nf, disp_each, method);
47 end
48
49 % Set output
50 varargout{1} = Txy;
51 varargout{2} = dev;
52 end
53
54 %--------------------------------------------------------------------------
55 % Function to run over channels
56 function [Txy,dev]= computeTF(fs, l, K, m, winType, winPsll, X, olap, order, nc, ffi, fi, nf, disp_each, method)
57
58 switch lower(winType)
59 case 'kaiser'
60 lwin = specwin(winType, l, winPsll);
61 otherwise
62 lwin = specwin(winType, l);
63 end
64
65 % Compute DFT coefficients
66 twopi = 2.0*pi;
67 p = 1i * twopi * m(fi)/l.*[0:l-1];
68 C = lwin.win .* exp(p);
69
70 if mod(fi,disp_each) == 0 || fi == 1 || fi == nf
71 utils.helper.msg(utils.const.msg.PROC1, 'computing frequency %04d of %04d: %f Hz', fi, nf, ffi);
72 end
73 % Loop over input channels
74 [Txy,dev] = in2out(l, X, K, olap, order, nc, method, C, lwin, fs);
75
76 end
77
78 %--------------------------------------------------------------------------
79 % Compute 1 input to multiple outputs
80 function [Txy,dev]= in2out(l, X, K, olap, order, nc, method, C, lwin, fs)
81
82 % if no errors are required the function returns zero but errors are not
83 % stored in the final ao
84 dev = 0;
85
86 switch lower(method)
87 case 'tfe'
88 % Core cross-DFT part in C-mex file
89 % We need cross-spectrum and Power spectrum
90 [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
91 Txy = conj(XY)/(XX);
92 if nsegs == 1
93 dev = Inf;
94 else
95 dev = sqrt((nsegs/(nsegs-1)^2)*(YY./XX).*(1 - (abs(XY).^2)./(XX.*YY)));
96 % dP = sqrt((k/(k-1)^2)*(Pyy./Pxx).*(1 - (abs(Pxy).^2)./(Pxx.*Pyy)));
97 end
98
99 case 'cpsd'
100 % Core cross-DFT part in C-mex file
101 [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
102 S1 = lwin.ws;
103 S2 = lwin.ws2;
104 Txy = 2.0*XY/fs/S2;
105 Var = 4.0*M2/fs^2/S2^2/nsegs;
106 if nsegs == 1
107 dev = Inf;
108 else
109 dev = sqrt(Var);
110 end
111
112 case 'mscohere'
113 % Core cross-DFT part in C-mex file
114 % We need cross-spectrum and Power spectrum
115 [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
116 Txy = (abs(XY).^2)./(XX.*YY); % Magnitude-squared coherence
117 if nsegs == 1
118 dev = Inf;
119 else
120 dev = sqrt((2*Txy/nsegs).*(1-Txy).^2);
121 end
122
123 case 'cohere'
124 % Core cross-DFT part in C-mex file
125 % We need cross-spectrum and Power spectrum
126 [XY, XX, YY, M2, nsegs] = ltpda_dft(X(1,:), X(2,:), l, C, olap, order);
127 Txy = XY./sqrt(XX.*YY); % Complex coherence
128 if nsegs == 1
129 dev = Inf;
130 else
131 dev = sqrt((2*abs(Txy)/nsegs).*(1-abs(Txy)).^2);
132 end
133 end
134
135 end