From help-octave-request at bevo dot che dot wisc dot edu Sat Dec 8 11:08:19 2001 Subject: Re: need help filtering a signal From: Vincent Stanford To: Philipp Schwaha CC: help-octave at bevo dot che dot wisc dot edu Date: Sat, 08 Dec 2001 12:11:48 -0500 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 -------------------------------------------------------------