+
−% 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)
+
−
+
−