Such filters are used in various fields where it is necessary for signals which are limited in bandwidth and are sampled at twice the Nyquist frequency to be recon- structed at any time in respect to the sampling times. These interpolation filters are also referred to as fractional delay filters or variable phase delay filters.

Fields of application include, but are not limited to, synchronisation in digital modems, phase-controlled array antennae (for example"null-steering"antennae), asyn- chronous sample frequency conversion of previously sampled signals ("upsampling" and"downsampling"), speech coding, signal pattern recognition and temporal localisa- tion of patterns (for example speech recognition), periodical or non-periodical signal pattern generation in any phase desired, physical modelling of musical instruments (for example vibrating strings), but also the tracking of signals in digital navigation receiv- ers (GPS, LORAN-C, EUROFIX), the time-continuous tuning of numerically con- trolled oscillators (NCO) and adaptive averaging devices (in LORAN-C), the direct determination of signal characteristics such as zero crossings and maxima/minima in, for example, seismic measurement data, the change in the aspect ratio (3/4 <-> 9/16) in the case of video processing, simulating multipath phenomena in the case of antenna signals in mobile receivers, the rapid determination of (local) maxima and minima in sampled data (for example in GPS data processing), time-discrete Hilbert transformers for software radio receivers, signal pattern recognition for temporal and/or positional localisation in, for example, the medical sector (ECG (heart rhythm analysis), EEG, etc) and the synchronisation of many asynchronous digital data streams into a single data stream (for example the synchronisation of compass data, gyroscope data, level data, roll data, pitch data, (GPS) position data, speed data, acceleration data, and data

regarding the change in acceleration into one navigational solution for aeroplanes), or the synchronisation of various asynchronous process data streams into one data stream for digital control systems having a plurality of inputs in all forms of industry (i. e. a simplification of Kalman filtering).

Such filters are used as an alternative for known mathematical interpolation methods for determining the best curve through a particular set of equidistant points. If the original signal, after filtering by an ideal low pass filter has been sampled at twice the Nyquist frequency, the original signal can be reconstructed perfectly by convoluting the sample values with the impulse responses associated with the ideal low pass filter.

The impulse response of the ideal low pass filter is a time-continuous sinc function where the period of the side lobes is equal to the inverse of twice the Nyquist fre- quency. For the ideal case, even one convolution operation will already take an in- finitely long time for such a filter of infinite order M. In practice, however, the only possibility is to work with causal impulse responses and a limited order of the inter- polation filter, the so-called finite impulse response filters (FIR).

The frequency response of a fractional delay filter, implemented as a finite impulse response filter, is generally given by where hj represents the impulse response weights and M is the filter order.

Finite impulse response filters with fractional delay are disclosed, for example, by the article"Splitting the Unit Delay"by T. I. Laakso et al., in IEEE Signal Process- ing Magazine, Volume 13, No. 1, pp. 30-60, January 1996.

This publication discloses, for example, a design for a finite impulse response filter with fractional delay, where the infinite sinc impulse response is multiplied by a bell-shaped window function. The zero values outside the bell shape ensure that the ultimate impulse response will be constrained or finite. The computation of the filter coefficients for such a finite impulse response filter is complex (for each coefficient determination, a value of a computation-intensive bell-shaped window function must be multiplied by a value of a computation-intensive sinc function), and the error intro- duced by the filter is not minimise, moreover.

Another way, disclosed in this publication, of designing finite impulse response filters with fractional delay is the so-called Lagrange interpolation. Here, the filter coef- ficients hj are given by Here, too, many computation-intensive operations are required for each coeffi- cient, thus increasing the total operation time per interpolation of the filter.

Since the known filter designs require many, often complex computations for each filter coefficient of the finite impulse response filter, it is often not possible to achieve the accuracy, speed of adjustment to new values of the fractional delay d, bandwidth and/or maximum allowable total filter delay which are required for a par- ticular application.

It is an object of the present invention to provide a method for obtaining an in- terpolated value of a digital signal with a fractional delay d as a variable component of a delay of M/2 sampling times, with the aid of Mth order finite impulse response fil- tering, which requires the smallest possible number of operations per filter coefficient, so that highly accurate, rapidly adjustable filtering with fractional delay can take place at the highest possible bandwidth even for small values of the total filter delay M/2+d.

This object is achieved in a method of the type defined in the preamble, charac- terised in that the window function is a symmetric stepped window function and in that the method further comprises the following steps: -determining simplified filter coefficients (hj) wo as a function of the fractional delay d according to the formula and subsequently applying the finite impulse response filtering to the digital signal;

-multiplying the weighted result of the finite impulse response filtering by a common normalisation factor (F (d)) o, given by the formula:

which is a function of the simplified filter coefficients, in order to obtain the inter- polated signal value from the digital signal with fractional delay d.

This method has the advantage that each filter coefficient globally requires only three operations, namely an addition, a subtraction and a division. This allows the fabrication of finite impulse response filters which can be very rapidly adjusted to the fractional delay d and consequently have a small total filter operation time of the size of 4M operation times.

In one embodiment, the window weights of the symmetric stepped window function are chosen such that the resulting finite impulse response filter exhibits a fre- quency response having a ripple-shaped deviation of the normalised level (equiripple filter). In a further embodiment, the window weights of the symmetric stepped window function are chosen such that the resulting frequency response is maximally flat (Lagrange filter). This means that a choice can be made from a number of classes of finite impulse response filters by choosing the window weights without this having consequences in terms of the way in which the filter coefficients are calculated.

As a result of, in one embodiment of the method according to the invention, an additional common factor (l/4-d2) being added when the simplified filter coefficients (hj) Wo and the common normalisation factor (F (d)) wo are determined, the method can <BR> <BR> also be applied to a fractional delay d equal to i1/2. In the embodiment mentioned ear- lier this value was excluded to prevent a division by zero from determining the simpli- fied filter coefficients. To this end, in fact, two critical divisions are replaced by two multiplications, the total number of operations increases slightly.

In one embodiment of the method according to the present invention, two finite impulse response filtering operations are used to effect filtering which allows all fre- quencies to pass ("all pass filter"), the simplified filter coefficients being determined

solely from the window weights and the desired fractional delay d. The transfer func- tion of the all pass filter is given by Hence it can be deduced that both the multiplication by the common normalisa- tion factor (F (d)) wo and the addition of the factor (l/4-d2) can be dispensed with for this filter.

The above-described embodiments of the method according to the invention have the drawback that one of the operations in determining the simplified filter coeffi- cients and the common normalisation factor is a division. A division is an operation which, both when carried out by software and carried out by hardware, generally re- quires much more operation time than a multiplication operation.

In a further embodiment of the method according to the invention, the division in the course of determining the filter coefficients and the common normalisation factor is replaced by an approximation with the aid of one or more multiplications. Thus, the finite impulse response filtering with fractional delay d can be carried out more rapidly, compared with the previously described embodiment, with an acceptable approximation error.

In a further embodiment of the method according to the invention, the approximation for a main lobe of the impulse response and the approximation for side lobes of the impulse response take place separately in such a way that the resulting im- pulse response retains the exact values for the extreme values of the fractional delay (i. e. for d=0 and d=il/2).

In a first approximation, the finite impulse response filtering with a fractional delay d, where the simplified filter coefficients and the common normalisation factor are calculated with the aid of multiplications, can be described by a third-order poly- nomial (P0=3) of the fractional delay, the maximum approximation error in the impulse response being 4 10-3. The complexity of this embodiment can, based on the approxi- mation by a third order polynomial, be reduced in a further approximation to an em- bodiment using a second-order polynomial (P0=2) with a maximum approximation error remaining almost the same at 6 10-3. In an alternative approximation, the filtering

can be described by a fifth-order polynomial (PO=5), the maximum approximation error being 2-10-5. This, however, is at the cost of doubling the complexity of the filtering.

In a further embodiment of the method according to the invention use is made of the inverse of the finite impulse response filtering with fractional delay d as described above, to determine those values of the fractional delay at which the original signal has zero crossings and/or extremes. Not only the point in time and the corre- sponding amplitude, but also the conditions in which these phenomena occur can be determined accurately if the order of the polynomial is less than or equal to 3 (P0<3). In the case of a third-order polynomial approximation, zero crossings can be determined via a solution of a third-order equation and extreme values according to a simpler solution of a second-order equation. If the inverse of the second-order polynomial approximation of the fractional delay filter according to the present invention is used, zero crossings can be determined by solving a second-order equation, and extremes can be determined by solving a linear equation.

A second aspect of the present invention relates to a device for implementing the method according to any one of Claims 1 to 7 inclusive. To this end, the device is provided with delay elements, addition elements, multiplication elements and inverter elements for filtering a digital signal with a fractional delay d which is a fraction of the sampling period associated with the digital signal, predetermined window weights wk of a window function which is applied to the infinite impulse response of an ideal interpolator being implemented in multiplication elements, simplified filter coefficients (hj) Wo being implemented as a combination of addition elements, multiplication elements and inverter elements and a common normalisation factor (F (d)) wo being implemented as a multiplication ele- ment, in order to obtain the interpolated signal value from the digital signal with fractional delay d.

As this device is capable of hardware implementation, much more rapid opera- tion of the filter can be obtained, compared with a software implementation.

In a further embodiment of the device according to the present invention, the device is further provided with computing means in which the inverse of the finite im- pulse response filtering with fractional delay d is employed to determine the values of the fractional delay d at which the original signal has zero crossings and/or extremes.

A possible application of the method for filtering a digital signal in order to re- construct a value of a sampled signal at any time between two samples of the digital signal is a numerically controlled oscillator (NCO). Known NCOs (for example the type HSP45116 from Harris Corporation) comprise a K-bit phase accumulator and an electronic table containing 2K numbers having a word size of 16 bits. In the case of a 20-bit phase accumulator this leads to a table of 16 Mbit. With the aid of the phase accumulator, any frequency can be chosen in the range from 0 to fs/2 in steps of 2-K*fs.

A further aspect of the present invention relates to a numerically controlled oscillator for generating digitised band-limited periodic signals, for example sinusoidal signals, comprising a K-bit phase accumulator and an addressable table characterised in that the addressable table comprises twice the oversampling (=2*b) of groups of PO+1 polynomial coefficients which are addressed by the log (2*b) most significant bits of the K-bit phase accumulator, the addressed PO+1 polynomial coefficients being used for a POth order approximation of a filter with a fractional delay d, formed by PO multiplication elements and PO addi- tion elements, the fractional delay d being determined by the 2log (2*b) least significant bits of the K-bit phase accumulator. Since only 2*b groups of PO+I polynomial coeffi- cients need be stored in the table, a much smaller table is sufficient. In the case of, for example, a sine signal oversampled by a factor of 16 (2b=32) and a third-order poly- nomial approximation (P0=3) of a 5th order (M=5) fractional delay filter, only 128 numbers need be stored in the table. Simulations have shown that for a second-order (P0=2) polynomial approximation about 1000, and for a first-order (PO=1) polynomial approximation about 8000 numbers should be stored for comparable quality.

The present invention will now be explained in more detail with reference to examples and the appended drawings, in which Figure 1 shows a finite impulse response of a seventh-order filter, together with a window function which is applied to the infinite impulse response of the ideal low pass filter; Figure 2 shows a general schematic for the fractional delay filter, implemented using a finite impulse response filter; Figure 3 shows an example of a schematic of a fifth-order fractional delay filter, in which the range of the fractional delay d has been extended by the value ='4; Figure 4 shows an example of a schematic of a third-order fractional delay filter

which passes all frequencies to a same degree; Figure 5 gives the general schematic for a ninth order filter with 10 filter coeffi- cients according to an embodiment of the present invention; Figure 6 shows the schematic for a design of an integrated circuit for the first approximation (third-order polynomial) of the filter shown in Figure 5; Figure 7 shows the schematic for a design of an integrated circuit for the second approximation (fifth-order polynomial) of the filter shown in Figure 5; Figure 8 shows the schematic for a design of an integrated circuit of the filter shown in Figure 5 according to a second-order polynomial approximation, derived from the embodiment using a third-order polynomial approximation; Figure 9 shows the schematic for a design of an integrated circuit for determin- ing times of zero crossings and extremes and the corresponding amplitudes with the aid of the second-order polynomial approximation showing Figure 8; Figure 10 shows the schematic for a numerically controlled oscillator according to a further aspect of the present invention.

The frequency response of a fractional delay filter, which has been implemented as a finite impulse response, can generally be described as Figure 1 shows a finite impulse response of a seventh-order filter (a fragment of an infinite sinc function), together with the limiting window function which, in accor- dance with the method according to the invention, is applied to the infinite impulse response associated with the ideal interpolation filter. The symmetric window weights <BR> <BR> are given by W... W4, the symmetry being round the centre of gravity of the impulse response. The filter coefficients ho... h7 are indicated by the vertical lines at times 0, T, ..., 7T, where T is the sampling period. The total filter delay time is indicated in Figure 1 by the phase delay D and the fractional delay by the letter d.

To filter a digital signal with the aid of the method according to the invention, the symmetric stepped window function is first applied to the infinite impulse response of the ideal low pass filter to produce the Mth-order interpolation filter. Because of the symmetry of the window function it is possible to use the parameters i and k, thus

allowing the filter coefficients to be divided into two equal parts j= {i; M-i}. For the example in Figure 1, the values of t (s), j, k and i are shown below the graph. Conse- quently, the filter coefficients hj are derived by calculating the simplified filter coef- ficients (hj) wo according to the equation

whereupon the total result of the digital signals filtered therewith are multiplied by the common normalisation factor (F (d)) wo, which is calculated using the equation

The window weights Wk are determined as a function of the desired design type of the filter. For an equiripple filter, the window weights are given by

where the coefficients (hk) mr are determined via existing design software for equiripple filters.

For a maximally flat filter or Lagrange filter, the window weights are given by

Shown below as an example are the window weights which were calculated with the aid of a computer program for a finite impulse response filter where the phase delay ripple remains below about 0.1% and where the frequency response ripple remains below a value of 0.7%, based on the normalised frequency response. M is the filter order, and N=M+1 is the number of filter coefficients. N M Wl W2 W3 W4 W5 1000021 4 3 1 0.4 0 0 0 6 5 1 0.6 0.18 0 0 8 7 1 0.74 0.38 0.1 0 10 9 1 0. 84 0.57 0.57 0.11

The method according to the invention is also suitable for implementing a Lagrange filter. A Lagrange filter has a maximally flat filter response, but a smaller frequency range compared with the above-mentioned equiripple filter. The window weights in this case are given by the following table: N MW) W2 W3 W4 W5 2 1 1 0 0 0 0 4 3 1 1/3 0 0 0 6 5 1 1/2 1/10 0 0 8 7 1 3/5 1/5 1/35 0 10 9 1 2/3 2/7 1/14 1/126

Figure 2 shows a general schematic for the fractional delay filter, implemented by means of a finite impulse response filter. The schematic incorporates delay ele- <BR> <BR> ments, indicated by Z~t, multiplication elements, indicated by >, and a summation

element, indicated by +. The filter coefficients (hj) Wo and the common normalisation factor (F (d)) wo are calculated in the above-described manner. Once the filter order and the window weights have been chosen, the filter coefficients depend solely on the frac- tional delay d.

Figure 3 shows a schematic of an example of a fifth-order functional delay filter <BR> <BR> <BR> <BR> where the range of the fractional delay d has been extended by the value il ='/2. This is possible by introducing an additional common factor (l/4-d2) into the computation'of the filter coefficients. The total number of operations is increased slightly, because the implementation in practice of twice this additional factor in total requires four more multiplications plus two more additions. It should be noted in this context, that at the same time two critical divisions (critical because of division by zero) have been replaced by two safe multiplications.

Figure 4 shows an example of a third-order fractional delay filter which allows all frequencies to pass. These filters are used for specific applications requiring a com- pletely flat frequency response. The filter coefficients are calculated in the same way as for the previously described filters, whilst the multiplication by the common normal- sation factor (F (d)) wo can be dispensed with, as can the above-mentioned provision to avoid a division by zero. This last-mentioned advantage is a result of the fact that, compared with FIR filters, the fractional delay in all pass filters is increased by a factor of two, thus allowing the range of d to be halved and the critical division to be reliably avoided. Consequently, the third-order filter shown as an example in Figure 4 requires only three delay elements, eight multiplication elements, five addition elements and one subtraction element.

According to a further embodiment of the method according to the invention, the divisions carried out to calculate the simplified filter coefficients and the common normalisation factor are replaced by approximations using multiplications. In a further embodiment of the method according to the invention, the approximations for the side lobes of the impulse response (k > 1) and the approximations for the main lobes (k = 1) of the impulse response are performed separately. Hereinafter, the filter coefficients for the approximation will be indicated by h.

A first approximation (indicated by subscript 3 because of the third-order poly- nomial approximation) of the filter coefficients of the side lobes of the impulse response is given by

where subscript i represents the rank of the coefficients from the"left-hand"term of the impulse response to which only the"upper"part of the +, i. e. a minus sign, applies and where subscript k represents the rank of the coefficients from the"right-hand"term of the impulse response to which only the"bottom"part of the +, i. e. a plus sign, applies and where also <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> Wk<BR> <BR> <BR> <BR> ' (k- (O) is a constant and the common normalisation factor FM (0) is given by A second approximation (indicated by subscript 5 because of the fifth-order polynomial approximation) of the side lobes of the impulse response can be given by adding to the first approximation of the filter coefficients a term (I-Bk*d2), the calcula- tion of the filter coefficients consequently being given by Approximations can also be derived for the main lobes, resulting in a first approximation of the filter coefficients given by

where only the"upper"part of the +, i. e. a minus sign, applies to the coefficients from the left-hand term of the impulse response and where only the"bottom"part of the +, i. e. a plus sign, applies to the coefficient from the right-hand term of the impulse response and where also CCI = Cl-4.

A second approximation is provided by the addition of the factor (1-B1*d2), resulting in the formula The correction factors Bk can be calculated almost exactly for a specific value of the fractional delay, namely d = 0.361. At this value of the fractional delay, the first approximation error reaches a maximum as the filter order M increases.

Again, two situations can be distinguished, viz. the correction factors Bk (k > 1) for the side lobes of the impulse response and the correction factor B, for main lobes of the impulse response. This results in where the normalisation factor FM (0.361) is found by substituting d = 0.361 into the equation

The second approximation provides a more accurate approximation at the expense of a number of additional multiplication operations.

Figure 5 gives the general schematic for a ninth-order filter with 10 filter coeffi- cients. The section surrounded by a dotted line relates to the correction factors used for the second approximation. This form is suitable to be implemented, for example, in a software program for a digital signal processor.

A topology for both a software implementation and for the design of an inte- grated circuit can be found by rearranging the equations for the filter coefficients, so that a minimum number of multiplication operations using the fractional delay d are required.

The equations for the filter coefficients in the first approximation (Bk = 0) for the filter from Figure 5 are then given by a third-order polynomial in d: The schematic for a design of an integrated circuit which comprises delay ele- ments, multiplication elements, addition elements and inversion elements, for the first approximation (third-order polynomial), is shown in Figure 6. The maximum approxi- mation error in the impulse response of this finite impulse response filter with frac- tional delay is about 4 10-3.

The equations for the filter coefficients in the second approximation (Bk W 0) for the filter from Figure 5 are then given by a fifth-order polynomial in d:

The schematic for a design of an integrated circuit for the second approximation (fifth-order polynomial) is shown in Figure 7. The maximum approximation error in the impulse response for this finite impulse response filter with fractional delay is about 2 10-5.

From the formulae for the filter coefficients using the third-order polynomial approximation, filter coefficients can be derived for a second-order polynomial approximation. The filter coefficients for the side lobes of the filter using a second- order polynomial approximation (where subscript 2 is used to indicate the second-order polynomial approximation) are given by the formulae for 0 d'/2, and for-'/2 <_ d 5 0. For the main lobes, the filter coefficients are given by for 0 < d _<'/z, and

for-'/2 d 0. The maximum approximation error in the impulse response is less than 6-10-3 (in the equiripple embodiment), which is virtually no greater than in the case of the third-order polynomial approximation.

Figure 8 shows a schematic for a design of an integrated circuit of the filter shown in Figure 5 according to a second-order polynomial approximation, derived from the embodiment using a third-order polynomial approximation. The equations for the filter coefficients in this approximation for the filter from Figure 8 are then given by a second-order polynomial in d: for positive d: and for negative d: The schematic shown in Figures 6,7 and 8 can serve as a basis both for soft- ware implementations of the filters and for the design of integrated circuits.

In a further embodiment of the present invention, the inverse of the finite response filtering with fractional delay d is used for determining the values of the frac- tional delay d, at which the original signal has zero crossings and/or extremes. In a

fractional delay filter, the fractional delay parameter d (within the total phase delay) determines to what extent input samples are shifted in time by interpolation. For each clock pulse, a polynomial function in d is calculated. If the design of the fractional delay filter is such that the polynomial coefficients are provided directly, this filter is suitable to serve as an inverse fractional delay filter. The polynomial function and its <BR> <BR> <BR> <BR> derivative then directly enable the positions in time of zero crossings (dzero) and/or of extremes (deXtr) to be calculated. If the delay found at an extreme (deXcr) is introduced'as the input of a fractional delay filter, this directly gives the corresponding amplitude value as the output signal of the filter. Zero crossings can only occur between two suc- cessive sample values of opposite signs. To find the zero crossings it is first necessary, therefore, to determine where in the series of sample values this condition is met.

The exemplary schematics of Figures. 6,7 and 8 have been elaborated in such a <BR> <BR> <BR> <BR> way that the separate polynomial coefficients ao, a,, a2... are calculated directly. From these coefficients, the values of the fractional delay at zero crossings can be determined directly and without iterations via the inverse fractional delay filter concept, and after differentiation, extremes of the signal such as maxima and minima can be determined.

The calculations of these values of the fractional delay are the simpler, the lower the order of the polynomial approximation. The above-discussed second-order polynomial approximation is consequently the most suitable for use in the inverse filter with frac- tional delay. Zero crossings can therefore be found by solving the equation f2 (d) =ao+a, *d+a*d'=0 The extremes can be found by solving the equation 2 (d) =a, +2*a, *d=0 where the type of extreme can be found by a calculation of the second derivative if <BR> <BR> <BR> <BR> f2" (d) =2*a2, the extreme relating to a minimum if f2" (d) >0, a maximum if f2" (d) <0<BR> <BR> <BR> <BR> <BR> <BR> and a point of inflection if f2"(d) =0.

Such determinations of zero crossings can also be carried out with the aid of the third-order polynomial approximation, as discussed with reference to Figure 6. Figure 9 shows the schematic for a design of an integrated circuit for determining the times of

zero crossings and extremes and the corresponding amplitudes with the aid of the second-order polynomial approximation shown in Figure 8. This embodiment is provided with computation means 10 which perform the above-mentioned calculations of the zero crossings and extremes. If the fractional delay at which an extreme occurs (dextr) is directly fed back to the filter, the amplitude value at the extreme is supplied at the output of the filter.

Hereinafter a number of practical applications will be discussed for using the method or device according to the invention.

Hilbert transformers in the ideal case are based on infinite sinc functions. This opens the possibility of applying the concept of filtering with fractional delay according to the invention, resulting in a finite impulse response by using the windowing tech- nique with predetermined window weights. The ideal Hilbert transformer consists of two filters having a flat frequency response with output signals phase-shifted by pre- cisely 90 degrees for common real input signals. The two output signals form a single complex output signal. The windowing technique according to the invention cannot be applied directly to the two impulse responses of the ideal Hilbert transformer which are currently found in the literature, as there is an odd number of samples of the sinc im- pulse response within the even number of window weights which has to be used. This can be solved by causing the impulse responses of the two filters in the Hilbert trans- former to be mirror-antisymmetric and of the same odd order. The filter coefficients of <BR> <BR> <BR> <BR> <BR> the two filters are then found by substituting the value d ='/4 (resulting in the real im- pulse response) and d = _I/4 (resulting in the imaginary impulse response) into the inter- polation formulae according to the invention. Real bandwidth-limited signals (for example coming from an antenna) are thereby directly converted into complex band- width-limited signals, half-band filtering being used in the process. Both time errors (or phase errors) and the spurious signals introduced thereby, which occur in conventional Hilbert transformers, are thus eliminated.

A further application can be found for numerical controlled oscillators (NCO) like, for example, those used in software programs. Known NCOs, such as, for example, the type HSP45116 from Harris Corporation, consist of two parts, namely a digital K-bit phase accumulator linked to an electronic table which comprises stored sample values of a sine wave signal. The table contains 220 (about one million) numbers (for K=20) having a word size of 16 bits, (i. e. a total of 16 Mbit of data). By using the

phase accumulator it is possible to select any frequency value in the range from 0 to fs/2 at steps of 2-20 * fs.

The embodiments of the fractional delay filters, as discussed in the present application, all have the characteristic that they exhibit very low distortion both in phase delay and in the level of the frequency response for signal frequencies close to zero. As a result, these are eminently suitable for constructing all sinusoidal signal values by interpolation of a relatively small number of sample values of a constant sine wave which, for example, is oversampled by a factor of 16 (b=16). Starting from a set of at least 2*b=32 stored sample values of a sine wave (amply sufficient for charac- terising the entire period), it is possible to construct a numerical controlled oscillator which is able to generate a very accurate digital sine wave of any desired frequency.

Computer simulations have shown that the best results are achieved by using a third- order polynomial approximation (P0=3) and a filter order of M = 5. In this case, only 32 sine signal samples need be stored for a better quality, compared with the known NCO (error of 10-7 instead of error of 1.5* 10-5). If use is made of filters employing a lower-order polynomial approximation, more sine signal samples have to stored for comparable quality. Simulations have shown that about 250 sine signal samples should be stored for a second-order (P0=2) approximation and about 2000 with a first-order (PO=1) approximation.

If the approach is based on an interpolator with fractional delay using a third- order polynomial approximation, which is fed a periodic sine wave which, for example, is oversampled by a factor of 16 (b=16), a set of four polynomial coefficients {ao, al, a2, a3} is calculated at each sampling time. Owing to the periodicity of the digitised sine wave, not only the 32 samples of the sine wave but also the 32 related groups of four polynomial coefficients will be repeated periodically. This implies that instead of storing 32 sine wave samples it is better to store the 32 related groups of polynomial coefficients. As a result, a large number of operations of the filtering with fractional delay can be dispensed with. Only the calculation of the polynomial value at each sampling time remains. A possible set up for a numerical controlled oscillator is shown in Figure 10. The NCO comprises a K-bit phase accumulator 20 and an addressable table 22, for example in the form of a read-only memory (ROM). The K-bit phase accumulator 20 converts the desired frequency F into an addressing word and a delay word. The way in which this is done is assumed to be known to those skilled in the art

and will not be explained here in more detail. The addressable table 22 comprises 2*b sets of PO+1 polynomial coefficients, which are addressed by the addressing word which comprises the 2log (2*b) most significant bits of the K-bit phase accumulator 20.

The addressed PO+1 polynomial coefficients (in the examples shown ao, ai, a2, a3) are then used for a POth order approximation of a filter with a fractional delay d, formed by PO multiplication elements and PO addition elements. The fractional delay d is determined by the delay word which is formed by the K-21Og (2*b) least significant bits of the K-bit accumulator 20. Clearly, in this embodiment the ROM 22 can be much smaller than in the known NCO. In the known NCO, about one million numbers need to be stored. If use is made of a third-order (P0=3) polynomial approximation of a fifth order (M=5) filter with fractional delay, only 128 polynomial coefficients need to be stored to achieve comparable quality. The inverter 24 in the schematic shown in Figure 10 changes the binary representation of the fractional delay d to a two-comple- ment representation. The multiplexer 26 shown in Figure 10 has been added to allow the choice between generating a sine or a cosine signal by feeding the two most signifi- cant bits of the addressing word directly to the ROM 22, or via a second inverter 28 and an EXOR-gate 30. The frequency of the sine/cosine is defined, in the example shown, by the relationship: frequency = (F/2K) *f, where 0 F 2K- In the case of a signal oversampled by a factor of 16, a 5th order filter with fractional delay and a third-order polynomial approximation (b=16, M=5, P0=3 implying four polynomial coefficients), the polynomial coefficients are then given by: xn-5)-3C2(xn-1+xn-4)+C1(xn-2+xn-3)}a0(n)=#{5C3(xn+ a, (n) =-4au +2 (x+x) <BR> <BR> <BR> <BR> <BR> xn-5)+C2(xn-1+xn-4)+C1(-xn-2+xn-3)}a1(n)=#{C3(-xn+ <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> <BR> a3 (n) =-4a, +4 (x,,-x,,) where xn (not delayed) to Xn-5 (delayed by 5 clock pulses) inclusive form a set of six successive sample values from a set of 2*b sine wave sample values which are calcu- lated according to

x, = sin((n + ###)*#) n = {0,1, 2,..., (2 * b-1)} It will be obvious that instead of periodic digitised and oversampled sine wave shapes it is also possible to use other arbitrarily chosen periodic wave shapes which are limited in bandwidth and have been oversampled.

Many more fields than the two for which examples have been elaborated here- inabove lend themselves to the application of the method for filtering with a fractional delay according to the present invention, including applications requiring time-continu- ous adjustable delay lines, synchronisation of receivers by using inverse filtering with fractional delay, etc. All these elaborations of the general concept of filtering with frac- tional delay are deemed to be within the scope of the appended claims.