0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % function test_ao_firwhiten
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 % A test script for ao/firwhiten
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % M Hewitson 21-04-08
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 % $Id: test_ao_firwhiten.m,v 1.5 2009/08/26 18:27:06 luigi Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 mc
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 %% Make an AO
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 fs = 10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 nsecs = 1e3;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 % pole zero model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 pzm = pzmodel(2, [pz(0.1, 2) pz(0.5, 1)] , [pz(0.8) pz(1, 2)]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 % parameter list for ltpda_noisegen
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 pl = plist();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 pl = append(pl, param('nsecs', nsecs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 pl = append(pl, param('fs', fs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 % calling the noisegenerator
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 a = fngen(pzm, pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 % a = set(a, 'yunits', 'V');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 % Make series of sine waves
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 % s = ao(plist('waveform', 'sine wave', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 % 'A', [0.1 1 0.1], ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 % 'f', [3.3 0.42 1.23], ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 % 'phi', [0 0 0], ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 % 'fs', fs, 'nsecs', nsecs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 s = ao(plist('tsfcn', '0.05*sin(2*pi*1.23*t)', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 'fs', fs, 'nsecs', nsecs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 % s = ao(plist('waveform', 'sine wave', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 % 'A', [0.05], ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 % 'f', [1.23], ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 % 'phi', [0], ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 % 'fs', fs, 'nsecs', nsecs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 % Add signal and noise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 a = a + s;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 %% Whiten
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 Nfft = fs*100;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 [aw, ffs, nfs] = firwhiten(a, plist('Ntaps', 1024, 'Nfft', Nfft));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 % plot filter response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 iplot(resp(ffs))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 %% Remove end effects
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 spl = plist('split_type', 'times', 'times', [10 -20]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 a = timeshift(split(a, spl));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 aw = timeshift(split(aw, spl));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 %% Compare time-series
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 iplot(a,aw, plist('Arrangement', 'subplots', 'XRanges', {'All', [0 10]}))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 %% Compare spectra
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 pxx = pwelch(a, aw, plist('Nfft', fs*100));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 iplot(pxx)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 %% Help test
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 fs = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5, 'yunits','m'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 % filter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 pl = plist('type','bandpass', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 'order',3, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 'gain',1, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 'fc',[0.03 0.1], ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 'fs',fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 ft = miir(pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 % coloring white noise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 af = filter(a, ft);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 % firwhiten
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 [aw, ffs, nfs] = firwhiten(af, plist('Ntaps', 5000, 'Nfft', 1e5, 'BW', 5));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 awxx = aw.psd;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 afxx = af.psd;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 iplot(afxx,awxx)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 % END |