将模拟滤波器转变为数字滤波器

VibInfo

1.         subroutine aftodf(d,c,ln,iband,fln,fhn,b,a,ierror)
   
2. c----------------------------------------------------------------------
   
3. c Routine AFTODF: To convert normalized LP analog H(s) to digital H(z).
   
4. c   H(s)=D(s)/C(s),H(z)=B(z)/A(z).Filter order l is computed internally.
   
5. c   LN specifies coefficient array size. WORK(0:LN,0:LN) is a work array.
   
6. c   IF   IBAND=1:    lowpass   fln=normalized cutoff frequency
   
7. c             =2:    highpass  fln=normalized cutoff frequency
   
8. c             =3:    bandpass  fln=low  cutoff frequency
   
9. c                              fhn=high cutoff frequency
   
10. c             =4:    bandstop  fln=low  cutoff frequency
    
11. c                              fhn=high cutoff frequency
    
12. c   IF  IERROR=0:    no errors detected
    
13. c              1:    all zero transfer function
    
14. c              2:    biline: invalid transfer function
    
15. c              3:    filter order exceeds array size
    
16. c              4:    invalid filter type parameter (IBAND)
    
17. c              5:    invalid cutoff frequency
    
18. c       From Ref. [5] of Chapter 2 .      in chapter 7
    
19. c-----------------------------------------------------------------------
    
20.         dimension work(0:4,0:4),d(0:4),c(0:4),b(0:4),a(0:4)
    
21.         pi=4.*atan(1.)
    
22.         ierror=0
    
23.         if(iband.lt.1.or.iband.gt.4) ierror=4
    
24.         if(fln.le.0..or.fln.gt.0.5) ierror=5
    
25.         if(iband.ge.3.and.fln.ge.fhn) ierror=5
    
26.         if(iband.ge.3.and.fhn.gt.0.5) ierror=5
    
27.         if(ierror.ne.0) return
    
28.         do 1 i=ln,0,-1
    
29.            if(c(i).ne.0..or.d(i).ne.0.) go to 2
    
30. 1       continue
    
31.         ierror=1
    
32.         return
    
33. 2       m=i
    
34.         w1=tan(pi*fln)
    
35.         l=m
    
36.         if(iband.le.2) go to 3
    
37.         l=2*m
    
38.         w2=tan(pi*fhn)
    
39.         w=w2-w1
    
40.         w02=w1*w2
    
41. 3       continue
    
42.         ierror=3
    
43.         if(l.gt.ln) return
    
44.         go to (30,20,40,20) iband
    
45. c-------- substitution of 1/s to generate highpass (hp,bs) ------------
    
46. 20      continue
    
47.         do 25 mm=0,m/2
    
48.            tmp=d(mm)
    
49.            d(mm)=d(m-mm)
    
50.            d(m-mm)=tmp
    
51.            tmp=c(mm)
    
52.            c(mm)=c(m-mm)
    
53.            c(m-mm)=tmp
    
54. 25      continue
    
55.         if(iband.eq.4) go to 40
    
56. c-------- scaling s/w1 for lowpass,highpass ---------------------------
    
57. 30      continue
    
58.         do 35 mm=0,m
    
59.            d(mm)=d(mm)/(w1**mm)
    
60.            c(mm)=c(mm)/(w1**mm)
    
61. 35      continue
    
62.         go to 100
    
63. c-------- substitution of (s**2+w0**2)/(w*s)  bandpass,bandstop -------
    
64. 40      continue
    
65.         do 45 ll=1,l+1
    
66.            work(ll,1)=0.
    
67.            work(ll,2)=0.
    
68. 45      continue
    
69.         do 52 mm=0,m
    
70.            tmpd=d(mm)*(w**(m-mm))
    
71.            tmpc=c(mm)*(w**(m-mm))
    
72.            do 50 k=0,mm
    
73.               ls=m+mm-2*k
    
74.               tmp=spbfct(mm,mm)/(spbfct(k,k)*spbfct(mm-k,mm-k))
    
75.               work(ls+1,1)=work(ls+1,1)+tmpd*(w02**k)*tmp
    
76.               work(ls+1,2)=work(ls+1,2)+tmpc*(w02**k)*tmp
    
77. 50         continue
    
78. 52      continue
    
79.         do 55 ll=0,l
    
80.            d(ll)=work(ll+1,1)
    
81.            c(ll)=work(ll+1,2)
    
82. 55      continue
    
83. c---------- substitute (z-1)/(z+1) ------------------------------------
    
84. 100     continue
    
85.         call biline(work,d,c,ln,b,a,ierror)
    
86.         if(ierror.eq.0) return
    
87.         write(*,*)'   stop at routine biline,ierror=',ierror
    
88.         return
    
89.         end
    
90. c
    
91.         function spbfct(i1,i2)
    
92. c-------- generates (i1)!/(i1-i2)!=i1*(i1-1)*...*(i1-i2+1). -----------
    
93. c-------- note: 0!=1 and spbfct(i,i)=spbfct(i,i-1)=i!.      -----------
    
94.         spbfct=0.
    
95.         if(i1.lt.0.or.i2.lt.0.or.i2.gt.i1) return
    
96.         spbfct=1.
    
97.         if(i2.eq.0) return
    
98.         do 31 i=i1,i1-i2+1,-1
    
99.            spbfct=spbfct*i
    
100. 31      continue
     
101.         return
     
102.         end

复制代码

[ 本帖最后由 VibInfo 于 2006-8-6 07:27 编辑 ]