--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/MDC2/test_mdc2_filtimpresp.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,215 @@
+% A test script for to check the filters impulse response
+% 
+% 
+% L. Ferraioli 04-02-2009
+% 
+% $Id: test_mdc2_filtimpresp.m,v 1.3 2009/04/20 14:34:12 luigi Exp $
+% 
+
+%% General use variables and vectors
+
+fs = 1;
+f = logspace(-7,log10(fs/2),500);
+f = f.';
+Nsecs = 1e6; % number of seconds
+
+%% MDC2 Models
+
+b = ao(plist('built-in','mdc2r2_fd_ltpnoise','f',f));
+
+%%
+
+iplot(b)
+
+%%
+
+CSD = [b(1) b(2);conj(b(2)) b(3)];
+MSC = b(2).*conj(b(2))./(b(1).*b(3));
+
+%% some ploting
+
+iplot(CSD(1,1),CSD(2,2),MSC)
+iplot(CSD(1,1),CSD(2,2),CSD(1,2),CSD(2,1))
+
+%% Make white noise
+
+a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
+a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
+
+%% Noise generation
+
+plng = plist(...
+    'csd11', CSD(1,1), ...
+    'csd12', CSD(1,2), ...
+    'csd21', CSD(2,1), ...
+    'csd22', CSD(2,2), ...
+    'MaxIter', 90, ...
+    'PoleType', 2, ...
+    'MinOrder', 35, ...
+    'MaxOrder', 50, ...
+    'Weights', 2, ...
+    'Plot', false,...
+    'MSEVARTOL', 1e-2,...
+    'FITTOL', 1e-5,...
+    'UseSym', 0,...
+    'Disp', false);
+
+[ac1,ac2] = noisegen2D(a1,a2, plng);
+
+% extracting filters
+h11 = ac1.procinfo.find('Filt11');
+h12 = ac1.procinfo.find('Filt12');
+h21 = ac2.procinfo.find('Filt21');
+h22 = ac2.procinfo.find('Filt22');
+
+%% Calculate Impulse response with LTPDA filter
+
+% LTPDA filter initialize the filters before the filtering process
+
+imp = zeros(Nsecs,1);
+imp(1,1) = 1;
+imp = tsdata(imp,fs);
+imp = ao(imp);
+imp.setYunits('m');
+imp.setXunits('s');
+
+irh11 = filter(imp,h11);
+irh12 = filter(imp,h12);
+irh21 = filter(imp,h21);
+irh22 = filter(imp,h22);
+
+irh11.setName;
+irh12.setName;
+irh21.setName;
+irh22.setName;
+
+% iplot(irh11,irh12,irh21,irh22)
+%% Plot impulse response
+
+iplot(irh11,plist('Yscales','linear'))
+iplot(irh12,plist('Yscales','linear'))
+iplot(irh21,plist('Yscales','linear'))
+iplot(irh22,plist('Yscales','linear'))
+
+%% Impluse response with Matlab filter
+
+% No filter initialization
+
+idat = [1; zeros(999999,1)];
+
+Nfilt11 = numel(h11);
+num11 = h11(1).a;
+den11 = h11(1).b;
+irh11 = filter(num11,den11,idat);
+
+for ii = 2:Nfilt11
+  num11 = h11(ii).a;
+  den11 = h11(ii).b;
+  todat = filter(num11,den11,idat);
+  irh11 = irh11 + todat;
+end
+
+Nfilt12 = numel(h12);
+num12 = h12(1).a;
+den12 = h12(1).b;
+irh12 = filter(num12,den12,idat);
+
+for ii = 2:Nfilt12
+  num12 = h12(ii).a;
+  den12 = h12(ii).b;
+  todat = filter(num12,den12,idat);
+  irh12 = irh12 + todat;
+end
+
+Nfilt21 = numel(h21);
+num21 = h21(1).a;
+den21 = h21(1).b;
+irh21 = filter(num21,den21,idat);
+
+for ii = 2:Nfilt21
+  num21 = h21(ii).a;
+  den21 = h21(ii).b;
+  todat = filter(num21,den21,idat);
+  irh21 = irh21 + todat;
+end
+
+Nfilt22 = numel(h22);
+num22 = h22(1).a;
+den22 = h22(1).b;
+irh22 = filter(num22,den22,idat);
+
+for ii = 2:Nfilt22
+  num22 = h22(ii).a;
+  den22 = h22(ii).b;
+  todat = filter(num22,den22,idat);
+  irh22 = irh22 + todat;
+end
+  
+%%
+
+figure
+plot(irh11)
+
+figure
+plot(irh12)
+
+figure
+plot(irh21)
+
+figure
+plot(irh22)
+
+%% Filter state
+
+ac11 = filter(a1,h11);
+h11 = ac11.procinfo.find('Filters'); % prooagate histout
+ac12 = filter(a2,h12);
+h12 = ac12.procinfo.find('Filters');
+ac21 = filter(a1,h21);
+h21 = ac21.procinfo.find('Filters');
+ac22 = filter(a2,h22);
+h22 = ac22.procinfo.find('Filters');
+
+iplot(ac11)
+iplot(ac12)
+iplot(ac21)
+iplot(ac22)
+
+%% Filter state evolution
+
+Nrun = 20;
+
+for jj = 1:Nrun % filter state evolution
+  disp(num2str(jj))
+  ac11 = filter(a1,h11);
+  h11 = ac11.procinfo.find('Filters'); % prooagate histout
+  ac12 = filter(a2,h12);
+  h12 = ac12.procinfo.find('Filters');
+  ac21 = filter(a1,h21);
+  h21 = ac21.procinfo.find('Filters');
+  ac22 = filter(a2,h22);
+  h22 = ac22.procinfo.find('Filters');
+end
+%%
+iplot(ac11)
+iplot(ac12)
+iplot(ac21)
+iplot(ac22)
+
+%%
+
+ac1 = ac11 + ac12;
+ac2 = ac21 + ac22;
+
+iplot(ac1)
+iplot(ac2)
+
+%%
+
+ac1xx = ac1.lpsd(plist('Nfft',1e5,'order',2));
+ac2xx = ac2.lpsd(plist('Nfft',1e5,'order',2));
+
+%%
+iplot(CSD(1,1),ac1xx,CSD(2,2),ac2xx)
+
+