实现双线性Z变换

VibInfo

1.         subroutine biline(work,d,c,ln,b,a,ierror)
   
2. c----------------------------------------------------------------------
   
3. c Routine BILINE: To convert analog H(S) to digital H(Z) via bilinear
   
4. c                 transformation. H(S)=D(S)/C(S), H(Z)=B(Z)/A(Z)
   
5. c LN specifies the length of the coefficient arrays and filter order L
   
6. c is computed internally.   WORK is an internal array (2D)
   
7. c   IF  IERROR=0:    no errors detected in transformation
   
8. c             =1:    all zero transfer function
   
9. c             =2:    invalid transfer function; y(k) coef=0
   
10. c       From Ref. [5] of Chapter 2 .       in chapter 7
    
11. c----------------------------------------------------------------------
    
12.         dimension work(0:4,0:4),d(0:4),c(0:4),b(0:4),a(0:4)
    
13.         do 10 i=ln,0,-1
    
14.            if(c(i).ne.0..or.d(i).ne.0.)go to 20
    
15. 10      continue
    
16.         ierror=1
    
17.         return
    
18. 20      l=i
    
19.         do 30 J=0,L
    
20.            work(0,j)=1.
    
21. 30       continue
    
22.         tmp=1.
    
23.         do 40 i=1,l
    
24.            tmp=tmp*float(l-i+1)/float(i)
    
25.            work(i,0)=tmp
    
26. 40       continue
    
27.         do 60 i=1,l
    
28.            do 50 j=1,l
    
29.               work(i,j)=work(i,j-1)-work(i-1,j)-work(i-1,j-1)
    
30. 50         continue
    
31. 60      continue
    
32.         do 80 i=l,0,-1
    
33.            b(i)=0.
    
34.            atmp=0.
    
35.            do 70 j=0,l
    
36.               b(i)=b(i)+work(i,j)*d(j)
    
37.               atmp=atmp+work(i,j)*c(j)
    
38. 70          continue
    
39.            scale=atmp
    
40.            if(i.ne.0) a(i)=atmp
    
41. 80      continue
    
42.         ierror=2
    
43.         if(scale.eq.0.) return
    
44.         b(0)=b(0)/scale
    
45.         do 90 i=1,l
    
46.            b(i)=b(i)/scale
    
47.            a(i)=a(i)/scale
    
48. 90      continue
    
49.         do 100 i=l+1,ln
    
50.            b(i)=0.0
    
51.            a(i)=0.0
    
52. 100     continue
    
53.         ierror=0
    
54.         return
    
55.         end

复制代码