% Test script for ao/filter with matrix input
% 
% L. Ferraioli 02-12-08
% 
% $Id: test_ao_filtfilt.m,v 1.3 2011/05/12 13:23:35 luigi Exp $

%% Test a filter matrix with a single filter object

% get filter
pl = plist('type', 'lowpass', ...
                           'order', 2,         ...
                           'gain',  1.0,       ...
                           'fs',    10,        ...
                           'fc',    0.2);

filt = miir(pl);
fresp = resp(filt,plist('f',[logspace(-4,0,100)].'));

% put filter into a matrix
mfilt = matrix(filt);

% gen white noise data
b = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e4));

% do filtering with matrix
c = filter(b,mfilt);
sc = tfe(b,c);
sc.setName;

% do filtering with matrix inside plist
d = filtfilt(b,plist('filter',mfilt));
sd = tfe(b,d);
sd.setName;

% test rebuild
g = rebuild(c);
sg = tfe(b,g);
sg.setName;

h = rebuild(d);
sh = tfe(b,h);
sh.setName;

% check
iplot(sc,sd,sg,sh,fresp)

%% Test a filter with a single filterbank object

% get filter and response
pl = plist('type', 'lowpass', ...
                           'order', 2,         ...
                           'gain',  1.0,       ...
                           'fs',    10,        ...
                           'fc',    0.2);

filt = miir(pl);
fresp = resp(filt,plist('f',[logspace(-4,0,100)].'));

% generate a filterbank with zDomainFit
pl_fit = plist('FS',10,...
  'AutoSearch','on',...
  'StartPoles',[],...
  'StartPolesOpt','clin',...
  'maxiter',40,...
  'minorder',3,...
  'maxorder',25,...
  'weights',[],...
  'weightparam','abs',...
  'CONDTYPE','MSE',...
  'FITTOL',1e-4,... % check if MSE is lower than 1e-4
  'MSEVARTOL',1e-2,...
  'Plot','off',...
  'ForceStability','on',...
  'CheckProgress','off');

% Do fit
fbk_filt = zDomainFit(fresp, pl_fit);

% gen white noise data
b = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e4));

% do filtering with filterbank
c = filtfilt(b,fbk_filt);
sc = tfe(b,c);
sc.setName('direct filterbank')

d = filter(b,fbk_filt);
sd = tfe(b,d);
sd.setName('direct filter')

% % do filtering with filterbank inside plist
% d = filtfilt(b,plist('filter',fbk_filt));
% sd = tfe(b,d);
% sd.setName('plist filterbank')


% check
iplot(abs(sc),abs(sd).^2,abs(fresp).^2)


%% test filter with 2dim matrix

% test that filter trow an error when a N-dim matrix is input
pl1 = plist('type', 'lowpass', ...
                           'order', 2,         ...
                           'gain',  1.0,       ...
                           'fs',    10,        ...
                           'fc',    0.2);
                         
pl2 = plist('type', 'lowpass', ...
                           'order', 1,         ...
                           'gain',  1e-1,       ...
                           'fs',    10,        ...
                           'fc',    0.01);

filt1 = miir(pl1);
fresp1 = resp(filt1,plist('f',[logspace(-4,0,100)].'));
filt2 = miir(pl2);
fresp2 = resp(filt2,plist('f',[logspace(-4,0,100)].'));

mfilt = matrix(filt1,filt2);

% gen white noise data
b = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e4));

% do filtering with matrix
try
  c = filter(b,mfilt);
catch err
  disp(err.message)
end

try
  c = filter(b,plist('filter',mfilt));
catch err
  disp(err.message)
end

%% test filter with 2dim matrix

% test that filter trow an error when a N-dim matrix is input
pl1 = plist('type', 'lowpass', ...
                           'order', 2,         ...
                           'gain',  1.0,       ...
                           'fs',    10,        ...
                           'fc',    0.2);
                         
pl2 = plist('type', 'lowpass', ...
                           'order', 1,         ...
                           'gain',  1e-1,       ...
                           'fs',    10,        ...
                           'fc',    0.01);

filt1 = miir(pl1);
fresp1 = resp(filt1,plist('f',[logspace(-4,0,100)].'));
filt2 = miir(pl2);
fresp2 = resp(filt2,plist('f',[logspace(-4,0,100)].'));

% Fitting parameter list
pl_fit = plist('FS',10,...
  'AutoSearch','on',...
  'StartPoles',[],...
  'StartPolesOpt','clog',...
  'maxiter',40,...
  'minorder',3,...
  'maxorder',25,...
  'weights',[],...
  'weightparam','abs',...
  'CONDTYPE','MSE',...
  'FITTOL',1e-4,... % check if MSE is lower than 1e-4
  'MSEVARTOL',1e-2,...
  'Plot','off',...
  'ForceStability','on',...
  'CheckProgress','off');

% Do fit
mod = zDomainFit(fresp1,fresp2,pl_fit);

b = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e4));

try
  c = filter(b,mfilt);
catch err
  disp(err.message)
end

try
  c = filter(b,plist('filter',mfilt));
catch err
  disp(err.message)
end