[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: need help filtering a signal
From: |
Paul Kienzle |
Subject: |
Re: need help filtering a signal |
Date: |
Mon, 10 Dec 2001 13:14:23 -0500 |
Yes, it is easy to code in Octave. I have a version in octave-forge,
so no need for you to do it.
Paul Kienzle
address@hidden
On Sat, Dec 08, 2001 at 12:11:48PM -0500, Vincent Stanford wrote:
> 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
> -------------------------------------------------------------
>
-------------------------------------------------------------
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
-------------------------------------------------------------