Mercurial > hg > ltpda

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