REDUCTION OF DIGITAL FILTER DELAY

TECHNICAL FIELD

The present invention relates to reduction of delay in digital filters based on a desired magnitude response.

BACKGROUND

Requirements on a digital signal processing device to be used in a telephony- network normally aim at minimizing the signal propagation delay imposed by the device, to reduce the total round trip delay which otherwise would affect the conversation quality negatively, see [I]. It is common for signal processing devices (such as noise reduction devices, echo cancellers, etc.) to apply some kind of linear filtering, serially connected to the signal path, to improve the speech signal. Two examples are:

1. Noise reduction devices, where the common solution is to apply filtering so that the disturbing background noise is reduced while at the same time retaining the speech so as to improve intelligibility.

2. Echo cancellation devices, where typically a linear adaptive filter (connected in parallel with the echo path) models the echo path to enable subtraction of a replica of the echo and thus a reduction of the echo level. However, due to imperfections in the adaptive algorithms and due to time varying and possibly non-linear echo paths, residual echo is left after the subtraction. This residual echo is typically removed by a non-linear processor (NLP), which detects the presence of the residual echo and, if found, removes it (optionally replacing it with some kind of generated comfort noise). The NLP is usually implemented as a switch that either passes or suppresses the signal applied to it. How-

ever, the NLP may also be implemented as a serially connected linear filter, as described in [2].

In both these telecommunication applications it is vital that the voice en- hancement devices do not add to the round trip delay. Thus, it is desirable to minimize the delay of the signal processing filters. However, at the same time it is also desirable to maintain high stop-band rejection. Achieving high stop- band rejection is particularly important in echo cancellation applications, where a serially connected filter is required to suppress residual echo in cer- tain frequency bands (the ones where residual echo dominates over true near end signal), while passing some other bands (the ones dominated by the true near end signal).

The derivative of the phase characteristic of a linear filter with respect to fre- quency is a measure of the time delay (i.e. group delay or envelope delay) that signal frequency components undergo in passing through the filter. A minimum-phase filter characteristic implies a minimum delay function, while a maximum-phase filter characteristic implies that the delay characteristic is also maximum. A mixed-phase system has a delay in between the two extremes. One solution to decrease the delay, which is imposed by a linear filter, is to convert it to a minimum-phase filter, which by definition has a minimal group delay through the filter.

Suppose that we have a finite impulse response (FIR) filter prototype H{ω) with real coefficients. The magnitude square value of its frequency response is

Relationship (1) implies that if we replace a zero z _k of the filter by its inverse inside the unit circle, the magnitude characteristic of the filter does not change. Hence, a method for converting a finite impulse response (FIR)

mixed/ maximum phase filter to a minimum-phase filter is to mirror all zeros into the unit circle. However, this requires that the exact values of the zeros z _k need to be calculated and would also require a large amount of divisions, which is inefficient in a typical real-time application on a digital signal proc- essor (DSP). Mainly three workable other approaches exist for performing the conversion to a minimum-phase filter:

A. Spectral factorisation, e.g. using some approximating recursive algorithms as described in [3, 4]. However, a problem with spectral factori- sation is that it is computationally demanding and that it needs an increasing number of iterations to reliably form a minimum-phase filter of a demanding filter prototype (i.e. a complicated desired magnitude response with high dynamic range) with high required stop band attenuation.

B. Employing homomorphic filtering as described in, for example, [5]. The method essentially determines the real cepstrum of the magnitude response of the desired digital filter, transforms this into a corresponding complex cepstrum representing a minimum-phase filter, and de- termines the minimum-phase digital filter corresponding to the complex cepstrum. The problem with this method is that it requires long discrete Fourier transforms (DFTs) to attain high stop band rejection in the resulting minimum-phase filter and thus becomes computationally demanding.

C. The errors imposed due to the usage of a discrete Fourier transform (DFT) in method B have been analysed in [6], and a method is proposed to reduce these errors. The method basically trusts on the usage of long DFTs together with truncation of the resulting impulse re- sponses to minimise the effects of aliasing. This results in that the effective length of the impulse response will be half of the length of the

used DFT. When using short DFTs, this leads to difficulties in obtaining the desired magnitude response of the prototype filter.

SUMMARY

An objective of the present invention is to provide a method and apparatus for reducing the delay of a desired filter with a given magnitude response to approximately minimum-phase with shorter DFTs than in the prior art.

Another objective is an echo canceller including such an apparatus.

Still another objective is a noise suppressor including such an apparatus.

These objectives are achieved in accordance with the attached claims.

Briefly, the present invention is based on method B above, but shapes the coefficients to be included in the complex cepstrum corresponding to the minimum-phase filter by applying a smoothly decaying window function. This makes it possible to obtain estimates of minimum-phase filters from de- manding magnitude responses by using short DFTs. The use of short DFTs results in a low complexity solution compared to the prior art solutions, without sacrificing the requirements on high stop band rejection.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention, together with further objects and advantages thereof, may best be understood by making reference to the following description taken together with the accompanying drawings, in which:

Fig. 1 is a diagram illustrating the location of zeros of the magnitude response of a FIR filter;

Fig. 2 is a flow chart illustrating a prior art method for obtaining an estimate of a minimum-phase filter from a desired magnitude response;

Fig. 3 is a flow chart illustrating an embodiment of the method in accordance with the present invention;

Fig. 4(a)-(d) is a sequence of diagrams graphically illustrating the embodiment of Fig. 3; Fig. 5 is a diagram comparing the embodiment of Fig. 3 with the prior art method B;

Fig. 6 is a flow chart illustrating another embodiment of the method in accordance with the present invention;

Fig. 7 is a diagram comparing the embodiment of Fig. 3 with the em- bodiment of Fig. 6; and

Fig. 8 is a block diagram illustrating an embodiment of an apparatus for reducing digital filter delay in accordance with the present invention.

DETAILED DESCRIPTION

To better explain how the idea of the present invention is related to minimum-phase filter design, let us first reproduce some of the theory.

Suppose we have a finite impulse response (FIR) filter

h(n), O ≤ n ≤ N (2)

The system function of the filter h is given by its z-transform

H (ή = ∑h(n) Z-" (3)

Denote the logarithm of the magnitude response of the system function above as

(4)

The cepstrum, or more accurately the real cepstrum, of h is defined as the coefficients h(n) of the power series expansion of H(z), i.e.

given that the series is convergent. The coefficients of the power series expansion above, i.e. the real cepstrum, can be computed using the inverse z- transform of H{z)

h(n) = —j)H(z)z"- ^l dz (6) 2πy _c

or explicitly

Choosing the integration contour, C, to be the unit circle, (6) transforms to

h(n) = — ] H{e ^Jύ> )e ^jωn dω (8)

2π

Denote the zeros of H(z) by Ok. Now, it is well known that the zeros of

are located in conjugate reciprocal pairs, as illustrated in Fig. 1. Denote the zeros of H(z) that lie inside the unit circle by Ck and the zeros that lie outside the unit circle by bk. Then (9) can be rewritten as

Computing the logarithm of the square root of the above expression results in

H(z) = \nA ₊ -∑ln(l-c _k z- ⁱ ) ₊ ±∑\n{\~b _k z) (H)

* ^■ k=\ I Ar=I

Using the following power series representations

gives

Comparing the above with (5) and identifying coefficients of the same power of z it is evident that

From (14) it may be noted that the components originating from the zeros inside the unit circle only contribute to the cepstrum at positive n and that the components originating from zeros outside the unit circle only contribute to the cepstrum at negative n. A minimum-phase filter has all its zeros inside

the unit circle. Hence, to get the cepstrum corresponding to a minimum- phase filter one needs to set the cepstrum at negative n to zero.

Furthermore, from (14) it may also be noted that the real cepstrum has infinite duration although the original filter h(n) has a finite length.

A similar reasoning as above can be applied to the complex cepstrum (see [5]), i.e.

where it is noted that the magnitude (" | | ") taken in (7) has been removed. In the special case of a minimum-phase filter this leads to the following representation of the complex cepstrum

Comparing (16) with (14), it is evident that the complex cepstrum of a minimum-phase sequence can be restored from its real cepstrum as

In order to obtain the impulse response of the minimum-phase filter one needs to reverse the process, i.e. compute

(18) n =-oo

H _min (z)=ex _P {H _min (z)} (19)

The outlined method may be implemented by using the discrete Fourier transform (DFT) and its inverse (IDFT), as described in [5] and illustrated in Fig. 2. The discrete transforms are typically implemented by the fast and inverse fast Fourier transform (FFT and IFFT) of length M. Step S 1 determines the magnitude response J of a desired filter. Step S2 determines the real cepstrum corresponding to this magnitude response in accordance with

H(ω _k ) = 0,1,...,M-I (21) h(n) = IDFT{H{ω _k )), « = 0,1,..., M -I (22)

Here M must be sufficiently large to avoid aliasing (M= 1024 in the example given in [5]). Step S3 determines the complex cepstrum of the desired minimum-phase filter from the real cepstrum in accordance with

2- h(n), n = 1,2,..., M /2

*U(») = h{n), n = 0 (23)

0 « = M/2 + l,..., M -l

Note that due to the periodicity assumed by the IDFT, the cepstral components with negative n are (by convention) located between M/ 2+1 and M-I in equation (23). Finally step S4 determines the frequency response H^ω^ of the minimum-phase filter from this complex cepstrum in accordance with

k = 0,l,...,M-l (24 ⁾

The resulting filter h _min (n) is causal and gives an estimate of a minimum- phase filter, where the quality of the estimate depends on e.g. the involved FFT lengths.

When computing the minimum-phase impulse response, a problem arises in practice due to the fact that it is necessary to replace the infinite coefficient sequence h _n jn) with a finite sequence. This operation implies a truncation, which introduces errors in the computations. The prior art method B suggests that the problem can be relieved by selecting a long enough DFT. How- ever, as mentioned earlier, using long DFTs results in high computational complexity. The present invention takes another approach by noting that when truncating a function before computing its DFT, the resulting spectrum will be a convolution between the spectrum of interest and the spectrum of the window function used for truncation (Gibbs's phenomenon). Ac- cording to the present invention, truncation of the cepstrum is done by shaping the cepstrum with a shorter, smoothly decaying (in n, which is called "quefrency" in cepstrum analysis) window function instead of abruptly truncating it with a longer rectangular window.

Fig. 3 is a flow chart illustrating an embodiment of the method in accordance with the present invention. Steps Sl -S3 are the same as in the method in accordance with Fig. 3. However, in accordance with this embodiment of the present invention a new shaping step S3A has been inserted between steps S3 and S4. Step S3A applies a smoothly decaying window function to shape the coefficients to be included in the complex cepstrum. Thus, a window function w(n) is applied to h^ri) in (23) in accordance with

*«(») = w(«)W«) (25)

In step S4 h _shaped (n) replaces H _n Jn) in (24).

Fig. 4(a)-(d) is a sequence of diagrams graphically illustrating the embodiment of Fig. 3. Fig. 4 (a) illustrates a hypothetical real cepstrum. In practice a real cepstrum will of course not have this form, but this hypothetical real cepstrum is very convenient for illustrating the method in accordance with the present invention. Fig. 4(b) illustrates the complex cepstrum formed from the real cepstrum in accordance with (23). Fig. 4(c) illustrates a shaping window in the form of a Kaiser window with the Kaiser window parameter β = β (other values are of course also possible; typical values lie in the interval 1 < β < 10 . In Fig. 4(d) the shaping window is applied to the complex cepstrum to form the shaped complex cepstrum that is forwarded to step S4 in the method of Fig. 3.

Although the shaping window in this example has been applied to the complex cepstrum, it is appreciated that the window can equally well be applied directly to the real cepstrum coefficients before they are converted to complex cepstrum coefficients. Mathematically the result will be the same. Further- more, since taps M /2 + 1,...,M-I in the complex cepstrum are set to zero anyway to form a complex cepstrum that corresponds to a minimum-phase filter, actually only one half of the Kaiser window is required.

In the embodiment described above a Kaiser window has been used for shaping the sequence of coefficients A _1111n (n) . However, other smoothly decaying window functions, such as Hamming, Hanning, Blackman, Bartlett windows, are also feasible.

Fig. 5 is a diagram comparing the embodiment of Fig. 3 with the prior art method B in accordance with Fig. 2. In both cases a DFT of length M=64 was used to determine an estimated minimum-phase filter corresponding to the indicated desired magnitude response. As can be seen from Fig. 5 the present invention gives a stop band rejection in the order of 30 dB even with such a short DFT. The corresponding stop band rejection for the filter obtained by method B is only in the order of 10 dB. In fact, it has been found that in order to obtain a comparable stop band rejection with method B, a DFT having a length M equal to at least 512 is required. Thus, the significantly less complex

method of the present invention provides estimated minimum-phase filters with excellent stop band rejection.

From Shannon's sampling theorem it is known that replacement of the inte- grals in Fourier transforms by sums can be done without error if the spectrum of the signal is zero for all frequencies above half the sampling frequency. In the problem at hand, the logarithm of the system function H(_») = ln|H(ω) | must be sampled before computing the inverse Fourier transform. Hence, it should be ensured that a condition similar to that of the sampling theorem holds for ln|H(ω) | . Typically in echo cancellation and noise reduction applications, the dynamically updated desired magnitude response H(ω _k ) is computed in real time for a certain set of frequency bands based on some algorithm. The desired filter may need to pass some of the frequency bands and suppress other bands, so that the desired magnitude response is piecewise constant. A piecewise constant function is, however, not limited to contain only frequency components below a certain limit. Such a function is preferably smoothed before sampling to avoid aliasing. It should also be noted that even if the desired frequency response fulfils the sampling theorem, there is a possibility that computing the logarithm may introduce some components that do not fulfil the theorem. Hence, according to a further embodiment of the present invention the desired frequency response (or its logarithm) is smoothed before computing the cepstrum. The smoothing is accomplished by low pass filtering H(ω _k ) . This can be done for example by applying a frequency domain Hanning window to H{ω _k ) .

Fig. 6 is a flow chart illustrating an embodiment of the method in accordance with the present invention that includes smoothing in accordance with the previous paragraph. This method includes the same steps as the method of Fig. 3, but is supplemented by a smoothing step SlA between steps Sl and S2.

Fig. 7 is a diagram comparing the embodiment of Fig. 3 with the embodiment of Fig. 6. In this example, which is based on the same DFT length (M=64), the original magnitude response of Fig. 5 has been convolved (smoothed) with a Hanning window having the symmetrical kernel (1/4,1/2,1/4) , i.e.

H s.moothed (ω _k ) = ±H(ω _kA ) + ±H(ω _k ) + ±H(ω _M ) (26)

As can be seen in Fig, 7 the smoothing step improves the stop band rejection by about a further 15 dB compared to the example of Fig. 5.

From Fig. 7 it is also noted that the smoothing step widens the passband of the original magnitude response. This in turn results in a widened pass band for the minimum-phase filter obtained. If the passband width is of importance, widening of the passband can be counteracted by narrowing the passband of the original magnitude response one or several bins before smoothing.

In the description above a Hanning window has been described (due to its simple kernel) for smoothing the magnitude response. However, other smoothing windows, such as Hamming, Blackman, Bartlett, Kaiser windows, are also feasible.

Fig. 8 is a block diagram illustrating an embodiment of an apparatus for reducing digital filter delay in accordance with the present invention illustrated in an echo canceller environment. An estimator 10 estimates the desired filter magnitude, for example in accordance with the method described in [4]. Here

φ _x (ω _k ) represent an estimate of the power spectral density of a signal x(n) contaminated by echo and φ _e (co _k ) represents an estimate of the power spectral density of the echo signal. The determined magnitude response is passed through a smoothing filter 12 (optional, as indicated by the dashed lines), for example a Hanning filter. The (optionally) smoothed magnitude response is

forwarded to a logarithm block 14, which outputs the logarithm of its input signal. A frame containing M output signal samples from block 14 is forwarded to an inverse fast Fourier transformation block 26. The output signal from block 16 forms the real cepstrum, which in block 18 is transformed into a complex cepstrum associated with a minimum-phase filter. The complex cepstrum is shaped in a shaping window 20, for example a Kaiser window. The shaped complex cepstrum is forwarded to a fast Fourier transformation block 22. The output samples of block 22 are forwarded to an exponential function block 24. The final estimate of the minimum-phase filter from block 24 is stored in a filter block 26. This filter is used to suppress the echo contained in input signal x(n). Before signal x(n) is filtered it is typically transformed to the frequency domain by a fast Fourier transformation block 28. After filtering in filter 26, the output signal y(n) is obtained from an inverse fast Fourier transformation block 30.

As an alternative, the filtering of input signal x(n) may also be performed in the time domain, provided that the minimum-phase filter from block 24 is transformed into a corresponding impulse response in an IFFT block.

The different blocks of the echo canceller of Fig. 8 are typically implemented by a micro processor or a micro/ signal processor combination and corresponding software.

Although the present invention has been described in detail with reference to echo cancellation, the same principles may be used in other applications, for example noise reduction in telecommunication or audio systems. As is described in [4] this application may be based on a serially connected noise reduction filter with a desired magnitude response, the delay of which should be as small as possible. The essential difference is that instead of de- termining the desired magnitude response from an estimate φ _x (co _k ) of the power spectral density of the input signal and an estimate φ _e (co _k ) of the power spectral density of the echo signal, the desired magnitude response in

noise reduction is determined from an estimate of the power spectral density of the noisy input signal and an estimate of the power spectral density of the pure noise, as described in [4]. With these changes the block diagram of Fig. 8 may be used to implement a noise suppressor.

The present invention proposes a method for design of minimum-phase filters from a magnitude frequency response (with arbitrary phase) using only short DFTs. This results in a low complexity method, which outperforms all other known methods in terms of complexity and resulting stop band at- tenuation of the resulting minimum-phase filters. The method is particularly suitable for environments where a multitude of identical channels have to be processed (e.g. noise reduction, echo cancelling devices) at a minimum cost of computational operations and heat dissipation. The method allows for processing of more channels simultaneously than would have been possible using state of the art methods, given a limited amount of computational power. Alternatively, the reduced complexity provided by the present invention may be used to reduce heat dissipation and/ or improve battery utilization.

It will be understood by those skilled in the art that various modifications and changes may be made to the present invention without departure from the scope thereof, which is defined by the appended claims.

REFERENCES

[ 1] ITU-T G.131 , Talker echo and its control

[2] US 658107 Bl

[3] Kucera V., Discrete linear control, John Wiley & Sons, 1979

[4] WO 01/ 18960 Al

[5] A.V. Oppenheim and R.W. Schafer, Discrete-time signal processing, Prentice-Hall, Inter.Ed., 1989, Chapter 12

[6] Bysted, T.K., Aliasing in the complex cepstrum of linear-phase signals, International Conference on Information, Communication and Signal

Processing, 1997