Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_gap_problems.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 show_gap_problem
2 close all
3 mc
4 %% common params
5 ns = 1e4;
6 fs = 100;
7 ts = 1/fs;
8 tf = ns-ts;
9 Ns = ns*fs;
10 nfft = Ns;
11
12 %% no gaps
13 % p1 = pzmodel(1e-8,[pz(1e-5) pz(1e-5) pz(1e-5) pz(2e-5) pz(2e-5) pz(1) pz(1) pz(1)],...
14 % [ pz(2e-4) pz(2e-4) pz(1e-4) pz(1e-4) pz(1e-4) pz(1e-3) pz(1e-2)]);
15 % p1 = pzmodel(1e-8,[pz(1e-5) pz(2e-5) pz(3e-5) pz(4e-5) pz(5e-5) pz(0.9) pz(1) pz(1.1) pz(1.2) pz(1.3)],...
16 % [ pz(2e-4) pz(3e-4) pz(4e-4) pz(5e-4) pz(1e-4) pz(1e-3) pz(1e-2)]);
17 p1 = pzmodel(0.1,[pz(1e-5) pz(1e-5) pz(1e-5) pz(2e-5) pz(3e-5) pz(4e-5) pz(5e-5) pz(0.9) pz(1) pz(1.1) pz(1.2) pz(1.3)],...
18 [ pz(2e-4) pz(3e-4) pz(4e-4) pz(5e-4) pz(1e-4) pz(1e-3) pz(1e-2)]);
19 % iplot(resp(p1))
20 p2 = pzmodel(1e-9,[ pz(1e-5) pz(2e-5) pz(3e-5) pz(4e-5) pz(5e-5) pz(0.9) pz(1) pz(1.1) pz(1.2) pz(1.3)],...
21 [ pz(2e-4) pz(3e-4) pz(4e-4) pz(5e-4) pz(1e-4) pz(1e-3) pz(1e-2)]);
22
23 %% data generation
24 s1 = ssm(p1);
25 s1.modiftimestep(ts);
26 s1.isstable
27 a00hf = s1.simulate(plist('Nsamples', Ns,'noise variable names',{'U>1 '} ,'covariance' , 1, 'return outputs', {'Y>1 '} ));
28 a02hf = ao(p1,ns,fs );
29
30 %% plotting data generation results
31 p_hf = plist('Scale', 'ASD', 'Order',1,'Nfft', nfft);
32 ps_hf = plist('Scale', 'ASD', 'Order',1,'Nfft', floor(nfft/3));
33 % iplot([psd(a00hf,p_hf) resp(p1) psd(a02hf, p_hf)])
34
35 %% downsampling
36 dns = true;
37 k = 500;
38 if dns
39 opt = plist('fsout', fs/k);
40 a00 = a00hf.dsmean(opt);
41 a02 = a02hf.dsmean(opt);
42 ns = k*floor(ns/k);
43 fs = fs/k;
44 ts = ts*k;
45 tf = ns-ts;
46 Ns = floor(Ns/k);
47 nfft = Ns;
48 end
49
50 %% acceleration and spectrum
51 a00hf = a00hf.diff;
52 a00hf = a00hf.diff;
53 a02hf = a02hf.diff;
54 a02hf = a02hf.diff;
55 a00 = a00.diff;
56 a00 = a00.diff;
57 a02 = a02.diff;
58 a02 = a02.diff;
59 p = plist('Scale', 'ASD', 'Order',1,'Nfft', nfft);
60 ps = plist('Scale', 'ASD', 'Order',1,'Nfft', floor(nfft/3));
61 iplot([resp(p2) psd(a00hf,ps_hf) psd(a02hf, ps_hf) psd(a00,ps) psd(a02, ps)])
62
63 %% baseline : periodic gaps 1s / 249s / 400 cycles
64 kick1 = double(~(mod(0:ts:tf, 250)>min(245,250-ts-1) )).';
65 a1k = ao(tsdata(kick1,fs ));
66 a10 = a00.*a1k;
67
68 %% pb2 : not so periodic gaps : <dn> = 5s, <Tauk+ taud> = 1/500*sqrt(10/fs)
69 time2 = cumsum(ts*ones(1,Ns) + ts*0.1*randn(1,Ns));
70 kick2 = double(~(mod(time2, 250)>min(245,250-ts-1) )).';
71 a2k = ao(tsdata(kick2,fs ));
72 a20 = a00.*a2k;
73
74 %% pb3 : time series too short : 100 cycles
75 a30 = ao(tsdata(a10.data.y(1:floor(Ns/4)), fs));
76 p3 = plist('Scale', 'ASD', 'Order',1,'Nfft', floor(nfft/4));
77 ps3 = plist('Scale', 'ASD', 'Order',1,'Nfft', floor(nfft/20));
78
79 %% pb4 : longer gaps : 10s / 230s
80 kick1 = double(~(mod(0:ts:tf, 250)> min(200,250-ts-1) )).';
81 a4k = ao(tsdata(kick1,fs ));
82 a40 = a00.*a4k;
83 % iplot([psd(a00,p) psd(a10,p) psd(a20,p) psd(a30,p3) psd(a40,p)])
84 iplot([psd(a00,ps) psd(a10,ps) psd(a20,ps) psd(a30,ps3) psd(a40,ps)])
85
86 %% gap-fill all 4 examples
87 tic
88 a10g = a10.gapfillingoptim(p2);
89 t(1)=toc;
90 tic
91 a20g = a20.gapfillingoptim(p2);
92 t(2)=toc;
93 tic
94 a30g = a30.gapfillingoptim(p2);
95 t(3)=toc;
96 tic
97 try
98 a40g = a40.gapfillingoptim(p2);
99 t(4)=toc;
100 a40gworked = true;
101 catch
102 a40gworked = false;
103 end
104
105 %% plotting again
106 display('time')
107 display(t)
108
109 if a40gworked
110 pplot1= plist('LineColors',{'r' 'r' 'r' 'b' 'k' 'g ' 'r' 'b' 'k' 'g' } ,'LineStyles', {'-' '-' '-' '-' '-' '-' ':' ':' ':' ':'},'LineWidths', {2 2 1 1 1 1 1 1 1 1});
111 iplot([abs(resp(p2)) psd(a00,ps) psd(a10,ps) psd(a20,ps) psd(a30,ps3) psd(a40,ps) psd(a10g,ps) psd(a20g,ps) psd(a30g,ps3) psd(a40g,ps)], pplot1)
112 else
113 pplot1= plist('LineColors',{'r' 'r' 'r' 'b' 'k' 'g' 'r' 'b' 'k' } ,'LineStyles', {'-' '-' '-' '-' '-' '-' ':' ':' ':'},'LineWidths', {2 2 1 1 1 1 1 1 1});
114 iplot([abs(resp(p2)) psd(a00,ps) psd(a10,ps) psd(a20,ps) psd(a30,ps3) psd(a40,ps) psd(a10g,ps) psd(a20g,ps) psd(a30g,ps3) ], pplot1)
115 end
116
117
118