Mercurial > hg > ltpda

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