Fractional delay digital filter

This application claims priority under 35 U.S.C. §365(a) from International Application No. PCT/FR98/00615, filed Mar. 26, 1998, and published under PCT Article 21(2), on Oct. 8, 1998, in French, which is hereby incorporated by reference.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention concerns digital signal processing devices which perform fractional delay operations in a linear processing context.

The theory of such devices is well known. Conventionally, and as shown in

FIG. 1

, these devices comprise on the input side a row of finite impulse response interpolation filters each of which performs a fractional delay operation specific to each input signal.

Each finite impulse response interpolation filter is an approximation of an ideal fractional delay filter. The coefficients of the finite impulse response interpolation filters are calculated using either conventional methods which minimize the mean square error or the Lagrange method.

The delayed signals are then each subjected to specific linear processing, after which they are added. There are also similar devices using an arrangement which is the opposite of that just mentioned. In this case, the same input signal is duplicated to provide a plurality of signals identical to it, which signals are then subject to specific linear processing and a specific fractional delay operation. Such devices are used in voice coders employing fractional pitch resolution, for example.

2. Description of Related Art

FIGS. 1

to

3

show a prior art digital device in which M input signals u

1

, u

2

, . . . , u

M

are delayed before they are subjected to linear processing by a filter T

1

through T

M

and added by an adder S. This type of device is used in conventional applications, for example antenna processing, in which case the delays are chosen as a function of the antenna pointing angle.

The required delay r

k

is rarely an integer multiple of the sampling period. The delay r

k

therefore comprises an integer number of sampling periods to which a non-zero fraction τ of a period is added. Delays which are a multiple of the sampling period are not considered hereinafter because they are particularly simple to implement. The delays to be obtained are considered to consist of this additional fraction of a sampling period.

The device shown in

FIGS. 1

to

3

is fed with M input signals u

k

where k varies in the range from 1 to M. As shown in

FIG. 1

, each of the M input signals u

k

is delayed by a fractional delay r

k

specific to it.

A finite number of fractional delays r

k

is therefore required. A highest common factor r

o

of the fractional delays r

k

can be found. In practice an interpolation filter is implemented for each multiple of r

o

in the range from zero and the value of the sampling period T

e

. As shown in

FIG. 2

, D interpolation filters h

i

, denoted R

0

, R

1

, . . . , R

D−1

are obtained in this way, producing delays of value

\frac{i}{D} \times T_{e},

where i varies in the range from 0 to (D−1).

Because the same number L of coefficients is chosen for each filter, D filters h

i

(n) are obtained, where i varies in the range from 0 to (D−1) and n varies in the range from 0 to (L−1).

As shown in

FIG. 3

, each input signal u

k

, where k varies in the range from 1 to M, is submitted to an interpolation filter h

τk

(n), denoted RF

k

, which is one of the D filters h

i

(n) denoted R

i

in

FIG. 2

, i being in the range from 0 to (D−1).

The equations and the criteria used to calculate the coefficients h

i(n)

of the D interpolations R

0

to R

D−1

will now be summarized.

Let τ

k

denote the normalized fractional delay value, that is to say

τ_{k} = \frac{r_{k}}{T_{e}}

where r

k

is the required fractional delay for the signal and T

e

is the sampling period.

Let u

k

(n) be the input signal obtained by sampling an analog u

k

0

(t). Assuming that the sampling was done correctly, it is theoretically possible to obtain the delayed digital signal, defined as the samples of the analog signals delayed by u

k

(n), from u

k

0

(t−τ

k

T

e

). The skilled person knows that from the signal u

k

(n), all samples of which are known, the corresponding analog signal u

k

0

(t) can be found using the interpolation equation

u_{k}^{0} (t) = \sum_{m = - \infty}^{\infty} u_{k} (m) \frac{\sin (π (\frac{1}{T} - m))}{(\frac{1}{T_{e}} - m)} = \sum_{m = - \infty}^{\infty} u_{k} (m) sinc (\frac{1}{T_{e}} - m)

in which case the delayed analog signal r

k

is

u_{k}^{0} (t - r_{k}) = \sum_{m = - \infty}^{\infty} u_{k} (m) sinc (\frac{t - r_{k}}{T_{e}} - m) = \sum_{m = - \infty}^{\infty} u_{k} (m) sinc (\frac{t}{T_{e}} - τ_{k} - m)

The delayed digital signal u

k

ret

(n) is then obtained from:

u

k

ret

(

n

)=

u

k

0

(

nT

e

−r

k

)

that is:

u_{k}^{rct} (n) = \sum_{m = - \infty}^{\infty} u_{k} (m) sinc (n - τ_{k} - m) = h_{ideal, τ k} (n) * u_{k} (n)

where the symbol * symbolizes the convolution operation and where h

ideal,τk

(n)=sinc (n−

τk

).

The transfer function in the frequency domain of the filter h

ideal,τk

is

h_{ideal, τ k} (f) = ⅇ^{- j 2 π \frac{f}{f_{e}} τ_{k}} .

In theory, there is therefore an ideal filter h

ideal,τk

enabling the delayed sampled signal to be obtained from the sampled signal u

k

.

A filter of this kind cannot be implemented in practice because:

the summation is infinite, and

the filter is non-causal, i.e. it would be necessary to know the future samples to calculate the result.

The non-causal aspect means that it is only possible to approximate a version having an additional constant delay τ

0

denoted h

ideal,τ+τ0

(n).

In practice, and as shown in

FIG. 3

, a known process is used to obtain a filter h

τ

, having a finite number L of coefficients and constituting a good approximation of the filter h

ideal,τ+τ0

.

Finite impulse response filters are used in the most widespread applications. It is more difficult to design infinite impulse response filters and they require greater computation accuracy because the data is looped and computation errors can accumulate. Moreover, in a non-static environment, such as antenna processing, the long memory of the filters can be a problem: for example, it is difficult to change the pointing direction of an antenna suddenly.

All-pass filters, a subset of infinite impulse response filters, have the beneficial property of a constant modulus. The absence of amplitude error is unfortunately accompanied by an increased phase error compared to finite impulse response filters of the same complexity. To avoid these problems only finite impulse response filters are used.

Various methods are available for calculating the coefficients h

τ

(n) of a finite impulse response filter h

τ

approximating h

ideal,τ+τ0

.

When a finite number L of coefficients of the FIR filter has been chosen, a first method entails choosing a set of coefficients h

τ

(n) with n varying in the range from 0 to (L−1) which minimizes the mean square error ∫|E

τ

(f)|

2

df, where

E_{τ} (f) = \sum_{n = 0}^{L - 1} h_{τ} (n) ⅇ^{- j 2 π \frac{f}{f_{e}} τ} - H_{ideal, τ + τ_{0}} (f)

f

e

being the sampling frequency and the integration being formed over a wanted frequency band.

E

τ

(f) is often used to define cost functions used to assess the quality of a filter. This method is referred to as method M1. Its result is a truncated and delayed version of the ideal filter h

M1,τ

(n)=sinc(n−τ

0

−τ).

A second method entails first calculating the coefficients h

τ

(n) using the previous method and then multiplying them by an apodization window w(n). This method is referred to as method M2. Its effect is to smooth the frequency response of the filter.

If h

M1,τ

(n) are the coefficients of the filter obtained by the previous method, the coefficients h

M2,τ

(n) are obtained from the equation h

M2,τ

(n)=h

M1,τ

(n)·w(n) where, for example,

w (n) = 0.54 - 0.46 \cos (\frac{2 π n}{L - 1}) w (n)

is known as the Hamming window.

A third method entails adopting the set of coefficients h

τ

(n) which minimize the generalized mean square error C

τ

=∫W(f)|E

τ

(f)|

2

df, where W(f) is a frequency domain weighting function which can be modulated iteratively to obtain a satisfactory function h

M

3

τ

(n) and where the integration is performed over a wanted frequency band. This method is referred to as method M3. It yields a linear system of equations.

A fourth method entails adopting the coefficients h

τ

(n) for which

{(\frac{δ^{m} E_{τ} (f)}{δ f^{m}} &RightBracketingBar;)}_{f - 0} = 0,

where m varies in the range from 0 to (L−1). This method is referred to as method M4. The coefficients h

M4,τ

(n) are obtained from:

h_{M4, τ} (n) = \prod_{\underset{k \neq n}{k - 0}}^{L - 1} \frac{τ - k}{n - k}

Many examples and a comparison of infinite impulse response and finite impulse response filters can be found in an article in the IEEE signal processing magazine, January 1996, “Splitting the unit delay, tools for fractional delay filter design”, by Laakso, Välimäki, Karjalainen, Laine.

FIG. 3

shows that each input signal u

k

is then submitted to a delayed filter h

τ

i

(n), denoted RF

k

, where τ

k

is the normalized delay specific to the signal u

k

and is one of the values i/D where i varies in the range from 0 to (D−1).

The device shown in

FIGS. 1

to

3

therefore performs delay operations on M signals with interpolation filters RF

1

, . . . RF

M

of length L, which represents C

1

=M×L operations, where an operation is defined as a multiplication followed by an addition. Also, the memory needed to store the L coefficients of the D filters is C

2

=D×L words.

Such devices are not totally satisfactory. If the number of signals to be delayed is large, the computation load C1 and the memory capacity C2 needed to store the coefficients of the interpolation filters are so large that they make the devices extremely costly to implement and therefore compromise their economic viability.

BRIEF SUMMARY OF THE INVENTION

The aim of the invention is to propose a digital processing device which performs the same processing with the same accuracy but with greatly reduced computation load and memory requirement.

The above aims are achieved in accordance with the invention by a device for processing digital signals conveyed between a common port and M particular ports including a linear processing module characterized in that it includes particular filters which apply at the same time a phase-shift and a chosen amplitude emphasis at each of the M particular ports, said phase-shift being specific to each particular port and the amplitude emphasis being substantially the same for all the particular ports, and in that it includes at said common port a common inversion filter which performs the inverse amplitude acceleration through said chosen amplitude emphasis.

According to the invention, the common port can be an input port or an output port.

BRIEF DESCRIPTION OF THE DRAWINGS

Other features, aims and advantages of the present invention will become apparent on reading the following detailed description with reference to the accompanying drawings, which are given by way of non-limiting example and in which:

FIG. 1

is a functional block diagram of a prior art digital system applying specific fractional delay and linear processing operations to input signals;

FIG. 2

shows the same device, emphasizing the fact that the fractional delay filters use incremental delays;

FIG. 3

shows the same device as

FIG. 1

, emphasizing the frequency domain transfer functions of the ideal delay filters approximated in this prior art device;

FIG. 4

is as functional block diagram of a digital processing device of the present invention;

FIG. 5

shows the same device as

FIG. 4

, emphasizing approximated target filters;

FIG. 6

shows a device of the invention similar to the device shown in

FIGS. 4 and 5

, including phase-shift filters with four coefficients;

FIG. 7

shows maximum differences between Moduli for several fractional delays;

FIG. 8

shows phase delays each including the phase of a filter and the sampling frequency;

FIG. 9

shows phase errors, and constitute another representation of the phase-shift;

FIG. 10

represents a prior art device for providing specific output signals, each resulting from linear processing followed by fractional delays followed by respective filters; and

FIG. 11

shows a digital processing device according to the present invention, in which ay single input signal is processes to produce several output signals similar to those obtained with the device of FIG.

10

.

DETAILED DESCRIPTION OF THE INVENTION

FIGS. 4 and 5

show a device of the invention which processes M input signals u

1

, . . . , u

M

with the same result US as that produced by the device shown in

FIGS. 1

to 3 and with the same computation accuracy.

In the above device, each of the M input signals u

k

is passed first through a first phase-shift filter h

τ

d

(n) denoted FP

k

in

FIGS. 4 and 5

.

The phase-shift filters FP

1

, . . . FP

M

are finite impulse response filters of length L

1

, each approximating a target filter, and whose coefficients are calculated by a conventional method like the method M3 previously indicated. The target filters approximated in this way are not ideal delay filters as in the previous case. The target filters are defined in the frequency domain by

H_{target, τ} (f) = &LeftBracketingBar; H_{target} (f) &RightBracketingBar; ⅇ^{- j 2 π \frac{f}{f_{e}} (τ + τ_{0})}

where (τ+τ

0

) is the required non-integer normalized delay, being an additional constant delay equal to an integer number of sampling periods, and H

target

(f) depends only on the frequency.

Hereinafter, the expression “normalized fractional delay” refers to the non-integer part of the normalized delay actually introduced.

The function |H

target

(f)| is the same for all the phase-shift filters FP

1

to FP

M

. The function |H

target

(f)| therefore introduces the same input amplitude emphasis in each input filter.

Each input signal u

k

, after passing through the phase-shift filter h

τk

d

denoted FP

k

in

FIGS. 4 and 5

, corresponding to the required fractional delay τ

k

for that signal u

k

, is subjected to linear processing TT(k) by a linear processing filter T

k

specific to the input signal u

k

. The signals are then added.

Immediately downstream of this summation is an equalization filter of amplitude I, also referred to hereinafter as the common inversion filter. The role of the filter I is to perform the operation which is the inverse of the amplitude emphasis |H

target

(f)| introduced by the phase-shift filters h

τ

k

d

denoted FP

k

. In a preferred embodiment of the invention, this is done by approximating a filter with the frequency domain transfer function

H_{equal} (f) = \frac{1}{&LeftBracketingBar; H_{target} (f) &RightBracketingBar;} .

In many cases, the main role of the device is to put a signal back into phase or noise into antiphase, to reinforce the wanted signal and reduce the noise. In these cases, the interrelationship of the phases and of the moduli of the different branches is important.

However, zero overall distortion is not always necessary. Also, the equalization filter I can be included in another filter which is required for other purposes, such as equalizing a transfer function, for example that of the sensors.

Here the common inversion filter I is an FIR filter. Its coefficients can be calculated by one of the methods M1 through M3 previously described, or by any conventional filter design algorithm.

A linear filter like thecommon inversion filter I produces the same effects on the output signal US whether it is applied to each of those signals on the upstream side of the summation S or to the summed signal.

The same processing is therefore applied to the signals between M input ports P

1

to P

M

and the common output port C if a single common inversion filter I is placed after the summation S or if M filters identical to the common inversion filter I are placed upstream of the summation S.

We know that for a given input signal two linear filters in series supply the same output signal regardless of their order.

A common inversion filter like that shown in

FIGS. 4

to

6

and previously described therefore has the same effects on the output signal US of the entire device whether it is placed just upstream or just downstream of each of the specific linear processing operations T

1

to T

M

.

Clearly, then, the common amplitude inversion filter I, placed downstream of the summation S, produces the same results as if M identical filters were placed between each particular phase-shift filter h

τ

k

d

denoted FP

k

and the corresponding linear processing filter TT(k) denoted T

k

.

In the light of the above explanations, each combination of a phase-shifter filters FP

k

/and a common inversion filter I as shown in

FIG. 4

clearly therefore has the same effect on the output signal US as if each combination were contiguous at the location of each phase-shift filter FP

k

.

The transfer function of each of the contiguous pairs being equal to the product of the transfer functions of the phase-shift filter FP and the amplitude inversion filter I, each combination therefore constitutes a filter whose frequency domain transfer function is close to that of the ideal delay filter with transfer function

&LeftBracketingBar; H_{target, τ} (f) &RightBracketingBar; = ⅇ^{- j 2 π \frac{f}{f_{e}} (τ + τ_{0})} .

In one embodiment of the invention, the inversion filter I can also be an infinite impulse response filter. Because the equalization filter or amplitude inversion filter I does not vary with time, there is no “non-static” problem in this case.

If L

1

denotes the number of coefficients of the phase-shift filters FP

1

to FP

M

and L

2

denotes the number of coefficients of the equalization filter I, the device of the invention therefore performs C

1

=M·L

1

+L

2

operations, an operation being defined as a multiplication followed by an addition.

The device stores the coefficients of D filters of length L

1

and one filter of length L

2

, which represents a memory volume of C

2

=D·L

1

+L

2

words.

The number M of input signals is in practice very large.

For a given reference fractional delay filter h

τ

(n) the same phase-shift τ

1

to the sum of the signals.

If respective linear processing TT(1′), TT(1″), TT(1′″) is required for each input signal u

1

′, u

1

″, u

1

′″, it is then necessary to place the linear processing filters T

1

′, T

1

″, T

1

′″ upstream of the adder S

2

, applied to each input signal u

1

′, u

1

″, u

1

′″.

This produces a device in which there are no duplicated identical phase-shift filters. Thus the same operations are not performed at two different locations and the same set of coefficients is not stored twice.

A practical embodiment of a device of the kind shown diagrammatically in

FIG. 6

will now be described.

In this numerical example, the wanted input signal frequency band is 0-7 kHz. The sampling period f

e

is 16 kHz. The number D of filters to be implemented, also referred to as the resolution, is 20. The number of inputs is M=32. The length of the phase-shift filters is L1=4 coefficients and the length of the common inversion FIR filter is L

2

=11 coefficients.

The four coefficients of each phase-shifter filter are calculated by the generalized least squares algorithm, or method M3, and the seven coefficients of the amplitude inversion filter by a conventional linear phase finite impulse response filter design algorithm. The coefficients of the filters are normalized for calculation on a 16-bit fixed point DSP. The coefficients are as follows:

τ_{i} = \frac{i}{D}

h

i

(0)

h

i

(1)

h

i

(2)

h

i

(3)

0

0

14520

32767

14513

−19

1

0.05

13290

32697

15779

198

2

0.1

12089

32504

17067

449

3

0.15

10926

32189

18368

741

4

0.2

9806

31756

19672

1078

5

0.25

8737

31208

20968

1466

6

0.3

7723

30550

22246

1909

7

0.35

6769

29789

23494

2410

8

0.4

5878

28931

24703

2973

9

0.45

5052

27985

25861

3600

10

0.5

4293

26958

26958

4293

The values for i=11, . . . , 19 are calculated using a relationship of symmetry h

i

(n)=h

D−i

(L

1

−1−n) in the present case where the additional delay τ

0

is equal to (L

1

−2)/2.

This property will now be demonstrated:

An FIR filter of length L, can take non-zero values for n=0, . . . , L

1

−1. The point at the center of the impulse response is therefore (L

1

−1)/2. Delays close to this value are easier to obtain than others. This is why a preferred embodiment of the invention seeks to approximate a normalized ideal phase-shift of τ+τ

0

where τ

0

is fixed at (L

1

−2)/2 and the normalized fractional delay τ varies in the range from 0 to 1, which gives a range of normalized delays that varies in the range from (L

1

−2)/2 to L

1

/2 and is symmetrical about the center (L

1

−1)/2.

We also know that {tilde over (H)}(f)=H*(f)e

−j2πf(L1−1)

(where H*(f) is the complex conjugate of H(f)) if {tilde over (h)}(n)=h(L

1

−1−n).

H

τ

(f) approximates

\begin{matrix} &LeftBracketingBar; H_{target} (f) &RightBracketingBar; \cdot H_{ideal, τ + τ 0} (f) = &LeftBracketingBar; H_{target} (f) &RightBracketingBar; ⅇ^{- j2π f (τ + τ_{0})} \\ = &LeftBracketingBar; H_{target} (f) &RightBracketingBar; ⅇ^{- j2π f (τ + (L_{1} - 2) / 2)} \end{matrix}

in which {tilde over (H)}

τ

(f)=H*(f)e

−j2πf(L1−1)

approximates

\begin{matrix} &LeftBracketingBar; H_{target} (f) &RightBracketingBar; \cdot H_{ideal, τ + τ 0}^{*} (f) ⅇ^{- j2π f (L_{1} - 1)} & = & &LeftBracketingBar; H_{target} (f) &RightBracketingBar; \cdot ⅇ^{+ j2π f (τ + τ_{0})} ⅇ^{- j2π f (L_{1} - 1)} \\ = & &LeftBracketingBar; H_{target} (f) &RightBracketingBar; \cdot ⅇ^{- j2π f (L_{1} - 1 - (L_{1} - 2) / 2 - τ)} \\ = & &LeftBracketingBar; H_{target} (f) &RightBracketingBar; \cdot ⅇ^{- j2π f ((L_{1} - 2) / 2 + 1 - τ)} \\ = & &LeftBracketingBar; H_{target} (f) &RightBracketingBar; \cdot H_{ideal, (1 - τ) + τ0} (f) \end{matrix}

so that {tilde over (H)}

τ

(f) therefore approximates H

1−τ

(f) and {tilde over (H)}

1+τ

(f) approximates H

τ

(f), whence h

τ

is substantially equal to {tilde over (H)}

1+τ

(f), in other words: of length L, the inventors have found after much work that it is possible to produce a combination of a phase-shift filter FP

M

and an amplitude inversion filter I in which the phase-shift filter FP

k

has a number L

1

of coefficients which is less than L, or even very small compared to L, and such that the computation accuracy of the combination of filters is similar to that of the reference filter h

τ

(n).

A smoothed emphasis function |H

target

(f)| which has low values at half the sampling frequency is preferably chosen for this purpose.

These functions |H

target

(f)| are advantageously chosen from functions of the form

&LeftBracketingBar; H_{target} (f) &RightBracketingBar; = X + Y (1 + \cos (\frac{2 π f}{f_{e}})),

the ratio X/Y preferably being in the range from 0 to 0.3.

To normalize the value of |H

target

(f)| at f=0, X and Y preferably satisfy the condition X+2Y=1.

A preferred choice of the function |H

target

(f)| is

&LeftBracketingBar; H_{target} (f) &RightBracketingBar; = 0.06 + 0.47 (1 + \cos (\frac{2 π f}{f_{e}})) .

With this function |H

target

(f)|, the coefficients of the filters are advantageously calculated by method M3 and the function

W (f) = 1 for f < \frac{6}{16} f_{e};

W (f) = 5 for \frac{6}{16} f_{e} \leq f < \frac{7}{16} f_{e};

W (f) = 0.01 for f > \frac{7}{16} f_{e} .

Thus L

1

<L, or even L

1

<<L and M>>1, which implies:

M·L

1

+L

2

<<M·L

and

D·L

1

+L

2

<<D·L.

The computation load and the memory volume needed are therefore very small.

The computation load and the volume of memory needed can be made smaller if, as shown in

FIG. 4

, input signals u

1

′, u

1

″, u

1

′″ are to be delayed with the same fractional delay τ

1

.

An adder S

2

to sum the signals is then placed on the upstream side of a phase-shift filter FP

1

which applies h

τ

(n)=h

1−τ

(L

1

−1−n) regardless of the value of n in the range from 0 to L

1

−1.

This relationship of symmetry means that only the coefficients of half the necessary filters have to be stored in memory.

It is sufficient to store, as in the example explained here, the coefficients of the filters h

τ

corresponding of normalized delays (τ+τ

0

) in the range from

\frac{(L_{1} - 2)}{2} to \frac{(L_{1} - 1)}{2} .

The filters h

1

T corresponding to normalized delays 1−τ+τ

0

symmetrical to those delays τ+τ

0

about (L

1−

1)/2 are deduced using the relationship of symmetry h

τ

(n)=h

1−τ

(L

1

−1−n). In this way, for two phase-shift filters whose respective normalized fractional delays are symmetrical about a fractional delay equal to half the sampling period, the coefficients of one of the two filters are deduced from those of the other one using the abovementioned relationship.

For the finite impulse response amplitude inversion filter, we obtain for h

equal

(n) varying in the range from 0 to L

2

−1=10: −217, 747, −2914, 8225, −18117, 32767, −18117, 8225, −2914, 747, −217.

If the phase must not be strictly linear, a third order infinite impulse response filter can be used with L2=7 coefficients whose zeros are: −0.655 and 0.0016±0.0383i, and whose poles are: −0.7089 and −0.5891±0.2122i. This filter is calculated using the well-known Yule-Walker algorithm, for example.

The new structure is equivalent to a conventional structure with interpolation filters h

total

, where h

total

(n)=h

i

(n)*h

equal

(n), * being the convolution operation.

In the case of the above device, the computation load is therefore C

1

=M·L

1

+L

2

=32×4+7=135 operations and the memory required is C

2

=D·L

1

+L

2

=20×4+7=87 words.

FIG. 6

shows the structure of this filter in the form of a functional block diagram.

In the case of antenna processing, the delays R

i

with i varying in the range from 1 to M, for example, are chosen as a function of the antenna pointing angle.

The delays are then normalized and quantized. This yields the integer parts di and the fractional parts r

i

of the normalized delays:

r_{nq \cdot i} = \frac{round (D \cdot Ri \cdot f_{e})}{D},

d

i

=floor(r

nq·i

) and r

i

=r

nq·i

−d

i

, where the round and floor fractions are respectively the integer approximation and integer part functions.

The resulting digital device has been compared with a prior art filter like that shown in

FIGS. 1

to

3

, in which the coefficients of the fractional delay filters were calculated using:

method M1 and L=20 (C1=640 operations, C2=400 words),

method M2 with a Hamming window and L=18 (C1=576, C2=360),

method M2 with a Blackman window and L=24 (C1=768, C2=480),

method M3 with L=12 and W(f)=1 in the wanted band 0-7 kHz (C1=384, C2=240),

method M3 with L=13 and W(f)=1 in the wanted band 0-7 kHz (C1=416, C2=260),

method M4 with L=30 (C1=960, C2=600).

The results obtained with a device of the invention are also shown, in which device the phase-shift filters have six coefficients and the filters are calculated using the same methods as for the filter of the invention with four coefficients previously described.

The ideal filters have a phase φ which increases linearly with the frequency f. The fraction

\frac{- φ}{(2 π f)},

referred to as the “phase delay”, is therefore constant in the ideal case.

FIG. 7

shows the maximum differences between the moduli, i.e. max

τ

1

τ

2

(||H

τ

1

(f)|−|H

τ

2

(f)||) for D=20 fractional delays.

It enables the quality of the phase-shift filters of the new method to be assessed. The amplitude inversion filter does not modify the difference between the moduli, i.e. the same results are obtained for H

total

(f)=H

τ

d

(f)H

equal

(f) and for H

τ

d

(f) alone. This figure can therefore be used to assess the quality of the phase-shift filters H

τ

d

.

The order of the other filters has been chosen to obtain a maximum amplitude error of around 1 dB, which corresponds to a factor of 1.12. In the worst case scenario, there is therefore a risk of having a residue due to the amplitude error of 0.12 (−18 dB) instead of cancellation for destructive interference.

It can be shown that the performance of all the filters is excellent for low frequencies up to 6 kHz, except for method M1. The maximum error is always at the maximum frequency of 7 kHz. The two examples of the new method have the lowest maximum error. The overall performance is satisfactory for the filter with four coefficients and excellent for the filter with six coefficients, which shows the better results from all the filters.

FIG. 8

shows the phase delays

- \frac{φ_{τ} (f)}{2 π f} f_{e},

where φ

τ

(f) is the phase of the filter h

τ

d

and f

e

is the sampling frequency. It is difficult to assess the quality of the phase from this figure because errors at low frequencies are emphasized by the denominator and because it is not easy to compare two methods. It is included here to show that a filter of even length can be exact for the phase for integer and integer-and-a-half delays (referred to as linear phase filters) whereas an odd length filter cannot be exact for an integer-and-a-half delay. This is the main reason why filters of odd length are rarely used.

FIG. 8

also shows that at the maximum frequency the delays tend towards integer values except for the integer-and-a-half delay. This behavior is due to the fact that the discrete Fourier transform of a real sequence is always real at the maximum frequency. It follows from this that the phase is zero or π or undefined because the modulus is zero. This explains why the phase of the filters departs from the target value at high frequencies.

FIG. 9

shows the phase error [(φ

τ

(f)−φ

τ,ideal

(f)], in other words this is another representation of the phase-shift. It is well suited to comparing the quality of the phase of the filters. The greatest phase errors still occur at the maximum wanted frequency of 7 kHz. The filter with four coefficients is subject to slightly greater phase error than most other filters. On the other hand, the new method has the smallest two maximum errors.

FIG. 9

can also be used to assess the accuracy with which the delay must be quantified. There is no benefit in storing a large number of filters if the filters are insufficiently accurate. At the maximum frequency, a delay of one sample corresponds to a phase-shift of 180°. Given the maximum phase error of approximately 3°, it is therefore reasonable to store more than 60 filters so that the delay quantizing error is less than the error introduced by the filter. Performance is degraded if the available memory is unable to store this many filters. This situation can arise in typical real time environments, for example in respect of the internal memory of a DSP. In this case, the new structure produces better results with low complexity.

Lagrange filters have excellent performance at low frequencies but the number of coefficients increases rapidly with the wanted frequency band. To obtain good behavior in the 0-7 kHz band, an excessive number of coefficients is required.

In the worst case scenario (excluding Lagrange), the maximum phase difference is equal to 7°, which yields the same residual error of 0.12 (−18 dB) as that due to the amplitude error.

Thus even the critical situation of destructive interference can be processed with only four operations per input and four coefficients to be stored for each fractional delay. Comparison with conventional interpolation filters (methods M1 to M3) show that the following can be obtained simultaneously:

a reduction of the computation load by a factor in the range from 3 to 8,

a reduction in the memory requirement by a factor in the range from 3 to 8,

with a smaller maximum error and good overall performance.

Because the interpolation is often the most costly processing operation, the method is particularly suitable for forming lobes for antenna processing. The figures with L1=6 show that this method is also very effective for applications that demand a very high accuracy.

Note that the above two examples of filters cover most existing applications of fractional delay filters.

It is rare to need delay filters having less than four coefficients and there is rarely a need for an accuracy better than that obtained with the second example, for which L

1

=6.

The remainder of the description covers one embodiment of a digital filter of the invention and conforming to that described above, in which the parameters L

1

, L

2

and D are respectively equal to 4, 19 and 60. This example is given by way of illustration only, in the form of a program written in C. The coefficients of the filter routine 1.C were calculated in advance using the previously mentioned functions |H

target

(f)| and W(f).

The purpose of the prior art device shown in

FIG. 10

is to provide M specific output signals V

k

(with k varying in the range from 1 to M) each resulting from linear processing TT(k) performed by a filter T

k

followed by a fractional delay formed by a filter RF

k

, based on the same input signal VE.

This problem arises in fractional pitch resolution voice coders, for example. Several voice coders use the following equations to determine the pitch P:

\begin{matrix} x = (x (n), x (n - 1), \dots, x (n - N + 1)), \\ y = (x (n - P), x (n - P - 1), \dots, x (n - P - N + 1)) \\ P = argmax \frac{⟨ x, y ⟩}{{&LeftDoubleBracketingBar; y &RightDoubleBracketingBar;}^{2}} \end{matrix}

To obtain a fractional resolution of the pitch P, the signal x must be delayed by a fractional delay.

As shown in

FIG. 10

, the device includes on the input side a device D which first duplicates the input signal VE in the form of M signals VE′ identical to it, after which each of the signals VE′ is passed through a linear processing filter T

k

applying linear processing TT(k), followed by a fractional delay filter RF

k

introducing a fractional delay r

k

.

FIG. 11

shows a digital processing device of the present invention in which a single input signal VE is processed to produce M output signals V

1

to V

M

similar to those obtained with the device from FIG.

10

. The components are similar to those of the device shown in

FIGS. 4

to

6

, that is to say phase-shift filters h

τ

k

d

(n) denoted FP

1

, FP

2

, . . . , FP

M

, where k varies in the range from 1 to M, linear processing filters T

k

, where k varies in the range from 1 to M, and an amplitude inversion filter I.

The amplitude inversion filter I applies amplitude pre-emphasis to the input signal VE. To be more precise, the transfer function of the filter I has a non-constant modulus and a zero phase. The resulting signal VE″ is then duplicated by a device D to produce M identical signals VE′″. Each of the signals VE′″ is subjected to specific linear processing TT(k) (with k varying in the range from 1 to M) by a filter T

k

and is then passed through a specific phase-shift filter h

τ

k

d

denoted FP

k

.

Each phase-shift filter FP

k

is an approximation of a target filter whose transfer function has the following features:

Its phase is close to that of the ideal fractional delay filter introducing the delay r

k

required for the output signal, i.e. its phase is close to −2πr

k

f.

Its modulus is the inverse function of the modulus of the amplitude inversion I on the input side of the device. As in the device shown in

FIGS. 3 and 4

, the amplitude inversion filter clearly has the same effect as if it were applied to each duplicated signal just upstream of the phase-shift filters h

τ

k

d

denoted FP

k

, with k varying in the range from 1 to M.

In

FIG. 10

, as in

FIGS. 1 and 2

, the computation load is C

1

=M·L operations. For the same reasons, D filters of length L are stored, i.e. C

2

=D×L words in memory, whereas in the case of the device from

FIG. 11

the computation load is C

1

=M·L

1

+L

2

operations and the memory required is D·L

1

+L

2

words, where L

1

is the length of the phase-shift filters and L

2

is the length of the amplitude inversion filter.

In the same manner as previously, an amplitude pre-emphasis filter I/phase-shift filter FP

k

combination having the same performance as a fractional delay filter RF

k

of length L can be obtained in which the length L

1

of the phase-shift filter FP

k

is less than, or even very small compared to, the length L of the delay filter RF

k

. The device shown in

FIG. 11

then produces the same result with the same accuracy as the device from FIG.

10

and with a greatly reduced computation load and memory requirement.

Of course, the invention is not limited to the embodiments described. It encompasses any variant in accordance with the spirit of the invention.

In particular, the common amplitude emphasis filter I can also be used to equalize a transfer function, for example that of a sensor.

More generally, the invention is not limited to interpolation filters but can be used for other linear processing with either a large number of inputs or a large number of outputs.

Program in C

main.c program page 1

/* application to beamforming algorithm */

#include “cons.h”

#include <stdio.h>

#include <math.h>

/* antenna geometrie */

const float xposition[] = {0,.1,.2,.3,.4,.5,.6,.7,.8,.9};

const float xmin=0;

const float xmax=.9;

/* filter1 for phase shift introducing amplitude distortion */

/* filter2 for amplitude egalisation */

long filter1 (short*, int, struct delay*, float*);

int filter2(long*)

struct delay *calcdelay(float,int);

int main(int argc, char*argv[])

{

/************************************************/

/* declaration */

struct delay *del;

/* integer and fract. part*/

char *infile,*outfile,*parfile; /*

file names */

int i,m;

/* loop variables */

int delmax,base;

/* max. delay, address offset */

int neof,

/* not_end_of file */

long out1;

/* output filter 1 */

short out2;

/* output filter 2 */

float phi

/* steering angle */

float shaping[M];

/* microphone weighting */

short indata[INBUFFER] [M];

/*input buffer */

long outdata[L2];

/* buffer for second filter */

FILE *infpt,*outfpt,*parfpt;

/*file pointer */

/************************************************/

/* get command line parameters infile, outfile, parameter file*/

if (argc!=4)

{

printf(“use frac in_file out_file

parameter_file\n”)

return(−1 );

}

infile=argv[1];outfile=argv[2]; parfile=argv[3];

infpt=fopen(infile, “rb”) ;outfpt=fopen(outfile, “wb”);parfpt=fopen

(parfile,rb”);

if ((infpt==NULL) | (outfpt==NULL) | (parfile==NULL))

{

printf(“unable to open input or output or parameter file”);

return(−2);

}

/************************************************/

/* maximum delay */

delmax =(int) ((xmax-xmin) *Fs/c);

if (delmax+L1 >INBUFFER)

{

printf(“Input Buffer too small”);

return(−3);

}

main.c program page 2

/************************************************/

/* get multiplexed input data from infile and put into indata */

fread(indata,M*sizeof(short),(delmax+L1),infpt);

base=delmax+LI;

/************************************************/

/* main loop for each sample */

for (neof= 1;neof>0;)

{

/* read input data and steering parameters */

base=( (base+1)%INBUFFER);

fread(&indata[base] [0],M*sizeof(short),1,infpt)

neof=fread(&phi,sizeof(float),1,parfpt)

fread(shaping,sizeof(float),M,parfpt);

/* compute integer and fractionnal delay for each

microphone */

del=calcdelay(phi,delmax);

/* beamforming by delay-weighting-sum */

out1=filter1(&indata[0] [0],base,del,shaping);

/* the second filter shift loop */

for (i=O;i<L2−1;i++)

outdata [i]=outdata1 [i+1];

outdata{L2− 1)=out 1;

out2=filter2(outdata);

fwrite(&out2,sizeof(short),1,outfpt);

}

return(0);

}

filter1.c program page 1

/* performs delay - weighting - sum beamforming */

/* Wolfgang Tager, 28.10.96 */

#include “cons.h”

#include “math.h”

float abs(float); /* since default is int */

long filter1(short indata[INBUFFER] [M],int base,struct delay

*del,float shaping [M])

}

static int f1coeff[RESOL] [L1] =

{

{14520

,32767

,14513

,−19},

{14108

,32757

,14932

,50},

{13697

,32734

,15354

,122},

{13290

,32697

,15779

,198},

{12886

,32646

,16206

,277},

{12486

,32582

,16636

,361},

{12089

,32504

,17067

,449},

{11697

,32412

,17500

,541},

{11309

,32307

,17934

,639},

{10926

,32189

,18368

,741},

{10547

,32058

,18803

,848},

{10174

,31913

,19238

,960},

{9806

,31756

,19672

,1078},

{9444

,31586

,20106

,1202},

{9087

,31403

,20538

,1331},

{8737

,31208

,20968

,1466},

{8392

,31001

,21397

,1607},

{8055

,30781

,21823

,1755},

{7723

,30550

,22246

,1909},

{7398

,30308

,22666

,2069},

{7080

,30054

,23082

,2236},

{6769

,29789

,23494

,2410},

{6465

,29514

,23902

,2590},

{6168

,29228

,24305

,2778},

{5878

,28931

,24703

,2973},

{5596

,28625

,25095

,3174},

{5320

,28310

,25481

,3384},

{5052

,27985

,25861

,3600},

{4792

,27651

,26234

,3823},

{4539

,27309

,26600

,4054},

{4293

,26958

,26958

,4293},

{4054

,26600

,27309

,4539},

{3823

,26234

,27651

,4792},

{3600

,25861

,27985

,5052},

{3384

,25481

,28310

,5320},

{3174

,25095

,28625

,5569},

{2973

,24703

,28931

,5878},

{2778

,24305

,29228

,6168},

{2590

,23902

,29514

,6465},

{2410

,23494

,29789

,6769},

{2236

,23082

,30054

,7080},

{2069

,22666

,30308

,7398},

{1909

,22246

,30550

,7723},

{1755

,21823

,30781

,8055},

{1607

,21397

,31001

,8392},

{1466

,20968

,31208

,8737},

{1331

,20538

,31403

,9087},

{1202

,20106

,31586

,9444},

{1078

,19672

,31756

,9806},

{960

,19238

,31913

,10174},

{848

,18803

,32058

,10547},

{741

,18368

,32189

,10926},

{639

,17934

,32307

,11309},

{541

,17500

,32412

,11697},

{449

,17067

,32504

,12089},

{361

,16636

,32582

,12486},

{277

,16206

,32646

,12886},

{198

,15779

,32697

,13290}

{122

,15354

,32734

,13697}

{50

,14932

,32757

,14109}

};

filter1.c program page 2

long filtout,

/* integer calculation for fast exact results

no overflow for a single microphone

since sum_over_i(abs(f1 coeff[r] [i]))<65536 */

float sumout=O;

int m,i,firstel; /*loop counters, first element for fract. filter*/

/*********************************************/

/* delay - weighting - sum

*/

for (m=O;m<M;m++)

{

filtout=O,

firstel=base-del[m] .integer,

/* integer delay compensation */

for (i=O;i<L1;i++)

filtout+=f1coeff[del[m] .frac][i]*indata

[(firstel-

i+INBUFFER)%INBUFFER] [m];

sumout=sumout+shaping [m]*filtout;

}

/*********************************************/

/* normalisation and saturation if necessary

(does not occur if sum(abs(shaping))<M

*/

sumout=(sumout/M);

sumout=sumnout/(11<<16);

if (abs(surnout)>32767)

{

printf(“\nsumout filter 1 exceeds 16 bits: %f ”,sumout);

printf(“\nlimited to saturation value ”);

if (sumout>O)

sumout=32767;

else

sumout=−32767;

end

getchar();

}

return( (long) sumout );

}

calcdelay.c, filter2.c and cons.h programs

/* calculate integer and fractionnal delay for plane wave

*/

/* impinging on a linear array under an angle phi

*/

/* can be adapted to other arrays and/or nearfield equations

*/

/* wolfgang Tager, 28.10.96

*/

#include “cons.h”

#include <math.h>

#include <stdio.h>

struct delay *calcdelay(float phi,int delmax)

{

extern float xposition[];

extern float xmax;

static struct delay del [M];

int m;

float tmp;

for (m=0;m<M;m++)

{

tmp=delmax- (xmax-xposition[m])*cos(phi)*Fs/c + .5 /RESOL;

del[m] .integer=(int) tmp;

del *[m] .frac=(int) (RESOL* (tmp-del 8 m] .integer));

}

return del;

}

/* amplitude egalisation filter */

/* wolfgang Tager, 28.10.96 */

#include “cons.h”

int filter2(long* outdata)

{

static int f2coeff[L2]=

{ −30,66,−189,512,−1246,2760,−5649, 10760,−19311,32767,−19311,

10760,−5639,2760,−1246,512,−189,66,−30};

char shift=17;

/* sum(abs(f2coeff)) < 2{circumflex over ( )}17*/

int i;

float filtout=0;

for (i=0,i<L2;i++)

filtout+=f2coeff[i]*outdata[i]

return( (int) (filtout/ (I1<<shift) ) );

}

#define

M

10

#define

INBUFFER

1000

#define

RESOL

60

#define

L 1

4

#define

L2

19

#define

Fs

8000

#define

c

340

struct delay

{

int integer,

int frac;

};

Number	Name	Date	Kind
5654909	Sun	Aug 1997	A
5691929	Sun	Nov 1997	A
5926455	Allpress	Jul 1999	A

Fractional delay digital filter

Information

Patent Number

Date Filed

Date Issued

Inventors

Original Assignees

Examiners

Agents

CPC

US Classifications

Field of Search

US

International Classifications

Abstract

Description

Claims

Priority Claims (1)

PCT Information

US Referenced Citations (3)

Non-Patent Literature Citations (3)

Entry
Cetin A E et al.: “Digital Interpolation BeamForming Using Recursive Filters” Journal of the Acoustical Society of America; vol. 85, No. 1 Jan. 1, 1989.
Sangil Park: “Real Time Pitch Shifting Techniques” ICCE Jun. 5-7, 1991.
T.I. Laakso et al. “Splitting the Unit Delay, Tools for Fractional Delay Filter Design” IEEE Signal Processing Magazine, vol. 13, No. 1, Jan. 1996.