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