Mercurial > hg > ltpda

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