用频率抽样法设计FIR滤波器
主程序一C----------------------------------------------------------------------
cmain program H1DEFIR2: to test subroutine DEFIR2
CTo design FIR low-pass filter by frequency sampling method
cPlease link subroutine DEFIR2,FIRRES,AMPRES,PHARES,UNWRAP
C----------------------------------------------------------------------
dimension b(0:30),amp(0:255),phase(0:255)
complex h(0:255),bcmplx(0:30)
data l/30/,n/256/,iband/1/,iamp/0/
data fs/1.0/,fl/0.1/,trans/0.5/
call defir2(l,iband,fl,fh,bcmplx,trans,fs,ierror)
write(*,*)' ierror=',ierror
if(ierror.ne.0)stop
do 10 i=0,l
b(i)=real(bcmplx(i))
write(*,*)i,b(i)
10 continue
call firres(b,l,n,h)
call ampres(h,amp,n,fs,iamp)
call phares(h,phase,n,fs)
stop
end
主程序二
C----------------------------------------------------------------------
cmain program H4DEFIR2: to test subroutine DEFIR2
CTo design FIR band-stop filter by frequency sampling method
cPlease link subroutine DEFIR2,FIRRES,AMPRES,PHARES,UNWRAP
C----------------------------------------------------------------------
dimension b(0:40),amp(0:255),phase(0:255)
complex h(0:255),bcmplx(0:40)
data l/40/,n/256/,iband/4/,iamp/0/
data fs/1.0/,fl/0.2/,fh/0.3/,trans/0.5/
call defir2(l,iband,fl,fh,bcmplx,trans,fs,ierror)
write(*,*)' ierror=',ierror
if(ierror.ne.0)stop
do 10 i=0,l
b(i)=real(bcmplx(i))
write(*,*)i,b(i)
10 continue
open(3,file='h4hdfir2.dat',status='new')
do 20 k=0,l
write(3,*)k,b(k)
20 CONTINUE
close(3)
call firres(b,l,n,h)
call ampres(h,amp,n,fs,iamp)
call phares(h,phase,n,fs)
stop
end
子程序：
subroutine defir2(l,iband,fl,fh,b,trans,fs,ierror)
C-----------------------------------------------------------------------
CRoutine DEFIR2: To design FIR filter by frequency sampling method.
CFl:low cut-off frequency. Fh:high cut-off(For BP,BS). fl,fh,fs in Hz
c |--- | --- | --- |-- --
c | | || || | || |
c | | || || | || |
c --|------ -|-------- -|----------- --|--------------
c 0 fl 0fl 0 flfh 0 fl fh
C
C Digital filter coefficients are returned in b(l)
c H(z)=b(0)+b(1)z^(-1)+ ... +b(L)z^(-L)
cInput parameters:
c L+1: the length of FIR filter. L<200 and (L+1) must be the odd.
c iband:iband=1 lowpass FIR filter.
c iband=2 high pass FIR filter.
c iband=3 band pass FIR filter.
c iband=4 band stop FIR filter.
c trans: 0<= trans <1.0 , it is the transition point.
cOutput parameters:
c b: (L+1) dimensioned real array.the result is in b(0) to b(L).
c
c in Chapter 8
c-----------------------------------------------------------------------
complex b(0:l),h(0:200)
npass=0
pi=4.*atan(1.)
fln=fl/fs
fhn=fh/fs
do 5 i=0,l
h(i)=0.
5 continue
ierror=0
l=l+1
dly=float(l)/2.
lim=l/2
if(dly.eq.float(lim)) ierror=1
if(l.ge.201) ierror=2
if(iband.ne.4) goto 6
band=(fhn-fln)*l
if(band.ge.4.) goto 6
write(*,*)'Please increse the length L for Band-Stop FILTER'
ierror=3
6 if(iband.lt.1.or.iband.gt.4) ierror=4
if(fln.le.0..or.fln.gt.0.5) ierror=5
if(iband.ge.3.and.fln.ge.fhn) ierror=6
if(trans.lt.0.and.trans.ge.1.)ierror=7
if(ierror.ne.0) return
s=-(l-1)*pi/l
go to(10,20,30,10) iband
10 nlow=1
nhigh=fln*l-1
h(0)=1.0
do 15 i=nlow,nhigh
h(i)=cexp(cmplx(0.0,s*i))
h(l-i)=conjg(h(i))
15 continue
h(nhigh+1)=trans*cexp(cmplx(0.0,s*(nhigh+1)))
h(l-nhigh-1)=conjg(h(nhigh+1))
if(iband.eq.1)goto 100
nlow=fhn*l
goto 21
c
20 h(0)=0.0
nlow=fln*l
21 nhigh=lim
do 25 i=nlow,nhigh
h(i)=cexp(cmplx(0.0,s*i))
h(l-i)=conjg(h(i))
25 continue
h(nlow-1)=trans*cexp(cmplx(0.0,s*(nlow-1)))
h(l-nlow+1)=conjg(h(nlow-1))
goto 100
c
30 nlow=fln*l
nhigh=fhn*l
h(0)=0.0
do 35 i=nlow,nhigh
h(i)=cexp(cmplx(0.0,s*i))
h(l-i)=conjg(h(i))
35 continue
h(nhigh+1)=trans*cexp(cmplx(0.0,s*(nhigh+1)))
h(l-nhigh-1)=conjg(h(nhigh+1))
h(nlow-1)=trans*cexp(cmplx(0.0,s*(nlow-1)))
h(l-nlow+1)=conjg(h(nlow-1))
c
100 call cmpdft(h,b,l,1)
l=l-1
return
end
[ 本帖最后由 VibInfo 于 2006-8-8 07:17 编辑 ] 谢谢
页:
[1]