Mercurial > hg > ltpda
comparison m-toolbox/test/test_ao_firwhiten.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_ao_firwhiten | |
| 2 % A test script for ao/firwhiten | |
| 3 % | |
| 4 % M Hewitson 21-04-08 | |
| 5 % | |
| 6 % $Id: test_ao_firwhiten.m,v 1.5 2009/08/26 18:27:06 luigi Exp $ | |
| 7 % | |
| 8 | |
| 9 mc | |
| 10 | |
| 11 %% Make an AO | |
| 12 | |
| 13 fs = 10; | |
| 14 nsecs = 1e3; | |
| 15 | |
| 16 % pole zero model | |
| 17 pzm = pzmodel(2, [pz(0.1, 2) pz(0.5, 1)] , [pz(0.8) pz(1, 2)]); | |
| 18 | |
| 19 % parameter list for ltpda_noisegen | |
| 20 pl = plist(); | |
| 21 pl = append(pl, param('nsecs', nsecs)); | |
| 22 pl = append(pl, param('fs', fs)); | |
| 23 | |
| 24 % calling the noisegenerator | |
| 25 a = fngen(pzm, pl); | |
| 26 % a = set(a, 'yunits', 'V'); | |
| 27 | |
| 28 % Make series of sine waves | |
| 29 | |
| 30 % s = ao(plist('waveform', 'sine wave', ... | |
| 31 % 'A', [0.1 1 0.1], ... | |
| 32 % 'f', [3.3 0.42 1.23], ... | |
| 33 % 'phi', [0 0 0], ... | |
| 34 % 'fs', fs, 'nsecs', nsecs)); | |
| 35 | |
| 36 s = ao(plist('tsfcn', '0.05*sin(2*pi*1.23*t)', ... | |
| 37 'fs', fs, 'nsecs', nsecs)); | |
| 38 % s = ao(plist('waveform', 'sine wave', ... | |
| 39 % 'A', [0.05], ... | |
| 40 % 'f', [1.23], ... | |
| 41 % 'phi', [0], ... | |
| 42 % 'fs', fs, 'nsecs', nsecs)); | |
| 43 | |
| 44 % Add signal and noise | |
| 45 a = a + s; | |
| 46 | |
| 47 %% Whiten | |
| 48 Nfft = fs*100; | |
| 49 [aw, ffs, nfs] = firwhiten(a, plist('Ntaps', 1024, 'Nfft', Nfft)); | |
| 50 | |
| 51 % plot filter response | |
| 52 iplot(resp(ffs)) | |
| 53 | |
| 54 | |
| 55 %% Remove end effects | |
| 56 | |
| 57 spl = plist('split_type', 'times', 'times', [10 -20]); | |
| 58 a = timeshift(split(a, spl)); | |
| 59 aw = timeshift(split(aw, spl)); | |
| 60 | |
| 61 %% Compare time-series | |
| 62 iplot(a,aw, plist('Arrangement', 'subplots', 'XRanges', {'All', [0 10]})) | |
| 63 | |
| 64 %% Compare spectra | |
| 65 | |
| 66 pxx = pwelch(a, aw, plist('Nfft', fs*100)); | |
| 67 iplot(pxx) | |
| 68 | |
| 69 %% Help test | |
| 70 | |
| 71 fs = 1; | |
| 72 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5, 'yunits','m')); | |
| 73 | |
| 74 % filter | |
| 75 pl = plist('type','bandpass', ... | |
| 76 'order',3, ... | |
| 77 'gain',1, ... | |
| 78 'fc',[0.03 0.1], ... | |
| 79 'fs',fs); | |
| 80 ft = miir(pl); | |
| 81 | |
| 82 % coloring white noise | |
| 83 af = filter(a, ft); | |
| 84 | |
| 85 % firwhiten | |
| 86 [aw, ffs, nfs] = firwhiten(af, plist('Ntaps', 5000, 'Nfft', 1e5, 'BW', 5)); | |
| 87 | |
| 88 awxx = aw.psd; | |
| 89 afxx = af.psd; | |
| 90 | |
| 91 iplot(afxx,awxx) | |
| 92 | |
| 93 % END |