Mercurial > hg > ltpda

comparison m-toolbox/test/test_matrix_MultiChannelNoise.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % TEST matrix/MultiChannelNoise
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 23-04-2009 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION
9
10 '$Id: test_matrix_fromCSD.m,v 1.3 2009/04/23 16:03:55 luigi Exp $';
11
12 %% Loading spectra
13
14 load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
15
16 f = mpsd(:,1);
17 psd = mpsd(:,2);
18 fs = 10;
19
20 % 1dim model
21 mod1D = ao(plist('xvals', f, 'yvals', psd, 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1'));
22 mod1D.setName;
23
24 % 2dim model
25 csd11 = ao(plist('xvals', f, 'yvals', mpsd(:,2), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1 PSD'));
26 csd11.setName;
27 csd12 = ao(plist('xvals', f, 'yvals', mpsd(:,3), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH12 CSD'));
28 csd12.setName;
29 csd21 = ao(plist('xvals', f, 'yvals', conj(mpsd(:,3)), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH21 CSD'));
30 csd21.setName;
31 csd22 = ao(plist('xvals', f, 'yvals', mpsd(:,4), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH2 PSD'));
32 csd22.setName;
33
34 mod2D = [csd11 csd12;csd21 csd22];
35
36 %%
37 %--------------------------------------------------------------------------
38 % 1 Dimensional Noisegen Filter
39 %--------------------------------------------------------------------------
40
41 %% 1 dim noisegen filter in s domain - symbolic precision All pass
42
43 % plist for noise generation
44 plns = plist(...
45 'csd', mod1D, ...
46 'targetobj', 'parfrac', ...
47 'Nsecs', 1e4, ...
48 'fs', fs, ...
49 'MaxIter', 70, ...
50 'PoleType', 3, ...
51 'MinOrder', 15, ...
52 'MaxOrder', 30, ...
53 'Weights', 2, ...
54 'Plot', true,...
55 'MSEVARTOL', 1e-2,...
56 'FITTOL', 1e-3,...
57 'UseSym', 'on');
58
59 na = matrix(plns);
60
61
62
63 %% 1 dim noisegen filter in z domain - symbolic precision All pass
64
65 % plist for noise generation
66 plns = plist(...
67 'csd', mod1D, ...
68 'targetobj', 'miir', ...
69 'Nsecs', 1e4, ...
70 'fs', fs, ...
71 'MaxIter', 70, ...
72 'PoleType', 3, ...
73 'MinOrder', 10, ...
74 'MaxOrder', 30, ...
75 'Weights', 2, ...
76 'Plot', true,...
77 'MSEVARTOL', 1e-2,...
78 'FITTOL', 1e-3,...
79 'UseSym', 'on');
80
81 na = matrix(plns);
82
83 %%
84 %--------------------------------------------------------------------------
85 % 2 Dimensional Noisegen Filter
86 %--------------------------------------------------------------------------
87
88
89 %% 2 dim noisegen filter in s domain - symbolic precision All pass
90
91 % plist for noise generation
92 plns = plist(...
93 'csd', mod2D, ...
94 'targetobj', 'parfrac', ...
95 'MaxIter', 70, ...
96 'PoleType', 3, ...
97 'MinOrder', 25, ...
98 'MaxOrder', 40, ...
99 'Weights', 2, ...
100 'Plot', false,...
101 'MSEVARTOL', 1e-2,...
102 'FITTOL', 1e-4,...
103 'UseSym', 'on');
104
105 na = matrix(plns);
106
107 %% Check response
108
109 na11 = resp(na.objs(1,1),plist('f',f,'bank','parallel'));
110 na12 = resp(na.objs(1,2),plist('f',f,'bank','parallel'));
111 na21 = resp(na.objs(2,1),plist('f',f,'bank','parallel'));
112 na22 = resp(na.objs(2,2),plist('f',f,'bank','parallel'));
113
114 mna = matrix([na11 na12;na21 na22]);
115
116 ECSD = mna*conj(transpose(mna));
117
118 iplot(csd11,ECSD.objs(1,1))
119 iplot(csd12,ECSD.objs(1,2))
120 iplot(csd22,ECSD.objs(2,2))
121
122
123 %% 2 dim noisegen filter in z domain - symbolic precision All pass
124
125 % plist for noise generation
126 plns = plist(...
127 'csd', mod2D, ...
128 'targetobj', 'miir', ...
129 'fs', fs, ...
130 'MaxIter', 70, ...
131 'PoleType', 3, ...
132 'MinOrder', 25, ...
133 'MaxOrder', 40, ...
134 'Weights', 2, ...
135 'Plot', false,...
136 'MSEVARTOL', 1e-1,...
137 'FITTOL', 1e-4,...
138 'UseSym', 'on');
139
140 na = matrix(plns);
141
142 %% Check response
143
144 na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel'));
145 na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel'));
146 na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel'));
147 na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel'));
148
149 mna = matrix([na11 na12;na21 na22]);
150
151 ECSD = mna*conj(transpose(mna));
152
153 iplot(csd11,ECSD.objs(1,1))
154 iplot(csd12,ECSD.objs(1,2))
155 iplot(csd22,ECSD.objs(2,2))
156
157 %%
158 %--------------------------------------------------------------------------
159 % 2 Dimensional Noise Generation
160 %--------------------------------------------------------------------------
161
162 plng = plist('model',na,...
163 'Nsecs',2e5,...
164 'fs',fs);
165
166 no = matrix(plng);
167
168 %% Check psd
169
170 no1 = no.objs(1);
171 no2 = no.objs(2);
172
173 no1xx = no1.psd;
174 no2xx = no2.psd;
175
176 iplot(no1xx,csd11,no2xx,csd22)
177
178 %% Check cross-coherence
179
180 coh = csd12./sqrt(csd11.*csd22);
181 rcoh = real(coh);
182 icoh = imag(coh);
183
184 no12xx = cpsd(no1,no2,plist('nfft',1e4));
185 no1xx = no1.psd(plist('nfft',1e4));
186 no2xx = no2.psd(plist('nfft',1e4));
187
188 ecoh = no12xx./(sqrt(no1xx.*no2xx));
189 recoh = real(ecoh);
190 iecoh = imag(ecoh);
191
192 iplot(recoh,rcoh,plist('Yscales',{'All','lin'}))
193 iplot(iecoh,icoh,plist('Yscales',{'All','lin'}))
194
195
196
197