--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_noisegen_vpa.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,61 @@
+function test_noisegen_vpa()
+
+% TEST_NOISEGEN tests the Franklin noise-generator algorithm.
+% 
+% this is an example from 'Generation of Random time series with prescribed
+% spectra' by Gerhard Heinzel
+% 
+% $Id: test_noisegen_vpa.m,v 1.2 2008/09/19 14:03:50 ingo Exp $% 
+
+
+
+%% poles and zeros
+p1 = pz(1);
+p2 = pz(10,2);
+p3 = pz(1,4);
+
+z1 = pz(5);
+z2 = pz(50);
+% z2 = pz(3);
+gain =1; 
+%% convert to get a and b coefficients
+%pole zero model
+pzm = pzmodel(gain, [p1 p2 p3],[z1 z2]);
+% parameter list for ltpda_noisegen
+pl = plist();
+pl = append(pl, param('nsecs', 7200));
+pl = append(pl, param('fs', 100));
+pl = append(pl, param('pzmodel', pzm));
+pl = append(pl, param('ndigits', 30));
+%% calling the noisegenerator
+% [b, pl1, pl2]  = ltpda_noisegen(pl);
+tic
+[b, pl1, pl2]  = ao.ngsetup_vpa(pl);
+pl
+b = set(b, 't0', '2007-07-19 12:00:00');
+toc
+tic
+[b2, pl3, pl4] = ltpda_noisegen_vpa(pl1,pl2);
+b2 = set(b2, 't0', '2007-07-19 14:00:00');
+toc
+%% Make LPSD  
+% Window function
+w = specwin('Kaiser', 10, 150);
+% parameter list for lpsd
+pl = plist();
+pl = append(pl, param('Kdes', 100));
+pl = append(pl, param('Kmin', 2));
+pl = append(pl, param('Jdes', 200));
+pl = append(pl, param('Win', w));
+pl = append(pl, param('scale','ASD'));
+% LPSD
+alpsd  = ltpda_lpsd(b, pl); % matlab
+a_vpa  = ltpda_lpsd(b2, pl); % matlab
+
+% a2lpsd = ltpda_lpsd(b2, pl);
+% % 
+%% Plot
+iplot([b b2])
+iplot([alpsd a_vpa abs(resp(pzm))])
+
+