Mercurial > hg > ltpda

comparison m-toolbox/test/diagnostics/test_ltpda_arma_time.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 function test_ltpda_arma_time()
2 % A test script for the ltpda_arma_time function
3 %
4 % miquel 20-02-08
5 %
6 % $Id: test_ltpda_arma_time.m,v 1.1 2008/03/04 12:55:05 miquel Exp $
7 %
8
9 %% Make test AOs
10
11 fs=1;
12 nsecs=10000;
13
14 pl = plist();
15 pl = append(pl, param('nsecs',nsecs));
16 pl = append(pl, param('fs',fs));
17 pl = append(pl,param('tsfcn','(heaviside(t-1000.5)).*(1-(exp(-(t-1000.5)/300))+(exp(-(t-1000.5)/500)-1))'));
18
19 ain = ao(pl);
20
21 % Filter input time-series
22 pl = append(pl, param('type', 'highpass'));
23 pl = append(pl, param('ripple', 5e-1));
24 pl = append(pl, param('order', 1));
25 pl = append(pl, param('gain', 1e-1));
26 pl = append(pl, param('fc', 1e-3));
27
28 f = miir(pl);
29 [aout, fout] = filter(ain,plist(param('filter',f)));
30
31 % Add independent noise to both channels
32 pn1 = plist('waveform', 'noise','type','Uniform','fs', fs, 'nsecs', nsecs);
33 an1 = ao(pn1);
34
35 ain = ain + 1e-3*an1;
36
37 pn2 = plist('waveform', 'noise','type','Uniform','fs', fs, 'nsecs', nsecs);
38 an2 = ao(pn2);
39
40 aout = aout + 1e-3*an2;
41
42
43 %% Find ARMA parameters for the transfer function
44
45 % Parameter list for ltpda_arma_time
46 pl = append(pl, param('MaxIter',1e3));
47 pl = append(pl, param('MaxFunEvals',1e3));
48 pl = append(pl, param('TolX',1e-7));
49 pl = append(pl, param('TolFun',1e-7));
50 pl = append(pl, param('ARMANum',2));
51 pl = append(pl, param('ARMADen',1));
52
53 % Use ltpda_arma_time
54 ao_arma = ltpda_arma_time(ain, aout, pl);
55 %
56 disp(sprintf('\r! Filter parameters:\r'))
57 disp(sprintf('%d\t', fout.a,fout.b(2:length(fout.b))))
58 disp(sprintf('\r! Estimated parameters:\r'))
59 disp(sprintf('%d\t',ao_arma.data.y(1:length(ao_arma.data.y)-1)))
60
61
62 % Compute fitted output
63 p = find(pl, 'ARMANum');
64 q = find(pl, 'ARMADen');
65
66 f = miir([ao_arma.data.y(1:p)],[1 ao_arma.data.y(p+1:p+q)],find(pl,'fs'));
67 [afit, ffit] = filter(ain,plist(param('filter',f)));
68
69 % Plot
70 iplot(ain,aout,afit,plist('Legends', {'Input','Output','Fit'}))
71
72 % %% Plot history
73 % figure
74 % plot(ao_arma.hist)
75 %
76 % %% Reproduce from history
77 %
78 % % Write an m-file from AO
79 % ao2m(ao_arma, 'ltpda_test.m');
80 %
81 % % now run it
82 % clear all;
83 % a_out = ltpda_test;
84 %
85 % % iplot(a_out)
86 %
87 % figure
88 % plot(a_out.hist)
89
90 %% Delete test files
91 % delete('ltpda_test.m');
92