Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_whiten1D.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 % A test script for ao/whiten1D
2 %
3 % M Hewitson 10-11-08
4 %
5 % $Id: test_ao_whiten1D.m,v 1.14 2010/05/05 04:19:42 mauro Exp $
6 %
7
8 %% Make test data
9 fs = 10;
10 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 10000, 'yunits', 'N'));
11
12 % filter
13 pzm = pzmodel(1e-2, {0.01}, {0.1});
14 ft = miir(pzm,plist('fs',fs));
15
16 af = filter(a, ft);
17
18 %%
19
20 rsp = pzm.resp;
21 rsmm = resp(ft,plist('f',logspace(-4,log10(5),100)));
22 iplot(rsp,rsmm)
23
24 %%
25 pl_psd = plist('scale', 'PSD', 'order', 1, 'navs', 16);
26 axx = a.psd(pl_psd);
27 afxx = af.psd(pl_psd);
28
29 iplot(axx,afxx)
30
31 %% Whiten it with no model
32
33 pl = plist(...
34 'model', [], ...
35 'MaxIter', 30, ...
36 'PoleType', 1, ...
37 'MinOrder', 2, ...
38 'MaxOrder', 9, ...
39 'Weights', 2, ...
40 'Plot', false,...
41 'Disp', false,...
42 'MSEVARTOL', 1e-2,...
43 'FITTOL', 1e-1); % tolerance on MSE Value
44
45 aw = whiten1D(af,pl);
46
47 %% Whiten with a model
48
49 mdl = (abs(rsp).^2).*(2/fs); % this corresponds to the theoretical psd of the data
50 mdl.setYunits(axx.yunits);
51
52 pl = plist(...
53 'fs', fs, ...
54 'model', mdl, ...
55 'MaxIter', 50, ...
56 'PoleType', 2, ...
57 'MinOrder', 2, ...
58 'MaxOrder', 9, ...
59 'Weights', 2, ...
60 'Plot', false,...
61 'Disp', false,...
62 'MSEVARTOL', 1e-2,...
63 'FITTOL', 1e-2); % tolerance on MSE Value
64
65 awm = whiten1D(af, pl);
66
67
68 %%
69
70 iplot(aw, awm)
71
72 %%
73
74 awxx = aw.psd(pl_psd);
75 awmxx = awm.psd(pl_psd);
76
77 iplot(axx, afxx, awxx, awmxx);
78
79 %% Testing flatness
80
81 svec = [axx.data.y afxx.data.y awxx.data.y awmxx.data.y];
82 fc = utils.math.spflat(svec);
83
84 %% Working with multiple inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85
86 fs = 10;
87 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 10000, 'yunits', 'V'));
88 a2 = ao(plist('tsfcn', 'randn(size(t)).''', 'fs', fs, 'nsecs', 10000, 'yunits', 'T'));
89
90 % filter
91 pzm = pzmodel(1e-2, {0.01}, {0.1});
92 ft = miir(pzm,plist('fs',fs));
93
94 af1 = filter(a1, ft);
95 af2 = filter(a2, ft);
96
97 pl = plist(...
98 'model', [], ...
99 'MaxIter', 30, ...
100 'PoleType', 2, ...
101 'MinOrder', 2, ...
102 'MaxOrder', 9, ...
103 'Weights', 2, ...
104 'Plot', false,...
105 'Disp', false,...
106 'MSEVARTOL', 1e-2,...
107 'FITTOL', 1e-2); % tolerance on fit residuals spectral flatness
108
109 [aw1,aw2] = whiten1D(af1,af2,pl);
110
111 %%
112
113 iplot(aw1,aw2)
114
115 %% help test
116
117 % Generate white noise
118 fs = 1;
119 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5, 'yunits','m'));
120
121 % filter
122 ft = miir(plist('type','bandpass','order',3,'gain',1,'fc',[0.03 0.1],'fs',fs));
123
124 % coloring white noise
125 af = filter(a, ft);
126
127 % Whitening colored noise
128 pl = plist(...
129 'model', [], ...
130 'MaxIter', 30, ...
131 'PoleType', 2, ...
132 'MinOrder', 9, ...
133 'MaxOrder', 15, ...
134 'Weights', 2, ...
135 'Plot', false,...
136 'Disp', false,...
137 'MSEVARTOL', 1e-1,...
138 'FITTOL', 5e-2); % tolerance on fit residuals spectral flatness
139
140 aw = whiten1D(af,pl);
141
142 % Calculate psd of colored and whitened data
143 afxx = af.psd;
144 awxx = aw.psd;
145
146 % plotting
147 iplot(afxx,awxx)
148
149 % END