- Email: [email protected]
A simple program for a phaseless recursive digital filter
A SIMPLE PROGRAM FOR A PHASELESS RECURSIVE DIGITAL FILTER P.B. Pynsent
and R. Hanka*
ABSTRACT A brief comparison of several types of digital filters is made and a simple subroutine implementing a recursive digital phaseless low and high pass filter is described. The subroutine, KeyworL:
computer programs, digital filter, phaseless filtering
INTRODUCTION In biomedical engineering it is often found that signals are contaminated by the presence of noise. If the frequency spectra of the signal and noise are separable then it is possible to supress the noise with a suitable analogue or digital filter. The design of an analogue filter presents no unusual problems as the field is generally well known and details of design can be found, either in specialized literature or in several useful ‘recipes’ for the design of active filters’ The increasing use of digital computers in biomedical research makes a digital filter the obvious tool when the signal to be filtered consists of digitally sampled data which are to be processed by a computer anyway. The advantages of stability and the ease with which the frequency characteristics of a digital filter can be changed, make the digital filter even more attractive. The literature on digital filters is reasonably extensive and comprehensive yet digital filters do not seem as popular in practice as might be expected. Perhaps the potential difficulties involved in converting a few equations into a computer program are considered too great and a more familiar analogue filter is preferred.
DIGITAL
written in standard Fortran, is suitable for any computer with a Fortran IV compiler.
FILTERS
A different type of digital filter is a convolution or regressive filter. In this filter, each output value yn is calculated from a finite number of elements of the input sequence x, by an equation of the form: P Yn
=
ak h-k
c k=l
where ak is a set of suitable
weighting
coefficients.
This type of filter has the advantage of pure linear phase characteristics but gives a rather poor control over the frequency response. A different type of digital filter is the autoregressive or recursive filter. In the algorithm of this filter, not only the input samples x, , but also the previous output values y n_ m , are included in the calculation of the next output value: P Yn
=
c k=l
t UkXn-k
+
1 m=l
bm yn -m
The main distinction between the last two types of filters is that the regressive technique can operate with linear phase characteristics but needs a large number of terms to realise a sharp cut-off, whereas the autoregressive filter uses substantially fewer terms but is generally accompanied by some phase distortion due to the z-plane poles and zeros placed inside the unit circle.
The most obvious digital filter is one operating in the frequency domain by use of the discrete Fourier transformation’ and is usually implemented through an FFT algorithm. This allows a very tight control of the frequency response although this control is often difficult to realise in practice. Furthermore, the most frequently used algorithms for FFT require sample size lengths equal to a power of two, a requirement which becomes more difficult to arrange with longer samples.
It is possible to design a recursive filter which has its poles and zeros on the unit circle in the z-plane, and thus introduces no phase distortion3. This gain is, however, offset by a poor frequency response. The advantage of this type of filter is its on-line potential where it could be used up to a sampling frequency of several kI-Iz.
School of Biology, University of East Anglia, Norwich, NR4 7TJ, UK. * School of ClinicalMedicine, University of Cambridge, Cambridge, UK.
Verburg and Strackee4 gave a filtering procedure based on an autoregressive technique, giving a sharp
0 1982 Buttenvorth & Co. (Publishers) 0141-5426/82/030252-03 $03.00
252
Ltd
J. Biomed. Eng. 1982, Vol. 4, July
Progv-amfor digitalfilter: P.B. Pynscnt and R. Hanka ~.*****....*.**...****.~..~.*.*~~..***.~.***..*.*~~...*“*..~
C C
RECURSIVE MODE = -1 MODE = +l
C
FILTER : LOW PASS : HIGH PASS
C C C C**~~~~~**......*.*~~~~*~~~...**..**~.~~~***..****..~*~~***C SUBROUTINE FILTER (FREP. T. X. I. H@DE. DINENSION X(11 LOGICAL ITRAP C-----CHECK FCR NYQUIST FREQUENCY ETC. IERROR=O P1.1.0/(2.0*T) IF ((FRE~.LE.P~).oR.(I.GE.~)) GOTO loo IERROR=l RETURN 100 ITRAP=.FALSE. Y1=2.O*FREP'3.141592? C-----COMPENSATE FREQUENCY FOR WARPING W1:(2.0/T)*ATAN(Wl'T/2.0) Y3.U1**3.0 P1=2.O/T P3.Pl9.3.0 F=2.0*Pl'Wl YT.(P3+W3I+F*(Pl+WlI DlJW1=(3.O'(W3-P3)+F*o)/YT DUN2=(3.O'(W3+P3)-F'(Wl+Pl~)/YT DUM3:((W3-P3)+F*(Pl-Ul))/YT DUWII.W3/YT IF (MODE.CT.01 DUHU.P3/YT DUW5.DUM4*3.0 101 ITRAP=.NOT.ITRAP C-----FILTER DATA RECURSIVELY DO 102 N- 1.3 102 X(N)=O.O YH3=0.0 YM2tO.O YMlzO.0 IF (NODE) 103.103.105 C-----LOU PASS SECTION CONTINUE 103 DO 104 N-4.1 HX.DUN4.(X(N)+X(N-3))+DUN5'(X(N-2)+X(N-O) Yl.YNl'DUMl Y2=YW2*DUM2 Y3=YM3*DUM3 X(N-3).YM3 YX3=YM2 YM2=YHl YNl=HX-(Yl+Y2+Y3) 104 CONTINUE COT0 10'7 C-----HIGH PASS SECTION CONTINUE 105 DO 106 N=4,1 HX=DUM4*(X(N)-X(N-3))+DUMS*(X(N-2)-X(N-I)) Yl=YMl'DUMl Y2+YM2*DUW2 Y3=YH3*DUM3 X(N-3I.YH3 YR?=YMZ YHZ.YMl YHl=HX-(Yl+YZ+Y3) CONTINUE 106 C-----EMPTY OUT LAST POINTS CONTINUE 107 X(1-2).YM3 X(I-l).YH2 X(I)=YHl C-----REVERSE DATA FOR PHASE CORRECTION
IERROR)
II-I J:I/2 DO 108 N-l.J YMl.X(N) X(N)=X(II) X(II)=YMl 11=11-l 100
CONTINUE IF (ITRAP) COT0 101 RETURN END
Figure 1 Listing of the FORTRAN and high pass filtering.
subroutine for low
frequency response, which does not introduce a phase shift. Their idea was to use the same filtering procedure twice; once in the forward time order while storing the resulting sequence yi, and then once more, this time backwards operating on the already once filtered sequence yi in the reverse direction. FILTERING
SUBROUTINE
Verburg’s approach which seemed most advantageous for the design of a phaseless filter with good
C C C C
frequency response was used to produce the subroutine for digital filtering of sampled signals in this paper. The listing of the subroutine is given in F$.w 1. The subroutine, written in standard Fortran IV, operates on a real array (X) which can be of any length (I). The filter is a 3rd order Butterworth filter which is used either in a high or low pass mode determined by the parameter MODE (1 = high-pass; -1 = low-pass). The desired cut-off frequency in Hz (FREQ) and the data sampling interval in seconds (T) must also be passed to the subroutine. The cut-off frequency of the filter is constrained in both directions. On one hand it cannot exceed the Nyquist frequency (half of the sampling frequency) and on the other hand it is related to the sampling frequency and the precision of the computer, which is of crucial importance as the algorithm relies on accurate evaluation of very small differences between large numbers. (The single precision of mainframes and usual mini-computers is sufficient but a care should be taken when the subroutine is used on microcomputers which tend to differ in the number of significant digits used in arithmetic operations). The subroutine first checks that the parameters passed to it are in the correct ranges and then compensates for frequency warping introduced by the bilinear transformation’ as the cut-off frequency approaches the Nyquist frequency. The routine then calculates the individual coefficients of the recursive formula before filtering the array X according to the type of filter selected by MODE. The data array X, is then reversed and the control either returns to the calling program or the routine continues to filter the reverse data, depending on whether the pass is the first or the second one. The subroutine is designed to include only low and high pass filters. Figure 2a shows the amplitude response of this filtering subroutine when used as a low-pass filter with a cut-off frequency of 2 Hz, operating on a signal with sampling frequency of 100 Hz, both after a forward pass only and when used in both directions. The second pass through the filter resulted in general doubling of attenuation which also represents an effective shift of the 3 dB cut-off frequency. This may have to be taken into account and the desired frequency of the filter adjusted accordingly. Most striking is the effect of the forward and reverse operation on the phase of the output signal (Figure 2b). Here it is perhaps worth emphasizing that although these responses correspond very closely with the theoretical expectations, the responses in Figure 2 and 3 were obtained experimentally. The possibility of the subroutine being used for the band-pass and the band-stop filters was also considered but during extensive tests over a number of installations it was found that the precision of several commonly used small and even large computers was not sufficient, even when using dual precision variables to guarantee a proper operation.
J. Biomed. Eng. 1982, Vol. 4, July 253
Program for digital filter: P.B. Pynsent and R. Hanka
80
40 t2 0
2 !z r!
-40
Q
-80
2 -120
-160
1
Frequency
Frequency
10
1
0.1
( Hz 1
( Hz 1
Figure 2 An example of a low pass filter (FREQ = 2.0, T = 0.01). Output using the subroutine as written; result when processing was carried out in one direction only. (a) the amplitude response; )b) the phase response.
----,
The low and high pass filters can, however, be combined to represent either a band-pass or a band-stop processing of the data if required. The frequency response of a band-pass filter with fr = 0.75 Hz and f2 = 2 Hz obtained in this way is shown in Figure 3. The subroutine operates on an array which holds all data and it is assumed that sufficient memory is available. Even relatively small computers now seem to have memories of sizes which only a few years ago would have been difficult to imagine. Should it, however, prove impossible to store the whole signal in the memory, a trivial modification of the subroutine will make it possible to operate on data stored on a suitable sequential device such as a disc. Data could then be read from the disc in blocks (sectors) and stored again for the subsequent reverse processing. -50
ACKNOWLEDGEMENT This work was partially Foundation.
0.1
supported
by the Nuffield
REFERENCES 1
2
3
254
Derrett, C.J., Tavner, P J., Active filters: some simple design methods for biological workers, Med. Biol. Eng. 1975, 13, 883. Yoganathan, A.P., Gupta, R., Corcoran, W.H., Fast Fourier transform in the analysis of biomedical data, Med. Biol. Eng. 1976, 14, 239. Lynn, P.A., Outline digital filters for biological signals: some fast designs for a small computer, Med. Biol. Eng, 1977, 15, 534.
J. Biomed. Eng. 1982, Vol. 4, July
1
Frequency
10 (Hz)
Figure 3 Amplitude response of a band pass filter obtained by combining a low pass filter (FREQ = 2.0, T = 0.01, MODE = -1) followed by a high pass filter (FREQ = 0.75, T = 0.01, MODE = 1).
4
5
Verburg, J., Strackee, J., Phaseless recursive filtering applied to chest wall displacements and velocities using accelerometers, Med. Biof. Eng. 1974, 12, 483. Golden, R.M., Kaiser, J.F., Design of wideband sampled data filters, Bell Syst. Tech. J. 1964, 43, 1533.
and R. Hanka*
ABSTRACT A brief comparison of several types of digital filters is made and a simple subroutine implementing a recursive digital phaseless low and high pass filter is described. The subroutine, KeyworL:
computer programs, digital filter, phaseless filtering
INTRODUCTION In biomedical engineering it is often found that signals are contaminated by the presence of noise. If the frequency spectra of the signal and noise are separable then it is possible to supress the noise with a suitable analogue or digital filter. The design of an analogue filter presents no unusual problems as the field is generally well known and details of design can be found, either in specialized literature or in several useful ‘recipes’ for the design of active filters’ The increasing use of digital computers in biomedical research makes a digital filter the obvious tool when the signal to be filtered consists of digitally sampled data which are to be processed by a computer anyway. The advantages of stability and the ease with which the frequency characteristics of a digital filter can be changed, make the digital filter even more attractive. The literature on digital filters is reasonably extensive and comprehensive yet digital filters do not seem as popular in practice as might be expected. Perhaps the potential difficulties involved in converting a few equations into a computer program are considered too great and a more familiar analogue filter is preferred.
DIGITAL
written in standard Fortran, is suitable for any computer with a Fortran IV compiler.
FILTERS
A different type of digital filter is a convolution or regressive filter. In this filter, each output value yn is calculated from a finite number of elements of the input sequence x, by an equation of the form: P Yn
=
ak h-k
c k=l
where ak is a set of suitable
weighting
coefficients.
This type of filter has the advantage of pure linear phase characteristics but gives a rather poor control over the frequency response. A different type of digital filter is the autoregressive or recursive filter. In the algorithm of this filter, not only the input samples x, , but also the previous output values y n_ m , are included in the calculation of the next output value: P Yn
=
c k=l
t UkXn-k
+
1 m=l
bm yn -m
The main distinction between the last two types of filters is that the regressive technique can operate with linear phase characteristics but needs a large number of terms to realise a sharp cut-off, whereas the autoregressive filter uses substantially fewer terms but is generally accompanied by some phase distortion due to the z-plane poles and zeros placed inside the unit circle.
The most obvious digital filter is one operating in the frequency domain by use of the discrete Fourier transformation’ and is usually implemented through an FFT algorithm. This allows a very tight control of the frequency response although this control is often difficult to realise in practice. Furthermore, the most frequently used algorithms for FFT require sample size lengths equal to a power of two, a requirement which becomes more difficult to arrange with longer samples.
It is possible to design a recursive filter which has its poles and zeros on the unit circle in the z-plane, and thus introduces no phase distortion3. This gain is, however, offset by a poor frequency response. The advantage of this type of filter is its on-line potential where it could be used up to a sampling frequency of several kI-Iz.
School of Biology, University of East Anglia, Norwich, NR4 7TJ, UK. * School of ClinicalMedicine, University of Cambridge, Cambridge, UK.
Verburg and Strackee4 gave a filtering procedure based on an autoregressive technique, giving a sharp
0 1982 Buttenvorth & Co. (Publishers) 0141-5426/82/030252-03 $03.00
252
Ltd
J. Biomed. Eng. 1982, Vol. 4, July
Progv-amfor digitalfilter: P.B. Pynscnt and R. Hanka ~.*****....*.**...****.~..~.*.*~~..***.~.***..*.*~~...*“*..~
C C
RECURSIVE MODE = -1 MODE = +l
C
FILTER : LOW PASS : HIGH PASS
C C C C**~~~~~**......*.*~~~~*~~~...**..**~.~~~***..****..~*~~***C SUBROUTINE FILTER (FREP. T. X. I. H@DE. DINENSION X(11 LOGICAL ITRAP C-----CHECK FCR NYQUIST FREQUENCY ETC. IERROR=O P1.1.0/(2.0*T) IF ((FRE~.LE.P~).oR.(I.GE.~)) GOTO loo IERROR=l RETURN 100 ITRAP=.FALSE. Y1=2.O*FREP'3.141592? C-----COMPENSATE FREQUENCY FOR WARPING W1:(2.0/T)*ATAN(Wl'T/2.0) Y3.U1**3.0 P1=2.O/T P3.Pl9.3.0 F=2.0*Pl'Wl YT.(P3+W3I+F*(Pl+WlI DlJW1=(3.O'(W3-P3)+F*o)/YT DUN2=(3.O'(W3+P3)-F'(Wl+Pl~)/YT DUM3:((W3-P3)+F*(Pl-Ul))/YT DUWII.W3/YT IF (MODE.CT.01 DUHU.P3/YT DUW5.DUM4*3.0 101 ITRAP=.NOT.ITRAP C-----FILTER DATA RECURSIVELY DO 102 N- 1.3 102 X(N)=O.O YH3=0.0 YM2tO.O YMlzO.0 IF (NODE) 103.103.105 C-----LOU PASS SECTION CONTINUE 103 DO 104 N-4.1 HX.DUN4.(X(N)+X(N-3))+DUN5'(X(N-2)+X(N-O) Yl.YNl'DUMl Y2=YW2*DUM2 Y3=YM3*DUM3 X(N-3).YM3 YX3=YM2 YM2=YHl YNl=HX-(Yl+Y2+Y3) 104 CONTINUE COT0 10'7 C-----HIGH PASS SECTION CONTINUE 105 DO 106 N=4,1 HX=DUM4*(X(N)-X(N-3))+DUMS*(X(N-2)-X(N-I)) Yl=YMl'DUMl Y2+YM2*DUW2 Y3=YH3*DUM3 X(N-3I.YH3 YR?=YMZ YHZ.YMl YHl=HX-(Yl+YZ+Y3) CONTINUE 106 C-----EMPTY OUT LAST POINTS CONTINUE 107 X(1-2).YM3 X(I-l).YH2 X(I)=YHl C-----REVERSE DATA FOR PHASE CORRECTION
IERROR)
II-I J:I/2 DO 108 N-l.J YMl.X(N) X(N)=X(II) X(II)=YMl 11=11-l 100
CONTINUE IF (ITRAP) COT0 101 RETURN END
Figure 1 Listing of the FORTRAN and high pass filtering.
subroutine for low
frequency response, which does not introduce a phase shift. Their idea was to use the same filtering procedure twice; once in the forward time order while storing the resulting sequence yi, and then once more, this time backwards operating on the already once filtered sequence yi in the reverse direction. FILTERING
SUBROUTINE
Verburg’s approach which seemed most advantageous for the design of a phaseless filter with good
C C C C
frequency response was used to produce the subroutine for digital filtering of sampled signals in this paper. The listing of the subroutine is given in F$.w 1. The subroutine, written in standard Fortran IV, operates on a real array (X) which can be of any length (I). The filter is a 3rd order Butterworth filter which is used either in a high or low pass mode determined by the parameter MODE (1 = high-pass; -1 = low-pass). The desired cut-off frequency in Hz (FREQ) and the data sampling interval in seconds (T) must also be passed to the subroutine. The cut-off frequency of the filter is constrained in both directions. On one hand it cannot exceed the Nyquist frequency (half of the sampling frequency) and on the other hand it is related to the sampling frequency and the precision of the computer, which is of crucial importance as the algorithm relies on accurate evaluation of very small differences between large numbers. (The single precision of mainframes and usual mini-computers is sufficient but a care should be taken when the subroutine is used on microcomputers which tend to differ in the number of significant digits used in arithmetic operations). The subroutine first checks that the parameters passed to it are in the correct ranges and then compensates for frequency warping introduced by the bilinear transformation’ as the cut-off frequency approaches the Nyquist frequency. The routine then calculates the individual coefficients of the recursive formula before filtering the array X according to the type of filter selected by MODE. The data array X, is then reversed and the control either returns to the calling program or the routine continues to filter the reverse data, depending on whether the pass is the first or the second one. The subroutine is designed to include only low and high pass filters. Figure 2a shows the amplitude response of this filtering subroutine when used as a low-pass filter with a cut-off frequency of 2 Hz, operating on a signal with sampling frequency of 100 Hz, both after a forward pass only and when used in both directions. The second pass through the filter resulted in general doubling of attenuation which also represents an effective shift of the 3 dB cut-off frequency. This may have to be taken into account and the desired frequency of the filter adjusted accordingly. Most striking is the effect of the forward and reverse operation on the phase of the output signal (Figure 2b). Here it is perhaps worth emphasizing that although these responses correspond very closely with the theoretical expectations, the responses in Figure 2 and 3 were obtained experimentally. The possibility of the subroutine being used for the band-pass and the band-stop filters was also considered but during extensive tests over a number of installations it was found that the precision of several commonly used small and even large computers was not sufficient, even when using dual precision variables to guarantee a proper operation.
J. Biomed. Eng. 1982, Vol. 4, July 253
Program for digital filter: P.B. Pynsent and R. Hanka
80
40 t2 0
2 !z r!
-40
Q
-80
2 -120
-160
1
Frequency
Frequency
10
1
0.1
( Hz 1
( Hz 1
Figure 2 An example of a low pass filter (FREQ = 2.0, T = 0.01). Output using the subroutine as written; result when processing was carried out in one direction only. (a) the amplitude response; )b) the phase response.
----,
The low and high pass filters can, however, be combined to represent either a band-pass or a band-stop processing of the data if required. The frequency response of a band-pass filter with fr = 0.75 Hz and f2 = 2 Hz obtained in this way is shown in Figure 3. The subroutine operates on an array which holds all data and it is assumed that sufficient memory is available. Even relatively small computers now seem to have memories of sizes which only a few years ago would have been difficult to imagine. Should it, however, prove impossible to store the whole signal in the memory, a trivial modification of the subroutine will make it possible to operate on data stored on a suitable sequential device such as a disc. Data could then be read from the disc in blocks (sectors) and stored again for the subsequent reverse processing. -50
ACKNOWLEDGEMENT This work was partially Foundation.
0.1
supported
by the Nuffield
REFERENCES 1
2
3
254
Derrett, C.J., Tavner, P J., Active filters: some simple design methods for biological workers, Med. Biol. Eng. 1975, 13, 883. Yoganathan, A.P., Gupta, R., Corcoran, W.H., Fast Fourier transform in the analysis of biomedical data, Med. Biol. Eng. 1976, 14, 239. Lynn, P.A., Outline digital filters for biological signals: some fast designs for a small computer, Med. Biol. Eng, 1977, 15, 534.
J. Biomed. Eng. 1982, Vol. 4, July
1
Frequency
10 (Hz)
Figure 3 Amplitude response of a band pass filter obtained by combining a low pass filter (FREQ = 2.0, T = 0.01, MODE = -1) followed by a high pass filter (FREQ = 0.75, T = 0.01, MODE = 1).
4
5
Verburg, J., Strackee, J., Phaseless recursive filtering applied to chest wall displacements and velocities using accelerometers, Med. Biof. Eng. 1974, 12, 483. Golden, R.M., Kaiser, J.F., Design of wideband sampled data filters, Bell Syst. Tech. J. 1964, 43, 1533.