From help-octave-request at bevo dot che dot wisc dot edu Mon Dec 10 12:17:17 2001 Subject: Re: need help filtering a signal From: Paul Kienzle To: Vincent Stanford , Philipp Schwaha Cc: help-octave at bevo dot che dot wisc dot edu 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 pkienzle at users dot sf dot net 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 -------------------------------------------------------------