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