--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/mtxiirresp.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,74 @@
+% MTXIIRRESP calculate iir resp by matrix product
+% 
+% INPUT:
+%       fil: a miir filter or a vector of filoters
+%       freq: the frequency vector
+%       fs: sampling frequency
+%       bank: filterbank type parameter ('parallel', 'serial' or 'none'). 
+%       If it is left empty, it provide a vector of responses at the
+%       output.
+%
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function rsp = mtxiirresp(fil,freq,fs,bank)
+  
+  T = 1/fs;
+  
+  Nfil = numel(fil);
+  
+  if Nfil>1
+    % check for different orders
+    sza = zeros(Nfil,1);
+    szb = zeros(Nfil,1);
+    for jj=1:Nfil
+      sza(jj) = numel(fil(jj).a);
+      szb(jj) = numel(fil(jj).b);
+    end
+    Ncl = max([sza;szb]);
+  else
+    sza = numel(fil.a);
+    szb = numel(fil.b);
+    Ncl = max([sza szb]);
+  end
+  
+  % init A and B
+  A = zeros(Nfil,Ncl);
+  B = zeros(Nfil,Ncl);
+  
+  % get filters coefficients
+  if Nfil>1
+    for ii=1:Nfil
+      A(ii,1:sza(ii)) = fil(ii).a;
+      B(ii,1:szb(ii)) = fil(ii).b;
+    end
+  else
+    A(1:sza) = fil.a;
+    B(1:szb) = fil.b;
+  end
+  % build Z matrix
+  Z = ones(size(A,2),numel(freq));
+  zz = exp(-2.*pi.*1i.*T.*freq);
+  for jj=2:size(Z,1)
+    Z(jj,:) = zz.^(jj-1);
+  end
+  
+  % get numerator and denominator
+  num = A*Z;
+  den = B*Z;
+  % get response
+  rspm = num./den;
+  if Nfil>1
+    switch lower(bank)
+      case 'parallel'
+        rsp = sum(rspm,1);
+      case 'serial'
+        rsp = prod(rspm,1);
+      otherwise
+        rsp = rspm;
+    end
+  else
+    rsp = rspm;
+  end
+
+
+end
\ No newline at end of file