CN111010144B - Improved two-channel IIR QMFB design method - Google Patents

Abstract

The invention discloses an improved two-channel IIR QMF group design method, which aims at the problem of minimum maximization of the phase of an all-pass filter and the phase of an overall distortion transfer function, and performs first-order Taylor expansion on an obtained nonlinear optimization target by determining a corresponding phase error and a proper weighted value, and converts the phase error and the proper weighted value into a linear problem so as to obtain the optimal coefficient of the all-pass filter. The biggest improvement of the invention is that only the performance of the all-pass filter or the analysis-synthesis filter is not considered any more, but the performance of the all-pass filter and the analysis-synthesis filter are considered comprehensively, so that each sub-filter through which a signal passes has a more approximately linear phase, the possibility of phase distortion is reduced, and the design expected target is achieved. Furthermore, the invention does not discuss the amplitude, but only the phase of the all-pass filter and the overall transfer function, and the amplitude is controlled entirely by the phase.

Description

Improved two-channel IIR QMFB design method

Technical Field

The invention belongs to the technical field of digital signal processing, and particularly relates to an improved two-channel IIR orthogonal mirror image filter bank design method based on an all-pass filter.

Background

Digital signal processing systems have been replacing analog systems in the last decades and today the digital signal processing system can be found in a variety of different devices, such as mobile communication systems, consumer and automotive electronics, or hearing aids. A decisive advantage of digital signal processing is that a system can be realized by digital signal processing and not or hardly by analog processing. An important component of many digital signal processing algorithms is the filter and its series connection as a filter bank. The filter bank is mainly used for spectral analysis of signals, as a cross-multiplexer, or processing sub-band frequency domain or time domain signals. The last-mentioned processing of subband signals requires an analysis-synthesis filter bank to obtain the reconstructed time-domain signal. The invention researches a Quadrature Mirror Filter Bank (QMFB).

Two-channel QMFB has been applied in more and more fields in recent years, such as subband coding of speech and image signals, wavelet-based design, and the like. Due to the wide application of QMFB, more and more attention is paid to its design.

X.zhang and h.iwakura proposed in 1994 a method of designing quadrature mirror filter banks using Digital Allpass Filters (DAFs). The QMFB, which consists of an all-pass filter and a pure delay section, has an approximately linear phase response, and the phase distortion can be minimized by using another additional DAFs as an equalizer. In the paper, a QMFB is designed by using a DAFs and a pure delay, and the method for designing the DAFs is an algorithm proposed in 1992, which can be simply described as a new method for designing the equal ripple phase response DAFs based on a eigenvalue problem, and obtains the optimal filter coefficient by calculating the maximum eigenvector and an iterative method. The algorithm obtains the all-pass filter by using an equal ripple FIR design method, but the calculation is extremely complex; and the algorithm does not take into account the phase error of the all-pass filter and the pure delay combination and the resulting stop-band error of the low-pass filter.

Ju-Hong Lee and Yi-Lin Shieh proposed in 2013 a method that takes into account the stopband error of a low-pass filter in the analysis filter. However, without considering the performance of the all-pass filter, the theory of dual-channel QMFB using two IIR DAFs yielded a Chebyshev approximation for the expected group delay response of the IIR DAFs and the magnitude response of the low-pass analysis filter. The paper takes group delay error and amplitude linear weighting of low-pass stop band as an objective function, which is a highly nonlinear problem, and the two-channel QFB obtained by solving an analytic solution by a frequency sampling and iterative approximation method and a weighted least square method has approximate linear phase response and has no amplitude distortion. However, the prototype filter coefficients of this method are the analytical solutions obtained directly by calculating an equation composed of an approximate matrix inversion and matrix multiplication, and the matrix inversion itself has a high computational complexity, and in addition, there may be no matrix inversion due to the problem of filter length.

Disclosure of Invention

Aiming at the defects of the prior art, the invention provides an improved QMFB design method of a two-channel IIR.

The invention discloses an improved two-channel IIR QMFB design method, which specifically comprises the following steps:

step one, according to design requirements, determining the number L of frequency points on a full frequency band and the order N of two all-pass filters₁And N₂Passband cut-off frequency omega of a low-pass filter of the analysis filter_pStopband cut-off frequency omega_sLet the iteration initial coefficient k be 0, the k-th all-pass filter coefficient a_i(k) 0, initial weight value W _i1, 2; wherein N is₁＝N₂+1；

Step two, determining the actual phase error and the whole distortion transfer function T (e) of the all-pass filter^jω) The actual phase error.

2.1. Determining an ideal phase of an all-pass filter

The ideal phase of the all-pass filter is at omega ∈ [0, pi ]]Satisfies the following conditions: when the phase of the IIR filter satisfies ω ═ 0, the phase is 0; when ω is pi, the phase satisfies θ (pi) — N pi; the all-pass filter is stable at this time; the ideal phases of all-pass filters are each theta_d1(ω)＝-N₁ω +0.25 ω and θ_d2(ω)＝-N₂ω-0.25ω。

2.2. Actual phase error of all-pass filter

Obtaining an expression of the all-pass filter, an actual phase expression of the all-pass filter and an actual phase error of the all-pass filter through expressions (1), (2) and (3), wherein omega ∈ [0,2 omega ]_p]，i＝1,2。

θ_ei(ω)＝θ_i(ω)-θ_di(ω) (3)

Wherein a is_i(n) represents a filter coefficient a_iN is 1,2, …, N_i；

2.3. Solving the overall distortion transfer function T (e)^jω) Actual phase error of theta_eT(ω)

The equations (4), (5) and (6) are the overall distortion transfer functions

Ideal phase, actual phase and actual phase error of where ω ∈ [0, ω_p]。

θ_dT(ω)＝θ_d1(2ω)+θ_d2(2ω)-ω (4)

θ_T(ω)＝θ₁(2ω)+θ₂(2ω)-ω (5)

θ_eT(ω)＝θ_T(ω)-θ_dT(ω) (6)

Wherein theta is₁(2ω)、θ₂(2ω)、θ_d1(2ω)、θ_d2(2. omega.) are represented by formula (7), formula (8), formula (9) and formula (10), respectively

θ_d1(2ω)＝-2N₁ω+0.5ω (9)

θ_d2(2ω)＝-2N₂ω+0.5ω (10)

Step three, solving to obtain the filter coefficient a of the all-pass filter under the kth iteration_i(k)。

3.1. Obtaining an objective optimization function representation as

min W₁||θ_e1(ω)||+W₂||θ_e2(ω)||+||θ_eT(ω)|| (11)

Theta in the formula (7)_e1(ω)、θ_e2(ω) and θ_eT(ω) in turn means an all-pass filter A₁(e^jω)、A₂(e^jω) And the overall distortion transfer function T (e)^jω) The actual phase error.

The objective optimization function is a maximum minimization problem and is also a nonlinear problem, and the step 3-2 can convert the nonlinear problem into a linear problem;

3.2. calculate theta from the k-1 iteration_ei(ω) with respect to the coefficient a_iFirst partial derivative of (k-1)

As shown in equation (12).

3.3. Computing the all-pass filter A at the kth iteration_i(e^jω) Phase error of (theta)_ei(ω) is represented by the formula (13).

In formula (13), △_iRepresenting the all-pass filter A at the kth iteration_i(e^jω) Increment of coefficient of, Δ_i＝a_i(k)-a_i(k-1)。

Calculating the integral distortion transfer function T (e) at the k iteration^jω) Actual phase error of theta_eTAnd (omega) is represented by formula (14).

3.4. The optimization objective function at the k-th iteration is expressed asSolving the convex optimization problem to determine coefficient increment △ of the kth iteration₁And △₂。

And step four, if the formula (16) does not hold, turning to step five. If equation (16) is satisfied, a is_i(k) As a final designed all-pass filter A_i(e^jω) And (4) ending the iteration.

In the formula (16), E^kIs the maximum value of the objective function in the k-th iteration, E^k-1Is the maximum of the objective function in the k-1 th iteration. μ is a set threshold.

E^k＝max(W₁|θ_e1(a₁(k),ω)|)+max(W₂|θ_e2(a₂(k),ω)|)+max(|θ_eT(a(k),ω)|) (17)

Step five, calculating the weighted value W according to the envelope of the group delay error_i

First, the group delay error of the k iteration is calculated

Wherein g is_diAll-pass filter A representing the kth iteration_i(e^jw) Actual group delay of τ_diRepresenting an all-pass filter A_i(e^jω) The desired group delay of; then calculate out

Envelope of

If not satisfied

Make it

Then order

k is k +1, and the step III is returned; otherwise directly order

And k is k +1, and the step three is returned. Wherein the threshold is set.

The invention has the beneficial effects that:

the invention does not only consider the performance of the all-pass filter or the performance of the analysis-synthesis filter, but comprehensively considers the performance of the all-pass filter and the analysis-synthesis filter, ensures that each sub-filter through which a signal passes has a more approximately linear phase, reduces the possibility of phase distortion, and achieves the design expected target. Compared with other methods, the method has the greatest improvement that the method does not discuss the amplitude, only discusses the phase of the all-pass filter and the whole transmission function, and the amplitude is completely controlled by the phase; this method ensures that the signal is free of phase distortion at any time period and a better performance index can be obtained in a limited number of iterations.

Drawings

FIG. 1 shows an all-pass filter A plotted for the coefficients in Table 1₁The phase error of (2).

FIG. 2 shows an all-pass filter A plotted with the coefficients of Table 1₂The phase error of (2).

Fig. 3 is a graph of the QMFB amplitude frequency response plotted against the coefficients in table 1.

Fig. 4 is a QMFB reconstruction error frequency response plot plotted against the coefficients in table 1.

Fig. 5 is a graph of the phase error of the reconstruction filter plotted against the coefficients in table 1.

Detailed Description

The invention is further described below with reference to the accompanying drawings.

The improved two-channel IIR orthogonal mirror image filter bank design method based on the all-pass filter comprises the following specific steps:

step one, according to design requirements, determining the number L of frequency points on a full frequency band and the order N of two all-pass filters₁And N₂Passband cut-off frequency omega of a low-pass filter of the analysis filter_pStopband cut-off frequency omega_sLet the iteration initial coefficient k be 0, the k-th all-pass filter coefficient a_i(k) 0, initial weight value W _i1, 2; wherein N is₁＝N₂+1；

Step two, determining the actual phase error and the whole distortion transfer function T (e) of the all-pass filter^jω) The actual phase error.

2.1. Determining an ideal phase of an all-pass filter

The ideal phase of the all-pass filter is at omega ∈ [0, pi ]]Satisfies the following conditions: when the phase of the IIR filter satisfies ω ═ 0, the phase is 0; when ω is pi, the phase satisfies θ (pi) — N pi; the all-pass filter is stable at this time; the ideal phases of all-pass filters are each theta_d1(ω)＝-N₁ω +0.25 ω and θ_d2(ω)＝-N₂ω-0.25ω。

2.2. Actual phase error of all-pass filter

Obtaining an expression of the all-pass filter, an actual phase expression of the all-pass filter and an actual phase error of the all-pass filter through expressions (1), (2) and (3), wherein omega ∈ [0,2 omega ]_p]，i＝1,2。

θ_ei(ω)＝θ_i(ω)-θ_di(ω) (3)

Wherein a is_i(n) represents a filter coefficient a_iN is 1,2, …, N_i；

2.3. Solving the overall distortion transfer function T (e)^jω) Actual phase error ofθ_eT(ω)

The equations (4), (5) and (6) are the overall distortion transfer functions

Ideal phase, actual phase and actual phase error of where ω ∈ [0, ω_p]。

θ_dT(ω)＝θ_d1(2ω)+θ_d2(2ω)-ω (4)

θ_T(ω)＝θ₁(2ω)+θ₂(2ω)-ω (5)

θ_eT(ω)＝θ_T(ω)-θ_dT(ω) (6)

Wherein theta is₁(2ω)、θ₂(2ω)、θ_d1(2ω)、θ_d2(2. omega.) are represented by formula (7), formula (8), formula (9) and formula (10), respectively

θ_d1(2ω)＝-2N₁ω+0.5ω (9)

θ_d2(2ω)＝-2N₂ω+0.5ω (10)

Step three, solving to obtain the filter coefficient a of the all-pass filter under the kth iteration_i(k)。

3.1. Obtaining an objective optimization function representation as

min W₁||θ_e1(ω)||+W₂||θ_e2(ω)||+||θ_eT(ω)|| (11)

Theta in the formula (7)_e1(ω)、θ_e2(ω) and θ_eT(ω) in turn means an all-pass filter A₁(e^jω)、A₂(e^jω) And the overall distortion transfer function T (e)^jω) The actual phase error.

The objective optimization function is a maximum minimization problem and is also a nonlinear problem, and the step 3-2 can convert the nonlinear problem into a linear problem;

3.2. calculate theta from the k-1 iteration_ei(ω) with respect to the coefficient a_iFirst partial derivative of (k-1)

As shown in equation (12).

3.3. Computing the all-pass filter A at the kth iteration_i(e^jω) Phase error of (theta)_ei(ω) is represented by the formula (13).

In formula (13), △_iRepresenting the all-pass filter A at the kth iteration_i(e^jω) Increment of coefficient of, Δ_i＝a_i(k)-a_i(k-1)。

Calculating the integral distortion transfer function T (e) at the k iteration^jω) Actual phase error of theta_eTAnd (omega) is represented by formula (14).

3.4. The optimization objective function representation is solved as a convex optimization problem in the k iteration, and coefficient increment △ of the k iteration is determined₁And △₂。

And step four, if the formula (16) does not hold, turning to step five. If equation (16) is satisfied, a is_i(k) As a final designed all-pass filter A_i(e^jω) And (4) ending the iteration.

In the formula (16), E^kIs the maximum value of the objective function in the k-th iteration, E^k-1Is the maximum of the objective function in the k-1 th iteration. μ is a set threshold.

E^k＝max(W₁|θ_e1(a₁(k),ω)|)+max(W₂|θ_e2(a₂(k),ω)|)+max(|θ_eT(a(k),ω)|) (17)

Step five, calculating the weighted value W according to the envelope of the group delay error_i

First, the group delay error of the k iteration is calculated

Wherein g is_diAll-pass filter A representing the kth iteration_i(e^jw) Actual group delay of τ_diRepresenting an all-pass filter A_i(e^jω) The desired group delay of; then calculate out

Envelope of

If not satisfied

Make it

Then order

k is k +1, and the step III is returned; otherwise directly order

And k is k +1, and the step three is returned. Wherein the threshold is set.

For the effectiveness of the invention, computer simulation simulations were performed on the invention.

Design requirements in simulation: order N of all-pass filter₁＝9，N₂The number of frequency points L in the full frequency band is 8N as 8₁+1, passband cut-off frequency ω_p0.4 pi, stop band cut-off frequency omega_s0.6 pi, minimizes reconstruction errors, and maximizes stop-band attenuation.

The filter coefficients of the final filter were obtained by 10 iterations with the design method of the present invention, as shown in table 1, the all-pass filter a is drawn by the coefficients in table 1₁The phase error of (2) is shown in FIG. 1, the coefficients of Table 1 are plotted for all-pass filter A₂Is shown in fig. 2. The corresponding QMFB amplitude-frequency response, reconstruction error frequency response and phase error response are shown in fig. 3, 4 and 5.

TABLE 1 all-pass filter A in QMFB designed by the method of the present invention₁And A₂Coefficient table (2)

And finally, calculating the maximum stop band attenuation (PSR), the maximum group delay error (MVPGD) of the passband, the maximum phase error (MVPR) of the QMFB reconstruction response, the maximum group delay error (MVGR) and the reconstruction response error (MVFBR) of the low-pass analysis filter by using the obtained filter coefficients. The calculation formula is as follows:

wherein H₀(e^jω) Analyzing the amplitude response of the filter for the obtained low-pass; gd is the ideal group delay of QMFB; t (e)^jω)

Is the overall distortion transfer function. The calculated indices are shown in table 2.

The QMFB obtained by the present invention is shown in Table 2 together with the QMFB index pairs obtained by the Lee and Shieh design methods.

TABLE 2 comparison of key indices of the present invention with Lee and Shieh methods

As can be seen from Table 2, the QMFB obtained by the design method of the present invention is significantly superior to the design methods of Lee and Shieh in all the above key indexes.

Claims (1)

1. The improved two-channel IIR QMFB design method is characterized by comprising the following steps:

step one, according to design requirements, determining the number L of frequency points on a full frequency band and the order N of two all-pass filters₁And N₂Passband cut-off frequency omega of a low-pass filter of the analysis filter_pStopband cut-off frequency omega_sLet the iteration initial coefficient k be 0, the k-th all-pass filter coefficient a_i(k) 0, initial weight value W_i1, 2; wherein N is₁＝N₂+1；

Step two, determining the actual phase error and the whole distortion transfer function T (e) of the all-pass filter^jω) Actual phase error of (2);

2.1. determining an ideal phase of an all-pass filter

The ideal phase of the all-pass filter is at omega ∈ [0, pi ]]Satisfies the following conditions: when the phase of the IIR filter satisfies ω ═ 0, the phase is 0; when ω is pi, the phase satisfies θ (pi) — N pi; the all-pass filter is stable at this time; all-throughThe ideal phases of the filters are respectively theta_d1(ω)＝-N₁ω +0.25 ω and θ_d2(ω)＝-N₂ω-0.25ω；

2.2. Actual phase error of all-pass filter

Obtaining an expression of the all-pass filter, an actual phase expression of the all-pass filter and an actual phase error of the all-pass filter through expressions (1), (2) and (3), wherein omega ∈ [0,2 omega ]_p]，i＝1,2；

θ_ei(ω)＝θ_i(ω)-θ_di(ω) (3)

Wherein a is_i(n) represents a filter coefficient a_iN is 1,2, …, N_i；

2.3. Solving the overall distortion transfer function T (e)^jω) Actual phase error of theta_eT(ω)

The equations (4), (5) and (6) are the overall distortion transfer functions

Ideal phase, actual phase and actual phase error of where ω ∈ [0, ω_p]；

θ_dT(ω)＝θ_d1(2ω)+θ_d2(2ω)-ω (4)

θ_T(ω)＝θ₁(2ω)+θ₂(2ω)-ω (5)

θ_eT(ω)＝θ_T(ω)-θ_dT(ω) (6)

Wherein theta is₁(2ω)、θ₂(2ω)、θ_d1(2ω)、θ_d2(2. omega.) are represented by formula (7), formula (8), formula (9) and formula (10), respectively

θ_d1(2ω)＝-2N₁ω+0.5ω (9)

θ_d2(2ω)＝-2N₂ω+0.5ω (10)

Step three, solving to obtain the filter coefficient a of the all-pass filter under the kth iteration_i(k)；

3.1. Obtaining an objective optimization function representation as

min W₁||θ_e1(ω)||+W₂||θ_e2(ω)||+||θ_eT(ω)|| (11)

Theta in formula (11)_e1(ω)、θ_e2(ω) and θ_eT(ω) in turn means an all-pass filter A₁(e^jω)、A₂(e^jω) And the overall distortion transfer function T (e)^jω) Actual phase error of (2);

the objective optimization function is a maximum minimization problem and is also a nonlinear problem, and step 3.2 can convert the nonlinear problem into a linear problem;

3.2. calculate theta from the k-1 iteration_ei(ω) with respect to the coefficient a_iFirst partial derivative of (k-1)

As shown in formula (12);

3.3. computing the all-pass filter A at the kth iteration_i(e^jω) Phase error of (theta)_ei(omega) is represented by formula (13);

in formula (13), △_iRepresenting the all-pass filter A at the kth iteration_i(e^jω) Increment of coefficient of, Δ_i＝a_i(k)-a_i(k-1)；

Calculating the integral distortion transfer function T (e) at the k iteration^jω) Actual phase error of theta_eT(ω) is represented by the formula (14);

3.4. the optimization objective function representation is solved as a convex optimization problem in the k iteration, and coefficient increment △ of the k iteration is determined₁And △₂；

Step four, if the formula (16) does not hold, turning to step five; if equation (16) is satisfied, a is_i(k) As a final designed all-pass filter A_i(e^jω) The iteration is ended;

in the formula (16), E^kIs the maximum value of the objective function in the k-th iteration, E^k-1Is the maximum of the objective function in the k-1 iteration; mu is a set threshold value;

E^k＝max(W₁|θ_e1(a₁(k),ω)|)+max(W₂|θ_e2(a₂(k),ω)|)+max(|θ_eT(a(k),ω)|) (17)

step five, calculating the weighted value W according to the envelope of the group delay error_i

First, the group delay error of the k iteration is calculated

Wherein g is_diAll-pass filter A representing the kth iteration_i(e^jw) Actual group delay of τ_diRepresenting an all-pass filter A_i(e^jω) The desired group delay of; then calculate out

Envelope of

If not satisfied

Make it

Then order

k is k +1, and the step III is returned; otherwise directly order

k is k +1, and the step III is returned; wherein the threshold is set.

Images