Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/eigcsd.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 % EIGCSD calculates TFs from 2D cross-correlated spectra.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION:
5 %
6 % Calculates TFs or WFs from 2dim cross correlated spectra. Input the
7 % elemnts of a cross-spectral matrix
8 %
9 % / csd11(f) csd12(f) \
10 % CSD(f) = | |
11 % \ csd21(f) csd22(f) /
12 %
13 % and output the frequency response of four innovation transfer
14 % functions or four whitening filters that can be used to color white
15 % noise or to whitening colored noise.
16 %
17 % CALL: [h11,h12,h21,h22] = eigcsd(csd11,csd12,csd21,csd22,varargin)
18 % [h11,h12,h21,h22] = eigcsd(csd11,csd12,[],csd22,varargin)
19 % [h11,h12,h21,h22] = eigcsd(csd11,[],csd21,csd22,varargin)
20 %
21 % INPUT:
22 %
23 % csd11, csd12, csd21 and csd22 are the elements of the cross spectral
24 % matrix
25 % Input also the parameters specifying calculation options
26 % 'USESYM' define the calculation method, allowed values are 0
27 % (double precision arithmetic) and 1 (symbolic toolbox variable
28 % precision arithmetic)
29 % 'DIG' define the digits used in VPA calculation
30 % 'OTP' define the output type. Allowed values are 'TF' output the
31 % transfer functions or 'WF' output the whitening filters frequency
32 % responses
33 % 'KEEPVAR' get a whitening filter that preserve the variance of
34 % the input data. Values are true or false
35 % 'VARS' first and second channels variance
36 %
37 % OUTPUT:
38 %
39 % h11, h12, h21 and h22 are the four innovation TFs frequency responses
40 % or the four whitening filters frequency responses
41 %
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 % HISTORY: 02-10-2008 L Ferraioli
44 % Creation
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 %
47 % VERSION: $Id: eigcsd.m,v 1.8 2009/06/16 16:06:50 luigi Exp $
48 %
49 %
50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 function [h11,h12,h21,h22] = eigcsd(csd11,csd12,csd21,csd22,varargin)
52
53 % Init
54
55 S11 = csd11;
56 S22 = csd22;
57
58 if isempty(csd12) && isempty(csd21)
59 error(' You must input csd12 or csd21 ')
60 end
61
62 if isempty(csd12)
63 S12 = conj(csd21);
64 else
65 S12 = csd12;
66 end
67 if isempty(csd21)
68 S21 = conj(csd12);
69 else
70 S21 = csd21;
71 end
72
73 npts = length(S11);
74
75 % Finding parameters
76
77 % default
78 usesym = 0;
79 dig = 50;
80 otp = 'TF';
81 kv = false;
82 vars = [1 1];
83 % fs = 10;
84
85 if ~isempty(varargin)
86 for j=1:length(varargin)
87 if strcmp(varargin{j},'USESYM')
88 usesym = varargin{j+1};
89 end
90 if strcmp(varargin{j},'DIG')
91 dig = varargin{j+1};
92 end
93 if strcmp(varargin{j},'OTP')
94 otp = varargin{j+1};
95 end
96 if strcmp(varargin{j},'KEEPVAR')
97 kv = varargin{j+1};
98 end
99 if strcmp(varargin{j},'VARS')
100 vars = varargin{j+1};
101 end
102 end
103 end
104
105 % Defining calculation method
106 if usesym == 1
107 method = 'VPA';
108 else
109 method = 'NUM';
110 end
111
112 % Finding suppression
113
114 k = min(sqrt(S11./S22));
115 if k>=1
116 suppr = floor(k);
117 else
118 n=0;
119 while k<1
120 k=k*10;
121 n=n+1;
122 end
123 k = floor(k);
124 suppr = k*10^(-n);
125 end
126
127 supmat = [1 0;0 suppr];
128 % isupmat = [1 0;0 1/suppr];
129 isupmat = inv(supmat);
130
131 % check for keep variance option
132 if isequal(otp,'WF') && kv
133 % get mean of spectra
134 s1 = sqrt(vars(1));
135 s2 = sqrt(vars(2));
136 end
137
138 % Core Calculation
139
140 switch method
141 case 'NUM'
142
143 % initializing output data
144 h11 = ones(npts,1);
145 h12 = ones(npts,1);
146 h21 = ones(npts,1);
147 h22 = ones(npts,1);
148
149 % T = zeros(2,2,npts);
150 % D = zeros(2,2,npts);
151
152 for phi = 1:npts
153
154 % Appliing suppression
155 PP = supmat*[S11(phi) S12(phi);S21(phi) S22(phi)]*supmat;
156
157 [V,D] = eig(PP);
158
159 % Correcting the output of eig
160 Vp = fliplr(V);
161 Lp = rot90(rot90(D));
162
163 % Correcting the phase
164 [a,b] = size(PP);
165 for ii=1:b
166 Vp(:,ii) = Vp(:,ii).*(cos(angle(Vp(ii,ii)))-1i*sin(angle(Vp(ii,ii))));
167 Vp(ii,ii) = real(Vp(ii,ii));
168 end
169
170 % Definition of the transfer functions
171 switch otp
172 case 'TF'
173 HH = isupmat*[Vp(:,1) Vp(:,2)]*[sqrt(Lp(1,1)) 0;0 sqrt(Lp(2,2))];
174 case 'WF'
175 HH = [1/sqrt(Lp(1,1)) 0;0 1/sqrt(Lp(2,2))]*inv(Vp)*supmat;
176 % HH = [1/sqrt(Lp(1,1)) 0;0 1/sqrt(Lp(2,2))]*(Vp')*supmat;
177 end
178
179 if isequal(otp,'WF') && kv
180 h11(phi,1) = s1*HH(1,1);
181 h12(phi,1) = s1*HH(1,2);
182 h21(phi,1) = s2*HH(2,1);
183 h22(phi,1) = s2*HH(2,2);
184 else
185 h11(phi,1) = HH(1,1);
186 h12(phi,1) = HH(1,2);
187 h21(phi,1) = HH(2,1);
188 h22(phi,1) = HH(2,2);
189 end
190
191 % T(:,:,phi) = conj(HH)*[S11(phi) S12(phi);S21(phi) S22(phi)]*HH.';
192 % D(:,:,phi) = conj(HH)*HH.' - [S11(phi) S12(phi);S21(phi) S22(phi)];
193
194 end
195
196 case 'VPA'
197 % Define the numerical precision
198 digits(dig)
199
200 % initializing output data
201 h11 = vpa(ones(npts,1));
202 h12 = vpa(ones(npts,1));
203 h21 = vpa(ones(npts,1));
204 h22 = vpa(ones(npts,1));
205
206 SS11 = vpa(S11);
207 SS12 = vpa(S12);
208 SS21 = vpa(S21);
209 SS22 = vpa(S22);
210
211 for phi = 1:npts
212
213 % Appliing suppression
214 PP = supmat*[SS11(phi) SS12(phi);SS21(phi) SS22(phi)]*supmat;
215
216 [V,D] = eig(PP);
217 Vp1 = V(:,1);
218 Vp2 = -1.*V(:,2);
219 Lp1 = D(1,1);
220 Lp2 = D(2,2);
221
222 % Definition of the transfer functions
223 switch otp
224 case 'TF'
225 HH = isupmat*[Vp1 Vp2]*[sqrt(Lp1) 0;0 sqrt(Lp2)];
226 case 'WF'
227 % HH = [1/sqrt(Lp1) 0;0 1/sqrt(Lp2)]*inv([Vp1 Vp2])*supmat;
228 HH = inv(isupmat*[Vp1 Vp2]*[sqrt(Lp1) 0;0 sqrt(Lp2)]);
229 end
230
231 if isequal(otp,'WF') && kv
232 h11(phi,1) = s1*HH(1,1);
233 h12(phi,1) = s1*HH(1,2);
234 h21(phi,1) = s2*HH(2,1);
235 h22(phi,1) = s2*HH(2,2);
236 else
237 h11(phi,1) = HH(1,1);
238 h12(phi,1) = HH(1,2);
239 h21(phi,1) = HH(2,1);
240 h22(phi,1) = HH(2,2);
241 end
242
243 end
244 h11 = double(h11);
245 h12 = double(h12);
246 h21 = double(h21);
247 h22 = double(h22);
248 end
249