Mercurial > hg > ltpda
comparison m-toolbox/classes/+utils/@math/mtxiirresp.m @ 0:f0afece42f48
Import.
| author | Daniele Nicolodi <nicolodi@science.unitn.it> |
|---|---|
| date | Wed, 23 Nov 2011 19:22:13 +0100 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:f0afece42f48 |
|---|---|
| 1 % MTXIIRRESP calculate iir resp by matrix product | |
| 2 % | |
| 3 % INPUT: | |
| 4 % fil: a miir filter or a vector of filoters | |
| 5 % freq: the frequency vector | |
| 6 % fs: sampling frequency | |
| 7 % bank: filterbank type parameter ('parallel', 'serial' or 'none'). | |
| 8 % If it is left empty, it provide a vector of responses at the | |
| 9 % output. | |
| 10 % | |
| 11 % | |
| 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 13 function rsp = mtxiirresp(fil,freq,fs,bank) | |
| 14 | |
| 15 T = 1/fs; | |
| 16 | |
| 17 Nfil = numel(fil); | |
| 18 | |
| 19 if Nfil>1 | |
| 20 % check for different orders | |
| 21 sza = zeros(Nfil,1); | |
| 22 szb = zeros(Nfil,1); | |
| 23 for jj=1:Nfil | |
| 24 sza(jj) = numel(fil(jj).a); | |
| 25 szb(jj) = numel(fil(jj).b); | |
| 26 end | |
| 27 Ncl = max([sza;szb]); | |
| 28 else | |
| 29 sza = numel(fil.a); | |
| 30 szb = numel(fil.b); | |
| 31 Ncl = max([sza szb]); | |
| 32 end | |
| 33 | |
| 34 % init A and B | |
| 35 A = zeros(Nfil,Ncl); | |
| 36 B = zeros(Nfil,Ncl); | |
| 37 | |
| 38 % get filters coefficients | |
| 39 if Nfil>1 | |
| 40 for ii=1:Nfil | |
| 41 A(ii,1:sza(ii)) = fil(ii).a; | |
| 42 B(ii,1:szb(ii)) = fil(ii).b; | |
| 43 end | |
| 44 else | |
| 45 A(1:sza) = fil.a; | |
| 46 B(1:szb) = fil.b; | |
| 47 end | |
| 48 % build Z matrix | |
| 49 Z = ones(size(A,2),numel(freq)); | |
| 50 zz = exp(-2.*pi.*1i.*T.*freq); | |
| 51 for jj=2:size(Z,1) | |
| 52 Z(jj,:) = zz.^(jj-1); | |
| 53 end | |
| 54 | |
| 55 % get numerator and denominator | |
| 56 num = A*Z; | |
| 57 den = B*Z; | |
| 58 % get response | |
| 59 rspm = num./den; | |
| 60 if Nfil>1 | |
| 61 switch lower(bank) | |
| 62 case 'parallel' | |
| 63 rsp = sum(rspm,1); | |
| 64 case 'serial' | |
| 65 rsp = prod(rspm,1); | |
| 66 otherwise | |
| 67 rsp = rspm; | |
| 68 end | |
| 69 else | |
| 70 rsp = rspm; | |
| 71 end | |
| 72 | |
| 73 | |
| 74 end |