'*************************************************************** '* Program to demonstrate the Butterworth filter * '* (removing frequencies greater then Fc) * '* * '* Basic version by J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ----------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Input data file (tfft.dat): * '* * '* (The test signal contains 3 frequencies: 50, 250, 500 hz) * '* * '* 1024 * '* 0.00000000000000E+0000 0.00000000000000E+0000 * '* 1.95503421309917E-0004 3.53914399999776E+0001 * '* 3.91006842619834E-0004 5.95684899999760E+0001 * '* 5.86510263929974E-0004 6.54621699999552E+0001 * '* 7.82013685239669E-0004 5.24038399999845E+0001 * '* ... ... * '* 1.98826979472187E-0001 2.77372500000183E-0001 * '* 1.99022482893497E-0001 -2.43361500000174E+0000 * '* 1.99217986314807E-0001 -4.84236799999780E+0000 * '* 1.99413489736116E-0001 -6.02247899999929E+0000 * '* 1.99608993157426E-0001 -5.45615399999951E+0000 * '* 1.99804496578736E-0001 -3.22824200000105E+0000 * '* 2.00000000000045E-0001 -2.96010699999982E-0003 * '* * '* Output file (tfilter.lst): * '* * '* Time Input Signal Filtered Signal * '* ------------------------------------------------------- * '* 0.0000 0.000000 0.000000 * '* 0.0002 35.391440 0.158005 * '* 0.0004 59.568490 1.278083 * '* 0.0006 65.462170 5.003994 * '* 0.0008 52.403840 12.809453 * '* 0.0010 26.229620 24.336649 * '* ... ... ... * '* 0.1988 0.277373 -2.224546 * '* 0.1990 -2.433615 -0.802179 * '* 0.1992 -4.842368 0.588941 * '* 0.1994 -6.022479 1.320860 * '* 0.1996 -5.456154 1.036090 * '* 0.1998 -3.228242 -0.213007 * '* 0.2000 -0.002960 -1.978450 * '* * '* Note: The cut off frequency Fc is 500 hz, the filter order * '* is 4 (A Fourier analysis of the filtered signal shows * '* that frequencies greater than Fc are satisfactorily * '* removed but with a slight time shift of the signal). * '*************************************************************** defint i-n defdbl a-h,o-z MACHEPS = 1e-12 'i,iorder,ndata : INTEGER 'dt,Fc,t,tbegin,tend : DOUBLE dim signal(1024),filtered(1024) dim C(5,10), D(2,10) 'NSections: INTEGER 'Tg,Xdc : DOUBLE cls 'open input and output file open "tfft.dat" for input as #1 open "tfilter.lst" for output as #2 'read number of input signal points in input file input #1, ndata 'read ndata couples T(i), Y(i) in input data file for i=1 to ndata input #1, temp, signal(i) if (i=1) then tbegin=temp if (i=ndata) then tend=temp next i close #1 'calculate sampling increment dt of input signal dt = (tend-tbegin)/(ndata-1) 'input cut off frequencys and order of filter (1 to 4) Fc=500# : iorder=4 : Xdc=0# 'no constant value 'call filtering procedures '1. Calculate the filter coefficients 'Butterworth(Fc,dt,iorder,C,NSections,Tg); gosub 3000 '2. Initialize filter memory 'Init(Xdc,C,NSections,D); gosub 2000 '3. Recursively call Butterworth filter FOR i=1 TO ndata 'Filter(filtered(i),signal(i),NSections,C,D); gosub 1000 NEXT i 'print input and filtered signals to output file t=tbegin print #2," Time Input Signal Filtered Signal " print #2," -------------------------------------------------------" for i=1 to ndata print #2,using " ##.#### ####.###### ###.######"; t; signal(i); filtered(i) t=t+dt next i close #2 print print " Results in file tfilter.lst." END 'of main program '*********************************************************************** '* Filtering a signal F(t) by Butterworth method * '* (removing frequencies greater then Fc) * '* ------------------------------------------------------------------- * '* Calling mode: Filter(Xs,Xd,NSections,C,D); * '* ------------------------------------------------------------------- * '* INPUTS: * '* ------- * '* signal(i): current value of input signal (double) * '* NSections: Number of required 2nd order sections (integer) * '* = n/2 for n even * '* = (n+1)/2 for n odd * '* C........: Table[1..5,1..NSections] of filter coefficients * '* calculated previously by BUTTERWORTH subroutine * '* D........: Table[1..2,1..NSections] of coefficients defining * '* the filter memory, initialized by INIT subroutine * '* ------------------------------------------------------------------- * '* OUTPUTS: * '* ------- * '* D..........: Table updated after the call to Filter subroutine * '* filtered(i): current value of filtered signal (double) * '* ------------------------------------------------------------------- * '* Reference * '* --------- * '* "Lawrence R.Rabiner et Bernard Gold * '* Theory and application of digital processing. * '* Prentice Hall Inc., EnglewoodclIFfs,NEW JERSEY,1975." * '* * '* Basic version by J-P Moreau, Paris * '* from Fortran version by J-P Dumont / Tuan Dang Trong * '*********************************************************************** 1000 'PROCEDURE Filter x=signal(i) FOR ii=1 TO NSections xerr=x+C(1,ii)*D(1,ii)+C(2,ii)*D(2,ii) y=C(5,ii)*(xerr +C(3,ii)*D(1,ii)+C(4,ii)*D(2,ii)) D(2,ii)=D(1,ii) D(1,ii)=xerr x=y NEXT ii filtered(i)=x RETURN '************************************************************************** '* INIT FILTER PROCEDURE * '* ---------------------------------------------------------------------- * '* The filter response is initialized to stationnary value for a constant * '* input signal value. * '* * '* Calling mode: INIT(Xdc,C,NSections,D); * '* ---------------------------------------------------------------------- * '* INPUTS: * '* ------ * '* Xdc......: constant input value (double) * '* C........: Table[1..5,1..NSections] of filter coefficients * '* calculated previously by BUTTERWORTH subroutine * '* NSections: Number of required 2nd order sections (integer) * '* = n/2 for n even (n=order of filter) * '* = (n+1)/2 for n odd * '* calculated previously by BUTTERWORTH subroutine * '* ---------------------------------------------------------------------- * '* OUTPUTS: * '* ------- * '* D........: Table[1..2,1..NSections] of coefficients defining * '* the filter memory, initialized by INIT subroutine * '************************************************************************** 2000 'PROCEDURE Init dc=Xdc FOR j=1 TO NSections D(2,j)=dc/(1#-C(1,j)-C(2,j)) D(1,j)=D(2,j) Csum=0# FOR ii=1 TO 4 Csum=Csum + C(ii,j) NEXT ii dc=C(5,j)*(dc+D(2,j)*Csum) NEXT j RETURN '********************************************************************** '* Calculates the Butterworth filter coefficients * '* ------------------------------------------------------------------ * '* Calling mode: Butterworth(Fc,Ts,n,C,NSections,Tg); * '* ------------------------------------------------------------------ * '* INPUTS: * '* ------ * '* Fc.......: Cut off frequency * '* dt.......: Sampling time of input signal * '* N........: Order of filter (1 to 4) * '* ------------------------------------------------------------------ * '* OUTPUTS: * '* ------- * '* C........: Table[1..5,1..NSections] of filter coefficients * '* calculated previously by BUTTERWORTH subroutine * '* NSections: Number of required 2nd order sections (integer) * '* = N/2 for N even * '* = (N+1)/2 for N odd * '* Tg.......: Group delay in seconds * '*********************************************************************} 3000 'PROCEDURE Butterworth Zero = 0# ONE = 1# TWO = 2# HALF = 0.5# Pi = 3.1415926535# N=iorder arg=Pi*dt*Fc If abs(arg) > 2#*Pi THEN m=INT(arg/2#/Pi) arg=arg-(m*2#*Pi) END IF Omega= Tan(arg) OmegaSq=Omega*Omega xn=N If xn/2 = INT(xn/2) THEN xmodn=Zero : temp=HALF else xmodn=HALF : temp=Zero end if xns2 = N/2 NSections=INT(xns2+xmodn) Tg=Zero If (N>1) THEN FOR i=1 TO xns2 Rep=Omega*COS(Pi*(i-temp)/N) Tg=Tg+dt*Rep/OmegaSq W0=TWO*Rep W1=ONE + W0+OmegaSq C(1,i)=-TWO*(OmegaSq-ONE)/W1 C(2,i)=-(ONE-W0+OmegaSq)/W1 C(3,i)=TWO C(4,i)=ONE C(5,i)=OmegaSq/W1 NEXT i end if If temp = Zero THEN C(1,Nsections)=(ONE-Omega)/(ONE+Omega) C(2,NSections)= Zero C(3,NSections)= ONE C(4,NSections)= Zero C(5,NSections)= Omega/(ONE+Omega) Tg= Tg+dt/(TWO*Omega) END IF RETURN ' from Butterworth 'End of file tfilters.bas