Mercurial > hg > ltpda
comparison m-toolbox/test/test_mfir_from_spec.m @ 0:f0afece42f48
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_mfir_from_spec() | |
2 | |
3 % TEST_MFIR_FROM_SPEC test script for mfir constructor from AO/fsdata. | |
4 % | |
5 % M Hewitson 14-05-07 | |
6 % | |
7 % $Id: test_mfir_from_spec.m,v 1.4 2008/09/19 14:03:50 ingo Exp $ | |
8 % | |
9 | |
10 Ntaps = 512; | |
11 | |
12 %% Make ao | |
13 | |
14 % pole zero model | |
15 pzm = pzmodel(1, [pz(0.1, 2) pz(0.5, 40)] , [pz(.03) pz(1.55, 25)]); | |
16 | |
17 % Compute time-series | |
18 Nsecs = 1e4; % Number of seconds desired | |
19 swin = specwin('Kaiser', 10, 150); % Blending window | |
20 x12 = fngen(pzm, plist('Win', swin, 'Nsecs', Nsecs)); % create time-series | |
21 | |
22 | |
23 %% Make spectrum | |
24 | |
25 pl = plist('Nfft', round(200*x12.data.fs)); | |
26 x12xx = pwelch(x12, pl); | |
27 | |
28 %% Make filter | |
29 | |
30 ffs = mfir(x12xx, plist('N', Ntaps, 'Win', specwin('Hamming', 10))); | |
31 | |
32 rfs = resp(ffs, plist('f', [x12xx.data.x].')); | |
33 | |
34 pl = plist('nsecs', Nsecs,'fs', x12.data.fs,'tsfcn', 'randn(size(t))'); | |
35 anoise = ao(pl); | |
36 | |
37 tic | |
38 anew = filter(4*anoise, ffs); | |
39 toc | |
40 pl = plist('Nfft', round(200*x12.data.fs)); | |
41 anewxx = pwelch(anew, pl); | |
42 | |
43 %% | |
44 iplot(x12xx, rfs, anewxx) | |
45 | |
46 iplot(x12xx./anewxx) | |
47 utils.plottools.yscale('lin') | |
48 | |
49 %% Make inverse filter | |
50 ffs = mfir(1./x12xx, plist('N', Ntaps, 'Win', specwin('Hamming', 10))); | |
51 | |
52 rfs = resp(ffs, plist('f', x12xx.data.x)); | |
53 iplot(rfs) | |
54 | |
55 x12wht = filter(x12, ffs); | |
56 pl = plist('Nfft', round(200*x12.data.fs)); | |
57 x12whtxx = pwelch(x12wht, pl); | |
58 | |
59 iplot(x12whtxx) | |
60 | |
61 return | |
62 | |
63 %% Write an m-file from AO | |
64 pth = get_curr_m_file_path(mfilename); | |
65 ao2m(x12l, [pth, 'test.m']); | |
66 | |
67 % now run it | |
68 clear all | |
69 test | |
70 | |
71 figure | |
72 | |
73 plot(a1) | |
74 figure | |
75 plot(a1.hist) | |
76 | |
77 % END | |
78 | |
79 | |
80 % END |