CN105719255B - Bandwidth varying linear-phase filter method based on Laplace structure - Google Patents

Abstract

The present invention discloses a filter design method suitable for Laplacian pyramid structure, which mainly solves the problem that existing filters with Laplacian pyramid structure cannot have linear phase at the same time except Haar filter. and orthogonality problems. Its technical scheme is: first set the filter length N, the passband cut-off frequency ω _p and the stopband start frequency ω _s , and the sampling factor M; then according to these parameters, call the firpm function in MatLab to generate the prototype filter p; Half of the coefficients of the prototype filter are optimized by calling the fmincon function, and the synthesis filter is obtained according to the optimized result; finally, the decomposition filter can be obtained according to the time-domain inversion relationship. By providing a reasonable orthogonal constraint condition, the invention not only enables the filter to have orthogonal and linear phase characteristics, but also has variable bandwidth, and provides wider application conditions for the filter of the Laplacian pyramid structure.

Description

Variable Bandwidth Linear Phase Filter Method Based on Laplace Structure

technical field

The invention belongs to the technical field of signal processing, and in particular relates to a design method of a linear phase filter with variable bandwidth, which can be used for image fusion, addition, compression, denoising and edge detection.

Background technique

The Laplacian pyramid structure LP was first proposed by Burt et al., and was originally used for image coding and compression. Later, people extended the Laplacian pyramid structure LP into a multi-scale and multi-resolution analysis tool, which is applied to many fields of image processing, such as image fusion, image denoising, etc.

Figure 1 shows the structure diagram of Laplacian pyramid. In Figure 1, for a given input signal x, the output c is an approximate approximation of x, and d is the residual between the original signal and its predicted signal p. In this structure, the analysis filter h and the synthesis filter g are usually a pair of low-pass filters, and the analysis matrix H and the synthesis matrix G constitute a set of transformation pairs.

In the traditional Laplacian pyramid structure LP design, both the decomposition filter h and the synthesis filter g are low-pass channels of the multi-channel complete reconstruction filter bank. The filter bank structure is shown in Figure 2 Show. In each channel, the filter and its up-down sampling factor form a forward and reverse transformation pair. Therefore, the filter of the low-pass channel in the filter bank is often directly applied to the Laplacian pyramid structure LP. Although this method provides help for the construction of the Laplacian pyramid structure LP, it also limits the performance of the filter in the Laplacian pyramid structure LP.

The Laplacian pyramid structure LP plays an important role in many fields of image processing, but there is still a key problem in its design, that is: the filters in the existing Laplacian pyramid structure LP are all taken from multiple channels fully reconstructs the low-pass channels in the filterbank. This makes the filter in the Laplacian pyramid structure LP subject to the conditions of fully reconstructing the filter bank itself, and some properties are restricted, including the very important linear phase property in image processing. In the design of a fully reconstructed filter bank, although the linear phase characteristic is easy to achieve in a two-channel filter bank, it is often difficult to achieve in a multi-channel fully reconstructed filter bank, resulting in a bandwidth of Except for the dual-channel filter bank, it is difficult for the fully reconstructed filter bank of other bandwidths to have linear phase characteristics. Therefore, in the Laplacian pyramid structure LP, this indirect filter design method taken from the fully reconstructed filter bank will cause the filter in the Laplacian pyramid structure LP to be fully reconstructed filtering It is difficult to achieve variable bandwidth and lack of linear phase at the same time due to the limitation of the device group.

Contents of the invention

The purpose of the present invention is to propose a kind of variable bandwidth linear phase filter method based on the Laplace structure for the above-mentioned deficiencies in the prior art, to get rid of the design constraint of completely reconfiguring the filter bank, so that the frequency band of the filter is variable, And it has linear phase property while having quadrature property.

Technical scheme of the present invention is realized like this:

1. Technical principles

The decomposition part of the Laplacian pyramid structure is shown in Figure 1. When H and G satisfy the H·G=I relationship, this channel in the LP structure constitutes a completely reconstructed subsystem, that is, it can realize the pair Complete reconstruction of the signal in space. According to the relationship of H·G=I, the conditions that the analysis filter h and the synthesis filter g need to meet can be deduced, which can be expressed as follows,

Among them, N is the length of the filter, M is the sampling factor, and the bandwidth of the filter is m is the number of shifts, the value is any integer, δ(m) is a Dirac sequence,

The above conditions enable a variable bandwidth filter to have both quadrature and linear phase characteristics.

2. Technical solution

According to above-mentioned principle, technical scheme of the present invention is as follows:

Technical solution 1: A low-pass filter design method based on a variable bandwidth linear phase of Laplace structure, including the design of the low-pass synthesis filter g _L and the design of the low-pass decomposition filter h _L :

(1) Set the length N _L of the low-pass synthesis filter g _L , and the passband cut-off frequency stop band start frequency Bandwidth is Where M is the sampling factor, M and N _L are both integers, M≥2, N _L is an integer multiple of M,

(2) According to the above parameters, call the firpm function in MatLab to generate the impulse response sequence p _L (n) of the low-pass prototype filter p _L of linear phase, n=0,1,2...N _L -1;

(3) Determine half of the coefficient q _L (k) of the low-pass prototype filter p _L :

When N _L is an odd number,

When N _L is even,

(4) Take q _L (k) as the initial value of the optimization function fmincon, and optimize according to the following optimization formula:

st<g _L0 (n), g _L0 (n-mM)>＝δ(m)

Among them, g _L0 (n) is the low-pass filter impulse response free variable, when N _L is an even number, g _L0 (n) = [q _L (k), fliplr (q _L (k))], when N _L When it is an odd number, g _L0 (n)=[q _L (k), fliplr(q _L (k-1)], fliplr is the sequence flip function in MatLab, G _L (e ^jw ) is the function of g _L0 (n) Frequency response, α is the weight, the value is 0<α<1, m is the number of shifts, the value is any integer, δ(m) is the Dirac sequence,

In the case of satisfying the constraint conditions, adjust the weight α, and when the value of the objective function φ _L reaches the minimum, obtain half of the coefficient Q _L (k) of the optimized low-pass reconstruction filter g _L ;

(5) According to the half coefficient Q _L (k) of the low-pass synthesis filter after optimization, obtain the impulse response sequence g _L (n) of the low-pass synthesis filter g _L :

When N _L is an odd number, g _L (n)=[Q _L (k), fliplr(Q _L (k-1))];

When N _L is an even number, g _L (n)=[Q _L (k), fliplr(Q _L (k))];

(6) According to the impulse response sequence g _L (n) of the low-pass synthesis filter and the impulse response sequence h _L (n) of the low-pass decomposition filter satisfy the time-domain inversion relation: h _L (n)=g _L (N _L -1-n), the impulse response sequence h _L (n) of the low-pass decomposition filter can be obtained from the impulse response sequence g _L (n) of the low-pass synthesis filter.

Technical solution 2: A design method of a band-pass filter based on a Laplace structure variable bandwidth linear phase, including the design of the band-pass synthesis filter g _B and the design of the band-pass decomposition filter h _B :

1) Set the length N _B of the band-pass synthesis filter g _B , and the pass-band cut-off frequencies are respectively The start frequency of the stop band is The sampling factor is M and the bandwidth is Where M and N _B are both integers, M≥2, N _B is an integer multiple of M,

2) According to the above parameters, call the firpm function in MatLab to generate the impulse response sequence p _B (n) of the bandpass prototype filter p _B of linear phase, n=0,1,2...N _B -1;

3) Determine half of the coefficients q _B (k) of the bandpass prototype filter p _B :

When N _B is an odd number,

When N _B is even,

4) Use q _B (k) as the initial value of the optimization function fmincon, and optimize according to the following optimization formula:

st<g _B0 (n), g _B0 (n-mM)>＝δ(m)

Among them, g _B0 (n) is the free variable of the impulse response of the band-pass filter. When N _B is an even number, g _B0 (n)=[q _B (k), fliplr(q _B (k))], when N _B When it is an odd number, g _B0 (n)=[q _B (k), fliplr(q _B (k-1)], fliplr is the sequence flip function in MatLab, G _B (e ^jw ) is the frequency response of g _B0 , α is the weight, the value is 0<α<1, m is the number of shifts, the value is any integer, δ(m) is the Dirac sequence,

In the case of satisfying the constraint conditions, adjust the weight α, when the value of the objective function φ _B reaches the minimum, get half of the coefficient Q _B (k) of the optimized low-pass synthesis filter;

5) According to the half coefficient Q _B (k) of the band-pass synthesis filter after optimization, obtain the impulse response sequence g _B (n) of the band-pass synthesis filter g _B :

When N _B is an odd number, g _B (n)=[Q _B (k), fliplr(Q _B (k-1))];

When N _B is an even number, g _B (n)=[Q _B (k), fliplr(Q _B (k))];

6) According to the impulse response sequence g _B (n) of the band-pass synthesis filter and the impulse response sequence h _B (n) of the band-pass decomposition filter satisfy the time-domain inversion relation: h _B (n)=g _B (N _B -1-n), the impulse response sequence h _B (n) of the band-pass analysis filter can be obtained from the impulse response sequence g _B (n) of the band-pass synthesis filter.

Technical solution 3: A high-pass filter design method based on Laplace structure variable bandwidth linear phase, including the design of the high-pass synthesis filter g _H and the design of the high-pass decomposition filter h _H :

(A) Set the length N _H of the high-pass synthesis filter g _H , the passband cut-off frequency is The stopband start frequency is Bandwidth is Where M is the sampling factor, M and N _H are both integers, M≥2, N _H is an integer multiple of M,

(B) According to the above parameters, call the firpm function in MatLab to generate the impulse response sequence p _H (n) of the high-pass prototype filter p _H of linear phase, n=0,1,2...N _H -1;

(C) Determine half of the coefficients q _H (k) of the high-pass prototype filter p _H :

When N _H is an odd number,

When N _H is an even number,

(D) Use q _H as the initial value of the optimization function fmincon, and optimize according to the following optimization formula:

st<g _H0 (n), g _H0 (n-mM)>＝δ(m)

Among them, g _H0 (n) is the high-pass filter impulse response free variable: when N _H is an even number, g _H0 (n)=[q _H (k),-fliplr(q _H (k))], when N _H When it is odd, g _H0 (n)=[q _H (k), fliplr(q _H (k-1)], fliplr is the sequence flip function in MatLab, G _H (e ^jw ) is the function of g _H0 (n) Frequency response, α is the weight, the value is 0<α<1, m is the number of shifts, the value is any integer, δ(m) is the Dirac sequence,

In the case of satisfying the constraint conditions, adjust the weight α, when the value of the objective function φ _H reaches the minimum, get half of the coefficient Q _H (k) of the optimized high-pass synthesis filter;

(E) According to half of the coefficients Q _H (k) of the optimized high-pass synthesis filter, the impulse response sequence g _H (n) of the high-pass synthesis filter g _H is obtained:

When N _H is an odd number, g _H (n)=[Q _H (k), fliplr(Q _H (k-1))],

When N _H is an even number, g _H (n)=[Q _H (k),-fliplr(Q _H (k))];

(F) According to the impulse response sequence g _H (n) of the high-pass synthesis filter and the impulse response sequence h _H (n) of the high-pass decomposition filter satisfy the time-domain inversion relation: h _H (n)=g _H (N _H - 1-n), the impulse response sequence h _H (n) of the high-pass decomposition filter can be obtained from the impulse response sequence g _H (n) of the high-pass synthesis filter.

Compared with the prior art, the present invention has the following advantages:

First, the filter designed by the present invention can have linear phase and quadrature characteristics at the same time.

Existing Daubechies filters, except Haar filters, only have quadrature characteristics, but not linear phase characteristics; Biorthogonal filters have linear phase characteristics, but do not have quadrature characteristics. All other filters designed by the fully reconstructed filter bank, except the Haar filter, have only one of the quadrature characteristic and the linear phase characteristic, but cannot have these two characteristics at the same time. The invention gets rid of the design constraints of the complete reconstruction filter bank, and directly designs the filter in the Laplacian pyramid structure LP. When the decomposition filter h and the reconstruction filter g meet the conditions: When , the filter can have both linear phase and quadrature characteristics;

Second, the bandwidth of the present invention is variable.

The existing Haar filter bandwidth is The filter bandwidth of the present invention can be That is, the spectral range is

Third, the design process of the present invention is simple.

The existing Haar filter is to design both the synthesis filter g and the analysis filter h, but the filter of the present invention only needs to design the synthesis filter g, and then according to the impulse response sequence h of the analysis filter h The relationship h(n)=g(N-1-n) between (n) and the impulse response sequence g(n) of the synthesis filter g, then the decomposition filter h can be obtained, which greatly reduces the complexity of the design.

Description of drawings

Fig. 1 is the structural diagram of the decomposition part of the existing Laplace structure LP;

Fig. 2 is the structural diagram of existing multi-channel complete reconstruction filter bank;

Fig. 3 is a design flowchart of the present invention;

Fig. 4 is that the bandwidth of the present invention based on the pull structure is Performance simulation diagram of a low-pass filter with a length of 24;

Fig. 5 is that the bandwidth of the present invention based on the pull structure is Performance simulation diagram of a bandpass filter with a length of 60;

Fig. 6 is that the bandwidth of the present invention based on the pull structure is Performance simulation diagram of a high-pass filter with a length of 60.

Detailed ways

The technical solutions and effects of the present invention will be further described in detail below in conjunction with the accompanying drawings.

The filter designed in the present invention is based on the Laplacian pyramid structure, as shown in FIG. 1 , the filter includes a decomposition filter and a synthesis filter. After the input signal x passes through the decomposition filter h, the up-sampling factor M and the synthesis filter g in sequence, the approximation signal p of the input signal x is obtained, and the residual between the input signal and the approximation signal p is d.

Example 1: Designing a bandwidth based on the Laplace structure is A low-pass filter of length 24.

This example includes two parts: the design of the low-pass synthesis filter g _L and the design of the low-pass analysis filter h _L. 1. Design a low-pass reconstruction filter g _L

Step 1, setting the parameters of the low-pass synthesis filter g _L .

Set length N _L =24, sampling factor M=4, passband cut-off frequency stop band start frequency Where r _L is the adjustment parameter of the low-pass transition band. By changing the size of r _L , the width of the transition band can be adjusted conveniently. In this example, r _L =0.4.

Step 2, determine the low-pass prototype filter p _L .

The low-pass prototype filter p _L is essentially a low-pass impulse response sequence p _L (n). In this example, the impulse response sequence p _L (n) of the linear-phase low-pass prototype filter p _L is obtained by calling the firpm function in MatLab. ,n=0,1,... N _L -1.

The firpm function is a function for designing a finite-length impulse response sequence FIR filter, which is introduced in the book "Digital Signal Processing", and the filters designed by the firpm function all have linear phase characteristics.

Step 3, determine the impulse response sequence q _L (k) of half of the low-pass prototype filter p _L.

Since the low-pass prototype filter p _L has a linear phase characteristic, that is, its impulse response sequence p _L (n) is centrosymmetric, in order to reduce the computational complexity, only half of the impulse response sequence of p _L (n) is optimized The impulse response sequence g _L (n) of the low-pass synthesis filter g _L can be obtained.

Take half of the impulse response sequence q _L (k) of the impulse response sequence p _L (n) of the low-pass prototype filter, where when N _L is an even number, When N _L is odd, For example, if N _L =4, p _L (n)=[p _L (0)p _L (1) p _L (2) p _L (3)], the impulse response sequence of half of p _L (n) is : In this example, N _L = 24, so half of the impulse response sequence p _L (n) of the low-pass prototype filter

Step 4, setting parameters of the optimization function fmincon.

Suppose the initial value of the optimization function fmincon is q _L (k), the constraint condition is <g _L0 (n), g _L0 (n-mM)>=δ(m), and the objective function is Among them: g _L0 (n) is the free variable of the low-pass filter impulse response, m is the number of shifts, the value is any integer, g _L0 (n-mM) is the m-time shift sequence of g _L0 (n), α is the weight, 0<α<1, G _L (e ^jω ) is the frequency response of g _L0 (n), δ(m) is the Dirac sequence,

The optimization function fmincon is a commonly used optimization function in MatLab;

The value of g _L0 (n) depends on the parity of N _L :

When N _L is an even number, g _L0 (n)=[q _L (k), fliplr(q _L (k))],

When N _L is an odd number, g _L0 (n)=[q _L (k), fliplr(q _L (k-1))],

fliplr is the sequence flip function in MatLab. For example, if q _L (k)=[q _L (0) q _L (1) q _L (2) q _L (3)], N _L =4, the sequence q _L ( k) Flip left and right, which can be expressed as fliplr(q _L (k))=[q _L (3) q _L (2) q _L (1) q _L (0)];

The constraint condition <g _L0 (n), g _L0 (n-mM)>=δ(m) indicates that g _L0 (n) is orthogonal to g _L0 (n-mM) sequence, namely

The objective function is to optimize the passband and stopband of the filter at the same time. The weight α determines the degree of restriction on the passband and stopband of the filter during optimization. The larger the α, the greater the restriction on the passband, that is, the smoother the passband. The greater the stop band fluctuation.

Step 5, half of the impulse response q _L (k) of the impulse response of the low-pass prototype filter determines half of the coefficients Q _L (k) of the low-pass synthesis filter g _L.

For half of the impulse response q _L (k) of the impulse response p _L (n) of the low-pass prototype filter, optimize the optimization function fmincon to adjust the value of the weight α. In this example, α = 0.8. When the constraints are satisfied Next, when the objective function φ _L is the smallest, the function return value is half of the coefficient Q _L (k) of the low-pass synthesis filter g _L.

Step 6, determine the impulse response sequence g _L (n) of the low-pass synthesis filter g _L.

Since the low-pass synthesis filter g _L is a linear phase filter, its impulse response sequence g _L (n) is centrosymmetric, so the low-pass synthesis filter can be determined by half the coefficient Q _L (k) of the low-pass synthesis filter g _L The impulse response sequence g _L (n) of the device g _L :

When N _L is an even number, g _L (n)=[Q _L (k), fliplr(Q _L (k))],

When N _L is an odd number, g _L (n)=[Q _L (k), fliplr(Q _L (k-1))].

2. Design the low-pass decomposition filter h _L

Step 7, determine the impulse response sequence h _L (n) of the low-pass decomposition filter.

Set the impulse response g _L (n) of the low-pass synthesis filter and the impulse response sequence h _L (n) of the low-pass decomposition filter to satisfy the time-domain inversion relationship: h _L (n) = g _L (N _L -1 -n), according to the time-domain inversion relationship, the impulse response sequence h _L (n) of the low-pass analysis filter can be obtained, that is, the low-pass analysis filter h _L .

For example, if g _L (n)=[g _L (0) g _L (1) g _L (2) g _L (3)], N _L =4, g _L (n) and h _L (n) satisfy the time domain Inverting the relation, h _L (n)=g _L (N _L -n-1)=[g _L (3) g _L (2) g _L (1) g _L (0)].

From the above steps, the bandwidth based on the Laplace structure can be obtained as A low-pass filter of length 24 with a spectral support range of The performance of this low-pass filter is shown in Figure 4. in:

Figure 4(a) is the time-domain impulse response of the low-pass synthesis filter g _L ,

Figure 4(b) is the frequency domain response of the low-pass synthesis filter g _L ,

Figure 4(c) is the time-domain impulse response of the low-pass decomposition filter _hL ,

Figure 4(d) is the frequency domain response of the low-pass analysis filter _hL .

From Fig. 4 (a) and 4 (c), it can be seen that this low-pass filter belongs to a linear phase filter, which verifies that the impulse response sequence g _L (n) of the low-pass synthesis filter designed in this example satisfies < g _L (n), g _L (n-mM)>=δ(m), so the low-pass synthesis filter has orthogonality, and because the impulse response g _L (n) of the low-pass synthesis filter and the low-pass analysis filter The impulse response sequence h _L (n) of the filter satisfies the time-domain inversion relation: h _L (n) = g _L (N _L -1-n), so the impulse response h _L (n) of the low-pass decomposition filter also satisfies <h _L (n), h _L (n-mM)>=δ(m), that is, the low-pass analysis filter has orthogonality. From Figure 4(b) and 4(d), it can be seen that the spectrum support range of the low-pass filter is That is, the bandwidth is

In Example 2, the bandwidth based on the Laplace structure is A bandpass filter of length 60.

This example includes two parts: the design of the band-pass synthesis filter g _B and the design of the band-pass analysis filter h _B. 1. Design a bandpass synthesis filter g _B

Step 1, setting the parameters of the band-pass synthesis filter g _B.

Let the length N _B =60, the sampling factor M=3, and the cut-off frequency of the passband be The start frequency of the stop band is Where r _B is the adjustment parameter of the band passing through the transition band. By changing the size of r _B , the width of the transition band can be adjusted conveniently. In this case, r _B =0.3.

Step 2, determine the prototype bandpass filter p _B .

The band-pass prototype filter p _B is essentially a band-pass impulse response sequence p _B (n). In this example, the impulse response sequence p _B (n) of the linear-phase band-pass prototype filter p _B is obtained by calling the firpm function in MatLab. ,n=0,1,...N _B -1.

Step three, determine the impulse response sequence q _B (k) of half of the bandpass prototype filter p _B.

Since the band-pass prototype filter p _B has a linear phase characteristic, that is, its impulse response sequence p _B (n) is symmetrical in center, in order to reduce the computational complexity, it is only necessary to optimize the impulse response sequence of half of p _B (n) The impulse response sequence g _B (n) of the bandpass synthesis filter g _B can be obtained.

Take half of the impulse response sequence q _B (k) of the impulse response sequence p _B (n) of the bandpass prototype filter, where when N _B is an even number, When N _B is odd, For example, if _NB = 4, p _B (n) = [p _B (0) p _B (1) p _B (2) p _B (3)], the impulse response sequence of half of p _B (n) is : In this example, N _B = 60, so half of the impulse response sequence p _B (n) of the bandpass prototype filter

Step 4, setting parameters of the optimization function fmincon.

Assuming that the initial value of the optimization function fmincon is q _B (k), the constraint condition is <g _B0 (n), g _B0 (n-mM)>=δ(m), and the objective function is Among them: g _B0 (n) is the free variable of the impulse response of the bandpass filter, m is the number of shifts, the value is any integer, g _B0 (n-mM) is the m shift sequence of g _B0 (n), α is the weight, 0<α<1, G _B (e ^jω ) is the frequency response of g _B0 (n), δ(m) is the Dirac sequence,

The value of g _B0 (n) depends on the parity of N _B :

When N _B is an even number, g _B0 (n)=[q _B (k), fliplr(q _B (k))],

When N _B is an odd number, g _B0 (n)=[q _B (k), fliplr(q _B (k-1))],

fliplr is the sequence flip function in MatLab. The constraint condition <g _B0 (n), g _B0 (n-mM)>=δ(m) indicates that g _B0 (n) and g _B0 (n-mM) sequence are orthogonal, that is

The objective function is to optimize the passband and stopband of the filter at the same time. The weight α determines the degree of restriction on the passband and stopband of the filter during optimization. The larger the α, the greater the restriction on the passband, that is, the smoother the passband. The greater the stop band fluctuation.

Step five, half of the impulse response q _B (k) of the impulse response of the band-pass prototype filter determines half of the coefficient Q _L (k) of the low-pass synthesis filter g _B.

By optimizing the optimization function fmincon, adjust the value of the weight α for half of the impulse response q _{B (k) of the impulse response p B (} n) of the bandpass prototype filter. In this example, α = 0.7, in the case of satisfying the constraints Next, when the objective function φ _B is the smallest, the function return value is half of the coefficient Q _B (k) of the band-pass synthesis filter g _B.

Step six, determine the impulse response sequence g _B (n) of the band-pass synthesis filter g _B.

Since the band-pass synthesis filter g _B is a linear phase filter, its impulse response sequence g _B (n) is centrosymmetric, so the band-pass synthesis filter can be determined by half the coefficient Q _B (k) of the band-pass synthesis filter g _B The impulse response sequence g _B (n) of the device g _B :

When N _B is an even number, g _B (n)=[Q _B (k), fliplr(Q _B (k))],

When N _B is an odd number, g _L (n)=[Q _B (k), fliplr(Q _B (k-1))].

2. Design the bandpass analysis filter h _B

Step seven, determine the impulse response sequence h _B (n) of the bandpass analysis filter.

According to the impulse response g _B (n) of the band-pass synthesis filter and the impulse response sequence h _B (n) of the band-pass analysis filter satisfy the time-domain inversion relationship: h _B (n) = g _B (N _B - 1-n), the impulse response sequence h _B (n) of the analysis filter can be obtained, that is, the band-pass analysis filter h _B .

For example, if g _B (n)=[g _B (0) g _B (1) g _B (2) g _B (3)], N _B =4, g _B (n) and h _B (n) satisfy the time domain Inverting the relation, h _B (n)=g _B (N _B -n-1)=[g _B (3) g _B (2) g _B (1) g _B (0)].

From the above steps, the bandwidth based on the Laplace structure can be obtained as A bandpass filter of length 60 with a spectral support range of The performance of this bandpass filter is shown in Figure 5. in:

Figure 5(a) is the time-domain impulse response of the band-pass synthesis filter g _B ,

Figure 5(b) is the frequency domain response of the bandpass synthesis filter g _B ,

Figure 5(c) is the time-domain impulse response of the bandpass decomposition filter _hB ,

Figure 5(d) is the frequency domain response of the band-pass analysis filter _hB .

It can be seen from Figure 5(a) and 5(c) that the bandpass filter is a linear phase filter, which verifies that the impulse response sequence g _B (n) of the bandpass synthesis filter designed in this example satisfies <g _B (n), g _B (n-mM)>=δ(m), so the band-pass synthesis filter has orthogonality, and because the impulse response g _B (n) of the band-pass synthesis filter and the band-pass analysis filter The impulse response sequence h _B (n) of the filter satisfies the time-domain inversion relation: h _B (n) = g _B (N _B -1-n), so the impulse response h _B (n) of the band-pass decomposition filter also satisfies <h _B (n), h _B (n-mM)>=δ(m), that is, the bandpass analysis filter has orthogonality. From Figure 5(b) and 5(d), it can be seen that the spectrum support range of the bandpass filter is That is, the bandwidth is

Example 3. Design the bandwidth based on the Laplace structure as A high-pass filter of length 60.

This example includes two parts: the design of the high-pass synthesis filter g _H and the design of the high-pass decomposition filter h _H. 1), design high-pass synthesis filter g _H

Step A, setting the parameters of the high-pass synthesis filter g _H .

Let length N _H =60, sampling factor M=3, passband cut-off frequency stop band start frequency Where r _H is a high-pass transition band adjustment parameter, by changing the size of r _H , the width of the transition band can be adjusted conveniently, and r _H =0.2 is taken in this example.

Step B, determine the high-pass prototype filter _pH .

The high-pass prototype filter p _H is essentially a high-pass impulse response sequence p _H (n), this example obtains the impulse response sequence p _H (n) of the linear phase high-pass prototype filter p _H by calling the firpm function in MatLab, n= 0,1,...,N _H -1.

Step C, determine the impulse response sequence q _H (k) of half of the high-pass prototype filter p _H.

Since the high-pass prototype filter p _H has a linear phase characteristic, that is, its impulse response sequence p _H (n) is centrosymmetric, in order to reduce the computational complexity, it is only necessary to optimize the impulse response sequence of half of p _H (n). The impulse response sequence g _H (n) of the high-pass synthesis filter g _H can be obtained.

Take half of the impulse response sequence q _H (k) of the impulse response sequence p _H (n) of the high-pass prototype filter, where when N _H is an even number, When N _H is odd, For example, if N _H = 4, p _H (n) = [p _H (0) p _H (1) p _H (2) p _H (3)], the impulse response sequence of half of p _H (n) is : In this example, N _H = 24, so half of the impulse response sequence p _H (n) of the high-pass prototype filter

Step D, setting parameters of the optimization function fmincon.

Assuming that the initial value of the optimization function fmincon is q _H (k), the constraint condition is <g _H0 (n), g _H0 (n-mM)>=δ(m), and the objective function is Among them: g _H0 (n) is the free variable of the high-pass filter impulse response, m is the number of shifts, and the value is any integer, g _H0 (n-mM) is the m-time shift sequence of g _H0 (n), α is weight, 0<α<1, G _H (e ^jω ) is the frequency response of g _H0 (n), δ(m) is the Dirac sequence,

The value of g _H0 (n) depends on the parity of N _H :

When N _H is an even number, g _H0 (n)=[q _H (k),-fliplr(q _H (k))],

When N _H is an odd number, g _H0 (n)=[q _H (k), fliplr(q _H (k-1))],

fliplr is the sequence flip function in MatLab. For example, if q _H (k)=[q _H (0) q _H (1) q _H (2) q _H (3)], N _H =4, the sequence q _H ( k) Flip left and right, which can be expressed as fliplr(q _H (k))=[q _H (3) q _H (2) q _H (1) q _H (0)];

The constraint condition <g _H0 (n), g _H0 (n-mM)>=δ(m) indicates that g _H0 (n) is orthogonal to g _H0 (n-mM) sequence, namely

The objective function is to optimize the passband and stopband of the filter at the same time. The weight α determines the degree of restriction on the passband and stopband of the filter during optimization. The larger the α, the greater the restriction on the passband, that is, the smoother the passband. The greater the stop band fluctuation.

Step E, determine half of the coefficients Q _H (k) of the high-pass synthesis filter g _H according to half of the impulse response q _H (k) of the high-pass prototype filter.

For half of the impulse response q _H (k) of the impulse response p _H (n) of the high-pass prototype filter, optimize the optimization function fmincon to adjust the value of the weight α. In this example, α = 0.8. When the constraints are satisfied , when the objective function φ _H is the smallest, the function return value is half of the coefficient Q _H (k) of the high-pass synthesis filter g _H.

Step F, determine the impulse response sequence g _H (n) of the high-pass synthesis filter g _H.

Since the high-pass synthesis filter g _H is a linear phase filter, its impulse response sequence g _H (n) is centrosymmetric, so the high-pass synthesis filter g can be determined by the half coefficient Q _{H (} k) of the high-pass synthesis filter g H The impulse response sequence g _H (n) of _H :

When N _H is an even number, g _H (n)=[Q _H (k),-fliplr(Q _H (k))],

When N _H is an odd number, g _H (n)=[Q _H (k), fliplr(Q _H (k-1))].

Step G, obtain the high-pass decomposition filter h _H according to the impulse response g _H (n) of the high-pass synthesis filter.

Since the impulse response g _H (n) of the high-pass synthesis filter and the impulse response sequence h _H (n) of the high-pass decomposition filter satisfy the time-domain inversion relationship: h _H (n) = g _H (N _H -1- n), so the impulse response sequence h _H (n) of the high-pass analysis filter can be obtained through the time-domain inversion relationship, that is, the high-pass analysis filter h _H .

For example, if g _H (n)=[g _H (0) g _H (1) g _H (2) g _H (3)], N _H =4, g _H (n) and h _H (n) satisfy the time domain Inverting the relationship, h _H (n) = g _H (N _H -n-1) = [g _H (3) g _H (2) g _H (1) g _H (0)].

From the above steps, the bandwidth based on the Laplace structure can be obtained as A high-pass filter of length 60 with a spectral support range of The performance of this high-pass filter is shown in Figure 6. in:

Figure 6(a) is the time-domain impulse response of the high-pass synthesis filter g _H ,

Figure 6(b) is the frequency domain response of the high-pass synthesis filter g _H ,

Figure 6(c) is the time-domain impulse response of the high-pass decomposition filter h _H ,

Figure 6(d) is the frequency domain response of the high-pass decomposition filter h _H.

It can be seen from Figure 6(a) and 6(c) that the high-pass filter is a linear phase filter, which verifies that the impulse response sequence g _H (n) of the high-pass synthesis filter designed in this example satisfies < g _H (n ), g _H (n-mM)>=δ(m), so the high-pass synthesis filter has orthogonality, and the impulse response g _H (n) of the high-pass synthesis filter and the impulse response sequence h of the high-pass decomposition filter _H (n) satisfies the time-domain inversion relationship: h _H (n) = g _H (N _H -1-n), so the impulse response h _H (n) of the high-pass decomposition filter also satisfies <h _H (n), h _H (n-mM)>=δ(m), that is, the high-pass decomposition filter has orthogonality. From Figure 6(b) and 6(d), it can be seen that the spectrum support range of the high-pass filter is That is, the bandwidth is

Claims (3)

1. A low-pass filter design method based on the variable bandwidth linear phase of the Laplace structure, including the design of the low-pass synthesis filter g _L and the design of the low-pass decomposition filter h _L : (1) Set the length N _L of the low-pass synthesis filter g _L , and the passband cut-off frequency stop band start frequency Bandwidth is Where M is the sampling factor, M and N _L are both integers, M≥2, N _L is an integer multiple of M, (2) According to the above parameters, call the firpm function in MatLab to generate the impulse response sequence p _L (n) of the low-pass prototype filter p _L of linear phase, n=0,1,2...N _L -1; (3) Determine half of the coefficient q _L (k) of the low-pass prototype filter p _L : When N _L is an odd number, When N _L is even, (4) Take q _L (k) as the initial value of the optimization function fmincon, and optimize according to the following optimization formula: st<g _L0 (n), g _L0 (n-mM)>＝δ(m) Among them, g _L0 (n) is the low-pass filter impulse response free variable, when N _L is an even number, g _L0 (n) = [q _L (k), fliplr (q _L (k))], when N _L When it is an odd number, g _L0 (n)=[q _L (k), fliplr(q _L (k-1)], fliplr is the sequence flip function in MatLab, G _L (e ^jω ) is the function of g _L0 (n) Frequency response, α is the weight, the value is 0<α<1, m is the number of shifts, the value is any integer, δ(m) is the Dirac sequence, In the case of satisfying the constraint conditions, adjust the weight α, when the value of the objective function φ _L reaches the minimum, obtain half of the coefficient Q _L (k) of the optimized low-pass synthesis filter g _L ; (5) According to the half coefficient Q _L (k) of the low-pass synthesis filter after optimization, obtain the impulse response sequence g _L (n) of the low-pass synthesis filter g _L : When N _L is an odd number, g _L (n)=[Q _L (k), fliplr(Q _L (k-1))]; When N _L is an even number, g _L (n)=[Q _L (k), fliplr(Q _L (k))]; (6) According to the impulse response sequence g _L (n) of the low-pass synthesis filter and the impulse response sequence h _L (n) of the low-pass decomposition filter satisfy the time-domain inversion relation: h _L (n)=g _L (N _L -1-n), the impulse response sequence h _L (n) of the low-pass decomposition filter can be obtained from the impulse response sequence g _L (n) of the low-pass synthesis filter. 2. A bandpass filter design method based on the variable bandwidth linear phase of the Laplace structure, including the design of the bandpass synthesis filter g _B and the design of the bandpass decomposition filter h _B : 1) Set the length N _B of the band-pass synthesis filter g _B , and the pass-band cut-off frequencies are respectively The start frequency of the stop band is The sampling factor is M and the bandwidth is Where M and N _B are both integers, M≥2, N _B is an integer multiple of M, 2) According to the above parameters, call the firpm function in MatLab to generate the impulse response sequence p _B (n) of the bandpass prototype filter p _B of linear phase, n=0,1,2...N _B -1; 3) Determine half of the coefficients q _B (k) of the bandpass prototype filter p _B : When N _B is an odd number, When N _B is even, 4) Use q _B (k) as the initial value of the optimization function fmincon, and optimize according to the following optimization formula: st<g _B0 (n), g _B0 (n-mM)>＝δ(m) Among them, g _B0 (n) is the free variable of the impulse response of the band-pass filter. When N _B is an even number, g _B0 (n)=[q _B (k), fliplr(q _B (k))], when N _B When it is an odd number, g _B0 (n)=[q _B (k), fliplr(q _B (k-1)], fliplr is the sequence flip function in MatLab, G _B (e ^jω ) is the frequency response of g _B0 , α is the weight, the value is 0<α<1, m is the number of shifts, the value is any integer, δ(m) is the Dirac sequence, In the case of satisfying the constraint conditions, adjust the weight α, when the value of the objective function φ _B reaches the minimum, get half of the coefficient Q _B (k) of the optimized low-pass synthesis filter; 5) According to the half coefficient Q _B (k) of the band-pass synthesis filter after optimization, obtain the impulse response sequence g _B (n) of the band-pass synthesis filter g _B : When N _B is an odd number, g _B (n)=[Q _B (k), fliplr(Q _B (k-1))]; When N _B is an even number, g _B (n)=[Q _B (k), fliplr(Q _B (k))]; 6) According to the impulse response sequence g _B (n) of the band-pass synthesis filter and the impulse response sequence h _B (n) of the band-pass decomposition filter satisfy the time-domain inversion relation: h _B (n)=g _B (N _B -1-n), the impulse response sequence h _B (n) of the band-pass analysis filter can be obtained from the impulse response sequence g _B (n) of the band-pass synthesis filter. 3. A high-pass filter design method based on the variable bandwidth linear phase of Laplace structure, including the design of the high-pass synthesis filter _g and the design of the high-pass decomposition filter h _H : (A) Set the length N _H of the high-pass synthesis filter g _H , the passband cut-off frequency is The stopband start frequency is Bandwidth is Where M is the sampling factor, M and N _H are both integers, M≥2, N _H is an integer multiple of M, (B) According to the above parameters, call the firpm function in MatLab to generate the impulse response sequence p _H (n) of the high-pass prototype filter p _H of linear phase, n=0,1,2...N _H -1; (C) Determine half of the coefficients q _H (k) of the high-pass prototype filter p _H : When N _H is an odd number, When N _H is an even number, (D) Use q _H as the initial value of the optimization function fmincon, and optimize according to the following optimization formula: st<g _H0 (n), g _H0 (n-mM)>＝δ(m) Among them, g _H0 (n) is the high-pass filter impulse response free variable: when N _H is an even number, g _H0 (n)=[q _H (k),-fliplr(q _H (k))], when N _H When it is an odd number, g _H0 (n)=[q _H (k), fliplr(q _H (k-1)], fliplr is the sequence flip function in MatLab, G _H (e ^jω ) is the function of g _H0 (n) Frequency response, α is the weight, the value is 0<α<1, m is the number of shifts, the value is any integer, δ(m) is the Dirac sequence, In the case of satisfying the constraint conditions, adjust the weight α, when the value of the objective function φ _H reaches the minimum, get half of the coefficient Q _H (k) of the optimized high-pass synthesis filter; (E) According to half of the coefficients Q _H (k) of the optimized high-pass synthesis filter, the impulse response sequence g _H (n) of the high-pass synthesis filter g _H is obtained: When N _H is an odd number, g _H (n)=[Q _H (k), fliplr(Q _H (k-1))], When N _H is an even number, g _H (n)=[Q _H (k),-fliplr(Q _H (k))]; (F) According to the impulse response sequence g _H (n) of the high-pass synthesis filter and the impulse response sequence h _H (n) of the high-pass decomposition filter satisfy the time-domain inversion relation: h _H (n)=g _H (N _H - 1-n), the impulse response sequence h _H (n) of the high-pass decomposition filter can be obtained from the impulse response sequence g _H (n) of the high-pass synthesis filter.