利用窗函数法设计FIR滤波器

VibInfo

主程序一

1. C----------------------------------------------------------------------
   
2. c  main program H1DEFIR1: to test subroutine DEFIR1
   
3. C  To design FIR low-pass filter by Fourier series method.
   
4. c  Please link subroutine DEFIR1,FIRRES,AMPRES,PHARES,UNWRAP,WINDOW
   
5. C----------------------------------------------------------------------
   
6.         dimension b(0:20),w(0:20),amp(0:255),phase(0:255)
   
7.         complex h(0:255)
   
8.         data lb/20/,n/256/,iband/1/,iwindow/5/,iamp/0/
   
9.         data fl/0.125/,fs/1.0/
   
10.         call defir1(lb,iband,fl,fh,fs,iwindow,b,w,ierror)
    
11.         write(*,*)'   ierror=',ierror
    
12.         if(ierror.ne.0)stop
    
13.         write(*,3)(k,b(k),k=0,lb)
    
14. 3       format('                b(',i2,')=',f12.7)
    
15.         call firres(b,lb,n,h)
    
16.         call ampres(h,amp,n,fs,iamp)
    
17.         call phares(h,phase,n,fs)
    
18.         stop
    
19.         end

复制代码

主程序二：

1. C----------------------------------------------------------------------
   
2. c  main program H4DEFIR1: to test subroutine DEFIR1
   
3. C  To design FIR band-stop filter by Fourier series method.
   
4. c  Please link subroutine DEFIR1,FIRRES,AMPRES,PHARES,UNWRAP,WINDOW
   
5. C----------------------------------------------------------------------
   
6.         dimension b(0:20),w(0:20),amp(0:255),phase(0:255)
   
7.         complex h(0:255)
   
8.         data fl/0.2/,fh/0.3/,fs/1.0/
   
9.         data lb/20/,n/256/,IBAND/4/,IWindow/5/,iamp/0/
   
10.         call defir1(lb,iband,fl,fh,fs,iwindow,b,w,ierror)
    
11.         write(*,*)'   ierror=',ierror
    
12.         if(ierror.ne.0)stop
    
13.         write(*,3)(k,b(k),k=0,lb)
    
14. 3       format('                b(',i2,')=',f12.7)
    
15.         call firres(b,lb,n,h)
    
16.         call ampres(h,amp,n,fs,iamp)
    
17.         call phares(h,phase,n,fs)
    
18.         stop
    
19.         end

复制代码

子程序：

1.       subroutine defir1(l,iband,fl,fh,fs,iwindow,b,w,ierror)
   
2. c-----------------------------------------------------------------------
   
3. C  Subroutine DEFIR1: To Design FIR Filter By Windowed Fourier Series.
   
4. C  Fl:low cut-off frequency. Fh:high cut-off(For BP,BS). fl,fh,fs in Hz
   
5. C  Digital filter coefficients are returned in b(l)
   
6. c                H(z)=b(0)+b(1)z^(-1)+ ... +b(L)z^(-L)
   
7. c  Input parameters:
   
8. c   L+1  : the length of FIR filter. (L+1) must be odd .
   
9. c   iband:  iband=1 low  pass FIR filter.
   
10. c           iband=2 high pass FIR filter.
    
11. c           iband=3 band pass FIR filter.
    
12. c           iband=4 band stop FIR filter.
    
13. c   fln,fhn: the edge frequency desired, in normalized form, 0.<=fln,fhn<=.5
    
14. c          |---        |   ---       |   ---           |--      --
    
15. c          |   |       |  |          |  |   |          |  |    |
    
16. c          |   |       |  |          |  |   |          |  |    |
    
17. c        --|------    -|--------    -|-----------    --|--------------
    
18. c          0  fln      0  fln        0 fln  fhn        0  fln  fhn
    
19. c             fhn=*       fhn=*
    
20. c        iband=1 LP    iband=2 HP      iband=3 BP        iband=4 BS
    
21. c  iwindow: window type desired.
    
22. c           iwindow=1: rectangular window ,   =2: triangular window ,
    
23. C                  =3: cosin window ,         =4: Hanning window ,
    
24. C                  =5: Hamming window ,       =6: Blackman window ,
    
25. c                  =7: Papoulis window .
    
26. c        w: L dimensioned real work array.
    
27. c  Output parameters:
    
28. c   b: L+1 dimensioned real array.the result is in b(0) to b(L).
    
29. C   ierror: ierror=0:      no errors detected
    
30. C                 =1:      invalid length  (L<=0)
    
31. C                 =2:      invalid window type
    
32. C                 =3:      invalid filter type
    
33. C                 =4:      invalid cut-off frequency
    
34. c       From Ref. [5] in chapter 2 .
    
35. c                                      in Chapter 8
    
36. c-----------------------------------------------------------------------
    
37.         dimension b(0:l),w(0:l)
    
38.         pi=4.*atan(1.)
    
39.         fln=fl/fs
    
40.         fhn=fh/fs
    
41.         do 1 i=0,l
    
42.            b(i)=0.
    
43. 1       continue
    
44.         ierror=0
    
45.         l=l+1
    
46.         dly=float(l)/2.
    
47.         lim=l/2
    
48.         if(dly.eq.float(lim)) ierror=1
    
49.         if(iwindow.lt.1.or.iwindow.gt.7) ierror=2
    
50.         if(iband.lt.1.or.iband.gt.4) ierror=3
    
51.         if(fln.le.0..or.fln.gt.0.5) ierror=4
    
52.         if(iband.ge.3.and.fln.ge.fhn) ierror=4
    
53.         if(iband.ge.3.and.fhn.ge.0.5) ierror=4
    
54.         if(ierror.ne.0) return
    
55.         dly=dly-.5
    
56.         lim=lim-1
    
57.         wc1=2.*pi*fln
    
58.         if(iband.ge.3) wc2=2.*pi*fhn
    
59.         call window(w,l,iwindow,ierron)
    
60.         l=l-1
    
61.         go TO (5,10,15,20) iband
    
62. c----------- lowpass FIR filter design --------------------------------
    
63. 5       continue
    
64.         do 6 i=0,lim
    
65.            s=i-dly
    
66.            b(i)=((sin(wc1*s))/(pi*s))
    
67.            b(i)=b(i)*w(i)
    
68.            b(l-i)=b(i)
    
69. 6       continue
    
70.         b((l)/2)=wc1/pi
    
71.         return
    
72. c--------- highpass FIR filter design ---------------------------------
    
73. 10      continue
    
74.         do 11 i=0,lim
    
75.            s=i-dly
    
76.            b(i)=((sin(pi*s)-sin(wc1*s))/(pi*s))
    
77.            b(i)=b(i)*w(i)
    
78.            b(l-i)=b(i)
    
79. 11      continue
    
80.         b((l)/2)=1.-wc1/pi
    
81.         return
    
82. c-------- bandpass FIR filter design ----------------------------------
    
83. 15      continue
    
84.         do 16 i=0,lim
    
85.            s=i-dly
    
86.            b(i)=((sin(wc2*s)-sin(wc1*s))/(pi*s))
    
87.            b(i)=b(i)*w(i)
    
88.            b(l-i)=b(i)
    
89. 16      continue
    
90.         b((l)/2)=(wc2-wc1)/pi
    
91.         return
    
92. c-------- bandstop FIR filter design ----------------------------------
    
93. 20      continue
    
94.         do 21 i=0,lim
    
95.            s=i-dly
    
96.            b(i)=((sin(wc1*s)+sin(pi*s)-sin(wc2*s))/(pi*s))
    
97.            b(i)=b(i)*w(i)
    
98.            b(l-i)=b(i)
    
99. 21      continue
    
100.         b((l)/2)=(wc1+pi-wc2)/pi
     
101.         return
     
102.         end

复制代码