Mercurial > hg > ltpda
diff 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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/m-toolbox/test/test_ao_firwhiten.m Wed Nov 23 19:22:13 2011 +0100 @@ -0,0 +1,93 @@ +% function test_ao_firwhiten +% A test script for ao/firwhiten +% +% M Hewitson 21-04-08 +% +% $Id: test_ao_firwhiten.m,v 1.5 2009/08/26 18:27:06 luigi Exp $ +% + +mc + +%% Make an AO + +fs = 10; +nsecs = 1e3; + +% pole zero model +pzm = pzmodel(2, [pz(0.1, 2) pz(0.5, 1)] , [pz(0.8) pz(1, 2)]); + +% parameter list for ltpda_noisegen +pl = plist(); +pl = append(pl, param('nsecs', nsecs)); +pl = append(pl, param('fs', fs)); + +% calling the noisegenerator +a = fngen(pzm, pl); +% a = set(a, 'yunits', 'V'); + +% Make series of sine waves + +% s = ao(plist('waveform', 'sine wave', ... +% 'A', [0.1 1 0.1], ... +% 'f', [3.3 0.42 1.23], ... +% 'phi', [0 0 0], ... +% 'fs', fs, 'nsecs', nsecs)); + +s = ao(plist('tsfcn', '0.05*sin(2*pi*1.23*t)', ... + 'fs', fs, 'nsecs', nsecs)); +% s = ao(plist('waveform', 'sine wave', ... +% 'A', [0.05], ... +% 'f', [1.23], ... +% 'phi', [0], ... +% 'fs', fs, 'nsecs', nsecs)); + +% Add signal and noise +a = a + s; + +%% Whiten +Nfft = fs*100; +[aw, ffs, nfs] = firwhiten(a, plist('Ntaps', 1024, 'Nfft', Nfft)); + +% plot filter response +iplot(resp(ffs)) + + +%% Remove end effects + +spl = plist('split_type', 'times', 'times', [10 -20]); +a = timeshift(split(a, spl)); +aw = timeshift(split(aw, spl)); + +%% Compare time-series +iplot(a,aw, plist('Arrangement', 'subplots', 'XRanges', {'All', [0 10]})) + +%% Compare spectra + +pxx = pwelch(a, aw, plist('Nfft', fs*100)); +iplot(pxx) + +%% Help test + +fs = 1; +a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5, 'yunits','m')); + +% filter +pl = plist('type','bandpass', ... + 'order',3, ... + 'gain',1, ... + 'fc',[0.03 0.1], ... + 'fs',fs); +ft = miir(pl); + +% coloring white noise +af = filter(a, ft); + +% firwhiten +[aw, ffs, nfs] = firwhiten(af, plist('Ntaps', 5000, 'Nfft', 1e5, 'BW', 5)); + +awxx = aw.psd; +afxx = af.psd; + +iplot(afxx,awxx) + +% END \ No newline at end of file