Re: need help filtering a signal

[Paul Kienzle]

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
-------------------------------------------------------------