Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/ndeigcsd.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 % NDEIGCSD calculates TFs from ND cross-correlated spectra.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION:
5 %
6 % Calculates TFs or WFs from Ndim cross correlated spectra. Input
7 % elemnts of a cross-spectral matrix in 2D are assumed to be:
8 %
9 % / csd11(f) csd12(f) \
10 % CSD(f) = | |
11 % \ csd21(f) csd22(f) /
12 %
13 %
14 % CALL: h = eigcsd(csd,varargin)
15 %
16 % INPUT:
17 %
18 % csd are the elements of the cross spectra matrix. It is a (n,n,m)
19 % matrix where n is the dimensionality of the system and m is the
20 % number of frequency samples
21 %
22 % Input also the parameters specifying calculation options
23 %
24 % 'OTP' define the output type. Allowed values are 'TF' output the
25 % transfer functions or 'WF' output the whitening filters frequency
26 % responses. Default 'TF'
27 % 'MTD' define the method for the calculation of the csd matrix of a
28 % multichannel system. Admitted values are 'PAP' referring to Papoulis
29 % [1] style calculation in which csd = TF*I*TF' and 'KAY' referring to
30 % Kay [2] style calculation in which csd = conj(TF)*I*TF.'.
31 % Default 'PAP'
32 %
33 % OUTPUT:
34 %
35 % h are the TFs or WFs frequency responses. It is a (n,n,m) matrix in
36 % which n is the dimensionality of the system and m is the number of
37 % frequency samples.
38 %
39 % REFERENCES:
40 %
41 % [1] A. Papoulis, Probability Random Variable and Stochastic Processes,
42 % McGraw-Hill, third edition, 1991.
43 % [2] S. M. Kay, Modern Spectral Estimation, Prentice Hall, 1988.
44 %
45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 % HISTORY: 23-04-2009 L Ferraioli
47 % Creation
48 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 %
50 % VERSION: $Id: ndeigcsd.m,v 1.4 2009/11/06 16:55:51 luigi Exp $
51 %
52 %
53 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
54 function h = ndeigcsd(csd,varargin)
55
56 [l,m,npts] = size(csd);
57 if m~=l
58 error('!!! The first two dimensions of csd must be equal. csd must be a square matrix frequency by frequency.')
59 end
60
61 % Finding parameters
62
63 % default
64 otp = 'TF';
65 mtd = 'PAP';
66
67 if ~isempty(varargin)
68 for j=1:length(varargin)
69 if strcmp(varargin{j},'OTP')
70 otp = upper(varargin{j+1});
71 end
72 if strcmp(varargin{j},'MTD')
73 mtd = upper(varargin{j+1});
74 end
75 end
76 end
77
78 % Finding suppression
79 suppr = ones(l,l);
80 for ii = 2:l
81 k = ones(l,1);
82 for jj = ii-1:-1:1
83 k(jj) = min(sqrt(csd(jj,jj,:)./csd(ii,ii,:)));
84 if k(jj)>=1
85 suppr(jj,ii) = floor(k(jj));
86 else
87 n=0;
88 while k(jj)<1
89 k(jj)=k(jj)*10;
90 n=n+1;
91 end
92 k(jj) = floor(k(jj));
93 suppr(jj,ii) = k(jj)*10^(-n);
94 end
95 end
96 % csuppr(ii) = prod(suppr(:,ii));
97 end
98 csuppr = prod(suppr,2);
99 supmat = diag(csuppr);
100 ssup = supmat*supmat.';
101
102 supmat = rot90(rot90(supmat));
103 isupmat = inv(supmat);
104
105 % Core Calculation
106
107
108 % initializing output dat
109
110 h = ones(l,m,npts);
111
112 for phi = 1:npts
113
114 % Appliing suppression
115
116 PP = csd(:,:,phi);
117 PP = supmat*PP*supmat;
118
119 % [V,D] = eig(PP,ssup);
120 [V,D,U] = svd(PP,0);
121 % [V,D] = eig(PP);
122
123 % Correcting the output of eig
124 % V = fliplr(V);
125 % D = rot90(rot90(D));
126
127 % % Correcting the output of eig
128 % Vp = fliplr(V);
129 % Lp = rot90(rot90(D));
130 %
131 % % Correcting the phase
132 % [a,b] = size(PP);
133 % for ii=1:b
134 % Vp(:,ii) = Vp(:,ii).*(cos(angle(Vp(ii,ii)))-1i*sin(angle(Vp(ii,ii))));
135 % Vp(ii,ii) = real(Vp(ii,ii));
136 % end
137 %
138 % V = Vp;
139 % D = Lp;
140
141
142 % Definition of the transfer functions
143
144 switch otp
145 case 'TF'
146 switch mtd
147 case 'PAP'
148 % HH = ssup*V*sqrt(D);
149 HH = isupmat*V*sqrt(D);
150 % HH = V*sqrt(D);
151 case 'KAY'
152 % HH = conj(ssup*V*sqrt(D));
153 HH = conj(isupmat*V*sqrt(D));
154 % HH = conj(V*sqrt(D));
155 end
156 case 'WF'
157 switch mtd
158 case 'PAP'
159 %HH = inv(ssup*V*sqrt(D));
160 HH = inv(isupmat*V*sqrt(D));
161 case 'KAY'
162 %HH = inv(conj(ssup*V*sqrt(D)));
163 HH = inv(conj(isupmat*V*sqrt(D)));
164 end
165 end
166
167
168 h(:,:,phi) = HH;
169
170 end
171
172
173 end