1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
|
//Prony's method for time-domain design of IIR Filters
//[b,a]=prony(h,nb,na)
//where b= coefficients of the numerator of the TF
// a=coefficients of the denominator of the TF
// h=impulse response of the digital filter
// nb=number of zeros
// na=number of poles
//Reference : T.W. Parks and C.S. Burrus, Digital Filter Design,
// John Wiley and Sons, 1987, p226.
//Author
//Debdeep Dey
//EXAMPLE 1:
//V=filter([1,1],[1,1,2],[1 zeros(1,31)]);//gives the impulse response of filter with the mentioned coefficients
//[b,a]=prony(V,1,2)
//
//OUTPUT:a = 1. 1. 2.//denominator coefficients
//b = 1. 1. //numerator coefficients
//EXAMPLE2:
//V=filter([1,2],[1,2,3,4],[1 zeros(1,31)]);
//[b,a]=prony(V,1,3)
//a = 1. 2. 3. 4. //denominato coefficients
//b = 1. 2. //numerator coefficients
function [b,a]=prony(h,nb,na)
K = length(h)-1;
M=double(nb);
N=double(na);
//zero-pad input if necessary
if K <= max(M,N) then
K = max(M,N)+1;
h(K+1) = 0;
end
c = h(1);
if c==0 //avoid division by zero
c=1;
end
H = toeplitz(h/c,[1 zeros(1,K)]);
//K+1 by N+1
if (K > N)
H(:,(N+2):(K+1)) = [];
end
//Partition H matrix
H1 = H(1:(M+1),:); //M+1 by N+1
h1 = H((M+2):(K+1),1); //K-M by 1
H2 = H((M+2):(K+1),2:(N+1)); //K-M by N
a = [1; -H2\h1].';
b = c*a*H1.';
endfunction
|
