Mercurial > hg > ltpda
view m-toolbox/test/MDC2/test_mdc2_filtimpresp.m @ 30:317b5f447f3e database-connection-manager
Update workspaceBrowser
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% 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)