Re: need help filtering a signal

[Vincent Stanford]

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