comparison m-toolbox/test/test_ao_noisegen1D.m @ 0:f0afece42f48

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/noisegen1D
2 %
3 % L. Ferraioli 10-11-08
4 %
5 % $Id: test_ao_noisegen1D.m,v 1.7 2010/05/04 13:30:35 luigi Exp $
6 %
7
8 %% Make white noise and compare results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9
10 % Description:
11 % 1) Generates a random series of data (white)
12 % 2) Build a pzm
13 % 3) Calculate equivalent one sided data psd from pzm response
14 % 4) Generates a miir object from pzm
15 % 5) Filer white data with miir to generate reference colored noise
16 % 6) Generates colored noise with noisegen1D
17 % 7) calculated psd of generated and reference data
18 % 8) check result by plotting
19
20 % 1)
21 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 10000));
22
23 % PSD model
24 % 2)
25 pzm = pzmodel(1, {0.01}, {0.1});
26 % 3)
27 mod = (abs(pzm.resp).^2).*(2/fs); % this corresponds to the theoretical one sided psd of the data
28 % 4)
29 fpzm = miir(pzm,plist('fs',fs));
30 % 5)
31 acr = filter(a,fpzm);
32
33 % 6) Noise generation
34 pl = plist(...
35 'model', mod, ...
36 'MaxIter', 10, ...
37 'PoleType', 3, ...
38 'MinOrder', 2, ...
39 'MaxOrder', 9, ...
40 'Weights', 2, ...
41 'Plot', false,...
42 'Disp', false,...
43 'MSEVARTOL', 1e-3,...
44 'FITTOL', 1e-5);
45
46 ac = noisegen1D(a, pl);
47
48 % 7)
49 acxx = ac.psd;
50 acrxx = acr.psd;
51 % 8)
52 iplot(acxx,acrxx,mod);
53
54 %% Noise generation from fsdata model object - output time series %%%%%%%%
55
56 % Description:
57 % 1) Generates a fsdata object to be used as psd model
58 % 2) Generates a random series of data (white)
59 % 3) Generates colored noise with noisegen1D
60 % 4) calculated psd of generated data
61 % 5) check result by plotting
62
63 % 1)
64 fs = 10; % sampling frequency
65 f = logspace(-6,log10(5),1000);
66 pl_mod1 = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
67 mod1 = ao(pl_mod1); % fsdata model object
68
69 % 2)
70 % generating white noise
71 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 10000));
72
73 % 3) Noise generation
74
75 pl1 = plist(...
76 'model', mod1, ...
77 'MaxIter', 30, ...
78 'PoleType', 2, ...
79 'MinOrder', 3, ...
80 'MaxOrder', 20, ...
81 'Weights', 3, ...
82 'Plot', false,...
83 'Disp', false,...
84 'MSEVARTOL', 1e-3,...
85 'FITTOL', 1e-5);
86
87 ac1 = noisegen1D(a1, pl1);
88
89 % 4)
90 acxx1 = ac1.psd;
91 % 5)
92 iplot(acxx1, mod1);
93
94 %% Noise generation from fsdata model object - output filter %%%%%%%%%%%%%
95
96 % Description:
97 % 1) Generates a fsdata object to be used as psd model
98 % 2) Generates coloring filter with noisegen1D
99 % 3) Generates white noise data
100 % 4) filter data
101 % 4) calculated psd of filtered data
102 % 5) check result by plotting
103
104 % 1) Generates a fsdata object to be used as psd model
105 fs = 10; % sampling frequency
106 f = logspace(-6,log10(5),1000);
107 pl_mod = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
108 mod = ao(pl_mod); % fsdata model object
109
110 % 2) noise coloring filter generation
111
112 pl1 = plist(...
113 'fs', fs, ...
114 'Iunits', '', ...
115 'Ounits', 'm', ....
116 'MaxIter', 30, ...
117 'PoleType', 2, ...
118 'MinOrder', 3, ...
119 'MaxOrder', 20, ...
120 'Weights', 3, ...
121 'Plot', false,...
122 'Disp', false,...
123 'MSEVARTOL', 1e-3,...
124 'FITTOL', 1e-5);
125
126 fil = noisegen1D(mod, pl1);
127
128 % 3) generating white noise
129 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5));
130
131 % 4) filter white noise
132 b = filter(a,fil);
133
134 % 4) psd estimation
135 bxx = b.psd(plist('navs',8,'olap',50,'order',1));
136 % 5)
137 iplot(bxx, mod);
138
139 %% Noise generation from fsdata model object - output filter 2 %%%%%%%%%%%%%
140
141 % Description:
142 % 1) Generates a fsdata object to be used as psd model
143 % 2) Generates coloring filter with noisegen1D
144 % 3) Generates white noise data
145 % 4) filter data
146 % 4) calculated psd of filtered data
147 % 5) check result by plotting
148
149 % 1) Generates a fsdata object to be used as psd model
150 fs = 10; % sampling frequency
151 f = logspace(-6,log10(5),1000);
152 pl_mod = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
153 mod = ao(pl_mod); % fsdata model object
154
155 % 2) noise coloring filter generation
156
157 pl1 = plist(...
158 'fs', fs, ...
159 'Iunits', '', ...
160 'Ounits', 'm', ....
161 'MaxIter', 30, ...
162 'PoleType', 2, ...
163 'MinOrder', 3, ...
164 'MaxOrder', 20, ...
165 'Weights', 3, ...
166 'Plot', false,...
167 'Disp', false,...
168 'MSEVARTOL', 1e-3,...
169 'FITTOL', 1e-5);
170
171 fil = noisegen1D(mod, mod, pl1);
172
173 % 3) generating white noise
174 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5));
175
176 % 4) filter white noise
177 b = filter(a,fil(1));
178
179 % 4) psd estimation
180 bxx = b.psd(plist('navs',8,'olap',50,'order',1));
181 % 5)
182 iplot(bxx, mod);
183
184 %% Noise generation from fsdata model object - output filter 3 %%%%%%%%%%%%%
185
186 % Description:
187 % 1) Generates a fsdata object to be used as psd model
188 % 2) Generates coloring filter with noisegen1D
189 % 3) Generates white noise data
190 % 4) filter data
191 % 4) calculated psd of filtered data
192 % 5) check result by plotting
193
194 % 1) Generates a fsdata object to be used as psd model
195 fs = 10; % sampling frequency
196 f = logspace(-6,log10(5),1000);
197 pl_mod = plist('fsfcn', '0.01./(0.01+f)', 'f', f);
198 mod = ao(pl_mod); % fsdata model object
199
200 % 2) noise coloring filter generation
201
202 pl1 = plist(...
203 'fs', fs, ...
204 'Iunits', '', ...
205 'Ounits', 'm', ....
206 'MaxIter', 30, ...
207 'PoleType', 2, ...
208 'MinOrder', 3, ...
209 'MaxOrder', 20, ...
210 'Weights', 3, ...
211 'Plot', false,...
212 'Disp', false,...
213 'MSEVARTOL', 1e-3,...
214 'FITTOL', 1e-5);
215
216 [fil1, fil2] = noisegen1D(mod, mod, pl1);
217
218 % 3) generating white noise
219 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5));
220
221 % 4) filter white noise
222 b = filter(a,fil1);
223
224 % 4) psd estimation
225 bxx = b.psd(plist('navs',8,'olap',50,'order',1));
226 % 5)
227 iplot(bxx, mod);
228
229 % END