Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/eigpsd.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 % EIGPSD calculates TFs from 2D cross-correlated spectra.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION:
5 %
6 % Calculates TFs or WFs from 2dim cross correlated spectra
7 %
8 % CALL: [h11,h12,h21,h22] = eigpsd(psd1,csd,psd2,varargin)
9 %
10 % INPUT:
11 %
12 % psd1 is the first power spectral density
13 % csd is the cross spectrum
14 % psd2 is the second power spectral density
15 % Input also the parameters specifying calculation options
16 % 'USESYM' define the calculation method, allowed values are 0
17 % (double precision arithmetic) and 1 (symbolic toolbox variable
18 % precision arithmetic)
19 % 'DIG' define the digits used in VPA calculation
20 % 'OTP' define the output type. Allowed values are 'TF' output the
21 % transfer functions or 'WF' output the whitening filters frequency
22 % responses
23 %
24 % OUTPUT:
25 %
26 % h11, h12, h21 and h22 are the four innovation TFs frequency responses
27 % or the four whitening filters frequency responses
28 %
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30 % HISTORY: 02-10-2008 L Ferraioli
31 % Creation
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 %
34 % VERSION: $Id: eigpsd.m,v 1.3 2008/10/24 14:40:20 luigi Exp $
35 %
36 %
37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38 function [h11,h12,h21,h22] = eigpsd(psd1,csd,psd2,varargin)
39
40 % Init
41
42 S11 = psd1;
43 S22 = psd2;
44 S12 = csd;
45 S21 = conj(S12);
46
47 npts = length(S11);
48
49 % Finding parameters
50
51 % default
52 usesym = 1;
53 dig = 50;
54 otp = 'TF';
55
56 if ~isempty(varargin)
57 for j=1:length(varargin)
58 if strcmp(varargin{j},'USESYM')
59 usesym = varargin{j+1};
60 end
61 if strcmp(varargin{j},'DIG')
62 dig = varargin{j+1};
63 end
64 if strcmp(varargin{j},'OTP')
65 otp = varargin{j+1};
66 end
67 end
68 end
69
70 % Defining calculation method
71 if usesym == 1
72 method = 'VPA';
73 else
74 method = 'NUM';
75 end
76
77 % Finding suppressio
78
79 k = min(sqrt(S11./S22));
80 if k>=1
81 suppr = floor(k);
82 else
83 n=0;
84 while k<1
85 k=k*10;
86 n=n+1;
87 end
88 k = floor(k);
89 suppr = k*10^(-n);
90 end
91
92 supmat = [1 0;0 suppr];
93 isupmat = [1 0;0 1/suppr];
94
95 % Core Calculation
96
97 switch method
98 case 'NUM'
99
100 % initializing output data
101 h11 = ones(npts,1);
102 h12 = ones(npts,1);
103 h21 = ones(npts,1);
104 h22 = ones(npts,1);
105
106 for phi = 1:npts
107
108 % Appliing suppression
109 PP = supmat*[S11(phi) S12(phi);S21(phi) S22(phi)]*supmat;
110
111 % Calculate eignevalues Matrix Lp
112 k = (4*PP(1,2)*PP(2,1))/((PP(1,1)-PP(2,2))^2);
113
114 Lp1 = (PP(1,1)+PP(2,2)+(PP(1,1)-PP(2,2))*sqrt(1+k))/2;
115 Lp2 = (PP(1,1)+PP(2,2)-(PP(1,1)-PP(2,2))*sqrt(1+k))/2;
116
117 % Calculate eigenvectors Matrix Vp, eigenvectors are on the columns
118 Vp1 = (-1/sqrt(1+(k/((1+sqrt(1+k))^2))))*[1;2*PP(2,1)/((PP(1,1)-PP(2,2))*(1+sqrt(1+k)))];
119 Vp2 = (1/sqrt(1+(((1-sqrt(1+k))^2)/k)))*[(PP(1,1)-PP(2,2))*(1-sqrt(1+k))/(2*PP(2,1));1];
120
121 % Definition of the transfer functions
122 switch otp
123 case 'TF'
124 HH = isupmat*[conj(Vp1) conj(Vp2)]*[sqrt(Lp1) 0;0 sqrt(Lp2)];
125 case 'WF'
126 HH = [1/sqrt(Lp1) 0;0 1/sqrt(Lp2)]*inv([conj(Vp1) conj(Vp2)])*supmat;
127 end
128
129 h11(phi,1) = HH(1,1);
130 h12(phi,1) = HH(1,2);
131 h21(phi,1) = HH(2,1);
132 h22(phi,1) = HH(2,2);
133
134 end
135
136 case 'VPA'
137 % Define the numerical precision
138 digits(dig)
139
140 % initializing output data
141 h11 = vpa(ones(npts,1));
142 h12 = vpa(ones(npts,1));
143 h21 = vpa(ones(npts,1));
144 h22 = vpa(ones(npts,1));
145
146 SS11 = vpa(S11);
147 SS12 = vpa(S12);
148 SS21 = vpa(S21);
149 SS22 = vpa(S22);
150
151 for phi = 1:npts
152
153 % Appliing suppression
154 PP = supmat*[SS11(phi) SS12(phi);SS21(phi) SS22(phi)]*supmat;
155
156 % Calculate eignevalues Matrix Lp
157 k = (4*PP(1,2)*PP(2,1))/((PP(1,1)-PP(2,2))^2);
158
159 Lp1 = (PP(1,1)+PP(2,2)+(PP(1,1)-PP(2,2))*sqrt(1+k))/2;
160 Lp2 = (PP(1,1)+PP(2,2)-(PP(1,1)-PP(2,2))*sqrt(1+k))/2;
161
162 % Calculate eigenvectors Matrix Vp, eigenvectors are on the columns
163 Vp1 = (-1/sqrt(1+(k/((1+sqrt(1+k))^2))))*[1;2*PP(2,1)/((PP(1,1)-PP(2,2))*(1+sqrt(1+k)))];
164 Vp2 = (1/sqrt(1+(((1-sqrt(1+k))^2)/k)))*[(PP(1,1)-PP(2,2))*(1-sqrt(1+k))/(2*PP(2,1));1];
165
166 % Definition of the transfer functions
167 switch otp
168 case 'TF'
169 HH = isupmat*[conj(Vp1) conj(Vp2)]*[sqrt(Lp1) 0;0 sqrt(Lp2)];
170 case 'WF'
171 HH = [1/sqrt(Lp1) 0;0 1/sqrt(Lp2)]*inv([conj(Vp1) conj(Vp2)])*supmat;
172 end
173
174 h11(phi,1) = HH(1,1);
175 h12(phi,1) = HH(1,2);
176 h21(phi,1) = HH(2,1);
177 h22(phi,1) = HH(2,2);
178
179 end
180 h11 = double(h11);
181 h12 = double(h12);
182 h21 = double(h21);
183 h22 = double(h22);
184 end