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))])