[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: need help filtering a signal
From: |
Vincent Stanford |
Subject: |
Re: need help filtering a signal |
Date: |
Sat, 08 Dec 2001 12:11:48 -0500 |
User-agent: |
Mozilla/5.0 (X11; U; Linux i686; en-US; rv:0.9.2.1) Gecko/20010901 |
Dear Philipp,
I am attaching a fortran code that I wrote a while back. It produces a
direct form Butterworth bandpass filter. It does a little trig to get
an appropriate set of poles and zeros and then multiplies them to get
the transfer function weights as polynomial coefficients in z. A filter
with a moderate order (like four to eight) will give you a very nice
bandpass filter that requires a minimum of arithmetic and has really
fantastic stop band attenuation, as well as a flat passband. Being an
IIR it does not have linear phase, but it is very cheap to run, so it is
well suited to real time operation. It just takes an order (must be
even) and the upper and lower half power points in normalized radian
frequency (nyquist=0.5) on stdin. It outputs the filter weights
(transfer function) to a stdout. It also outputs a frequency response
table to fort12. The transfer function is given as
X(z) 1 + x1*z^-1 + x2*z^-2 + x3*z^-3 + ....
H(z)= -------- = -----------------------------------------------
Y(z) 1 + y1*z^-1 + y2*z^-2 + y3*z^-3 + ....
so don't forget about the sign changes when you write the filter as:
y(t)=-(y1*y(t-1) + y2*y(t-2) + y3*y(t-3) + ....) + x(t) + x1*x(t-1) +
x2*x(t-2) + x3*x(t-3) + ....
the fortran code follows. Maybe you could translate it to octave and
give it to the community ;-). I think that there is a polynomial
multiply in octave so the code could be very simple in octave.
Vince Stanford
C=======================================================================
C Name butterworthBandpass
C-----------------------------------------------------------------------
C
C This program computes the poles and zeroes of a lowpass butterworth
C and a highpass butterworth filter with appropriate cuttoff frequencies
C for the bandpass desired. The poles and zeroes are then multiplied
C together to form the coefficents of a polynomial in z, which can
C be used as a direct form I realization of the band pass. The program
C requires that a filter of even order > 2 be generated. Apply the z=-z
C lowpass to highpass transform
C
C Author: Vince Stanford, NIST Smart Space Project
C
C Date: February 26, 2001
C
C=======================================================================
implicit none
C
C allows for fifty poles and fifty zeroes
C
complex center
complex point
complex poles(50)
complex zeroes(50)
complex xWeights(51)
complex yWeights(51)
complex scratch(51)
integer i
integer j
integer magPoints
integer dataOut /12/
integer order /0/
real centerFreq /0.0/
real gainCenter /0.0/
real gain /0.0/
real hiCutFreq /0.0/
real lowCutFreq /0.0/
real pi /3.141592654/
C
C linux std logical unit numbers
C
integer stdin
integer stdout
integer stderr
stdin=5
stdout=6
stderr=0
do i=1,50
poles(i)=(0.0,0.0)
zeroes(i)=(0.0,0.0)
xWeights(i)=(0.0,0.0)
yWeights(i)=(0.0,0.0)
scratch(i)=(0.0,0.0)
end do
C
C Read and parse low and high frequencies for errors.
C these must be in normalized radian frequency
C
write(stdout,*) '---------------------------'
write(stdout,*) 'Butterworth Bandpass filter'
write(stdout,*) '---------------------------'
write(stdout,*) 'filter order?'
read(stdin,*) order
if(order.gt.50) then
stop 'order > 50.'
elseif (order.lt.2) then
stop 'order < 2. '
elseif ((order/2)*2.ne.order) then
stop 'order not even.'
endif
write(stdout,*) 'lower half power point?'
read(stdin,*) lowCutFreq
write(stdout,*) 'higher half power point?'
read(stdin,*) hiCutFreq
if (lowCutFreq.gt.0.5) then
stop 'low cut freq. > 0.5. '
elseif (lowCutFreq.le.0.0) then
stop 'cut freq must be > 0.0. '
endif
if (hiCutFreq.gt.0.5) then
stop 'low cut freq. > 0.5. '
elseif (hiCutFreq.le.0.0) then
stop 'cut freq must be > 0.0. '
endif
if(lowCutFreq.gt.hiCutFreq) then
stop 'low cut freq > high cut freq'
endif
centerFreq=(lowCutFreq+hiCutFreq)/2.0
C
C compute poles and zeroes of lowpass filter
C
call butter(order/2,lowCutFreq,zeroes,poles)
hiCutFreq=0.5-hiCutFreq
call butter(order/2,hiCutFreq,zeroes(order/2+1),poles(order/2+1))
C
C apply lowpass to highpass tranform
C
do i=order/2+1,order
zeroes(i)=-zeroes(i)
poles(i)=-poles(i)
end do
write(stdout,*) '---------------'
do i=1, order
write(stdout,*) 'zero ',i,' at ',zeroes(i)
end do
write(stdout,*) '---------------'
do i=1, order
write(stdout,*) 'pole ',i,' at ',poles(i)
end do
C
C compute coefficients from the poles and zeroes
C
C zeroes
C
xWeights(1)=-zeroes(1)
xWeights(2)=(1.0,0.0)
if (order.gt.1) then
do i=2,order
scratch(1)=(0.0,0.0)
do j=1,i
scratch(j+1)=xWeights(j)
end do
do j=1,i+1
xWeights(j)=-zeroes(i)*xWeights(j)
end do
do j=1,i+1
xWeights(j)=xWeights(j)+scratch(j)
end do
end do
endif
write(stdout,*) '---------------'
do i=1,order+1
write(stdout,*) 'zero coefficient(',i,')=',real(xWeights(i))
end do
C
C poles
C
yWeights(1)=-poles(1)
yWeights(2)=(1.0,0.0)
if (order.gt.1) then
do i=2,order
scratch(1)=(0.0,0.0)
do j=1,i
scratch(j+1)=yWeights(j)
end do
do j=1,i+1
yWeights(j)=-poles(i)*yWeights(j)
end do
do j=1,i+1
yWeights(j)=yWeights(j)+scratch(j)
end do
end do
endif
write(stdout,*) '---------------'
do i=1,order+1
write(stdout,*) 'pole coefficent(',i,')=',real(yWeights(i))
end do
C
C compute gain at center due to poles and zeroes for
C unity gain normalization of the weights
C
gainCenter=1.0
center=exp(cmplx(0.0,pi*centerFreq))
do i=1,order
gainCenter=gainCenter*
& abs(center-zeroes(i))/abs(center-poles(i))
end do
write(stdout,*) '------------------------'
write(stdout,*) 'system gain=',gainCenter
write(stdout,*) '------------------------'
magPoints=1000
do j=1,magPoints
gain=1.0
point=exp(cmplx(0.0,pi*float(j)/float(magPoints)))
do i=1,order
gain=gain*
& abs(point-zeroes(i))/abs(point-poles(i))
end do
write(dataOut,*) gain/gainCenter
end do
stop 'normal termination: butterworth. '
end
C-------------------------------------------------------------
C
C Subroutine computes the poles and zeroes of a
C butterworth lowpass digital filter.
C
C Inputs:
C
C n = order of the filter required
C fc = required cutoff normalized radian frequency
C
C Outputs:
C
C alpha = complex array containing the transfer
C function zeroes in its first n locations.
C (all n zeroes lie at z=-1)
C
C beta = complex array containing the transfer
C function poles in its first n locations.
C the complex conjugate pairs of poles
C adjacent locations; if n is odd the real
C pole is in location 1
C---------------------------------------------------------------
subroutine butter(n,fc,alpha,beta)
complex alpha(n)
complex beta(n)
real pi /3.141592654/
wc=pi*fc
tan2=2.0*sin(wc)/cos(wc)
tansq=0.25*tan2**2
if (n.eq.1) go to 2
in=mod(n,2)
n1=n+in
n2=(3*n+in)/2-1
do m=n1,n2
a=pi*float(2*m+1-in)/float(2*n)
anum=1.0-tan2*cos(a)+tansq
u=(1.0-tansq)/anum
v=tan2*sin(a)/anum
i=(n2-m)*2+1
beta(i+in)=cmplx(u,v)
beta(i+in+1)=cmplx(u,-v)
end do
if(in) 3,3,2
2 beta(1)=cmplx(((1.0-tansq)/(1.0+tan2+tansq)),0.0)
3 do i=1,n
alpha(i)=(-1.0,0.0)
end do
return
end
Philipp Schwaha wrote:
hi!
what i need to do is fiter a signal to remove noise. what i wanted to try
first was a fft then cut all frequencies except the the important ones (e.g
100Hz, 200Hz ...), but when i did the fft i did not know how to do this,
since there are no frequencies associated with the values the function fft
returns.
this is the first time i try to do something using an fft.
then i discovered the function fftfilter, but i could not make it produce the
results i wanted.
how can i remove all frequencies from a certain signal, except a few selected
ones?
thanks
philipp
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------