0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % Test file for ltpda_timedomain fit. This does a noise subtraction from
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 % engineering model data.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % M Hewitson 18-03-07
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 % $Id: test_timedomainfit.m,v 1.2 2008/07/29 15:11:44 ingo 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 %% Load data into analysis objects
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 disp('*** Loading data....')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 x12 = ao('~/frommario/Data/x12.txt'); %x12
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 phi1 = ao('~/frommario/Data/phi1.txt'); %phi1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 eta1 = ao('~/frommario/Data/eta1.txt'); %eta1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 phi12 = ao('~/frommario/Data/phi12.txt'); %phi12
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 eta12 = ao('~/frommario/Data/eta12.txt'); %eta12
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 disp('*** Done.')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 % get the sample frequency for designing bandpass filter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 d1 = x12.data;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 fs = d1.fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 %% Load bandpass
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 pl = plist();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 pl = append(pl, param('type', 'bandpass'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 % pl = append(pl, param('fs', fs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 pl = append(pl, param('fc', [0.003 0.03]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 pl = append(pl, param('order', 2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 bp = miir(pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 %% Filter data and remove startup transients
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 % parameter list for split()
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 % to remove startup transient of filters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 spl = plist(param('times', [1000 2000]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 % parameter list for filter()
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 fpl = plist(param('filter', bp));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 % filter and split data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 x12f = split(filter(x12, fpl), spl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 phi1f = split(filter(phi1, fpl), spl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 eta1f = split(filter(eta1, fpl), spl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 phi12f = split(filter(phi12, fpl), spl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 eta12f = split(filter(eta12, fpl), spl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 %% Make Timedomainfit
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 coeffsAO = timedomainfit(x12f, phi1f, eta1f, phi12f, eta12f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 %% Make linear combination
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 x12_fit = lincom(phi1, eta1, phi12, eta12, coeffsAO);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 x12ns = x12 - x12_fit;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 %% Make LPSD of each
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 % Window function
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 w = specwin('Kaiser', 10, 150);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 % parameter list for lpsd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 pl = plist();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 pl = append(pl, param('Kdes', 100));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 pl = append(pl, param('Kmin', 2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 pl = append(pl, param('Jdes', 500));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 pl = append(pl, param('Win', w));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 pl = append(pl, param('order', 1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 % LPSD of filtered x12
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 x12fxx = lpsd(x12f, pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 x12xx = lpsd(x12, pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 x12nsxx = lpsd(x12ns, pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 %% Plot results
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 iplot([x12fxx x12xx x12nsxx])
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 legend('Bandpassed X12', 'X12', 'X12 with noise subtracted');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 plot(x12ns.hist)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 save(x12fxx, 'x12fxx.xml');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 save(x12xx, 'x12xx.xml');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 save(x12nsxx, 'x12nsxx.xml');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 return
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 %% Reproduce from history
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 % Write an m-file from AO
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 ao2m(x12nsxx, 'test.m');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 %% now run it
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 test
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 plot(a1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 figure
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 plot(a1.hist)
|