CN105719255A - Variable bandwidth linear phase filter method based on Laplacian pyramid structure - Google Patents

Abstract

The invention discloses a design method applicable to a filter of a Laplacian pyramid structure and mainly solves the problems that existing filters of the Laplacian pyramid structure, except a Haar filter, cannot have a linear phase and orthogonal properties at the same time.According to the technical scheme, the method comprises the steps of firstly, setting the length of the filter to be N, setting the cut-off frequency of a passband to be omegap, setting the starting frequency of a stop band to be omegas, and setting a sampling factor to be M; then according to the parameters, calling a firpm function in MatLab to generate a prototype filter p; optimizing half of parameters of the prototype filter by calling a fmincon function, and obtaining a synthesized filter according to the optimization result; finally, obtaining a decomposition filter according to the time domain overturning relation.According to the method, a reasonable orthogonal constrained condition is given, accordingly, the filter can have orthogonality and linear phase characteristics, moreover, the bandwidth is variable, and wider application conditions are provided for the filter of the Laplacian pyramid structure.

Description

Bandwidth varying linear-phase filter method based on Laplace structure

Technical field

The invention belongs to signal processing technology field, particularly to the linear-phase filter method for designing of a kind of adaptive-bandwidth, can be used for image co-registration, increase, compression, denoising and rim detection.

Background technology

Laplacian pyramid structure LP is proposed at first by Burt et al., is used primarily for picture coding and compression.Afterwards, laplacian pyramid structure LP is expanded to the instrument of a kind of multiple dimensioned, multiresolution analysis by people, is applied to the various fields of image procossing, such as image co-registration, image denoising etc..

Fig. 1 gives laplacian pyramid structure chart.In FIG, for given input signal x, output c is the close approximation of x, and d is the residual error between original signal and its prediction signal p.In this structure, resolution filter h and composite filter g is usually a pair low pass filter, analysis matrix H and synthetical matrix G and then constitutes one group of transfer pair.

In traditional laplacian pyramid structure LP design, resolution filter h and composite filter g the two wave filter all take the low frequency channel of multichannel Perfect reconstruction filter banks, and filter bank structure is as shown in Figure 2.In each passage, wave filter and up and down decimation factor constitute positive and negative transfer pair.Therefore, in bank of filters, the wave filter of low frequency channel is usually applied directly in laplacian pyramid structure LP.Although, the structure that this way is laplacian pyramid structure LP provides help, but also limit the performance of laplacian pyramid structure LP median filter.

Laplacian pyramid structure LP plays an important role in the various fields of image procossing, but there is also a key issue at its design aspect, it may be assumed that the wave filter in existing laplacian pyramid structure LP is taken from the low frequency channel in multichannel Perfect reconstruction filter banks.This makes the wave filter in laplacian pyramid structure LP be subject to the constraint of Perfect reconstruction filter banks itself, and some characteristic is restricted, including linear phase characteristic very important in image procossing.In the design of Perfect reconstruction filter banks, although linear phase characteristic is easy in two-channel PR filter banks to realize, but is often difficulty with in multichannel Perfect reconstruction filter banks, causes except bandwidth isTwo-channel PR filter banks beyond, the problem that the Perfect reconstruction filter banks of other bandwidth is difficult to have linear phase characteristic.So, in laplacian pyramid structure LP, it is taken from this indirectly filter design method of Perfect reconstruction filter banks, it will cause that the wave filter in laplacian pyramid structure LP is subject to the restriction of Perfect reconstruction filter banks and is difficult to reach the deficiency of bandwidth varying and linear phase simultaneously.

Summary of the invention

Present invention aims to the deficiency of above-mentioned prior art, a kind of bandwidth varying linear-phase filter method based on Laplace structure is proposed, to break away from the design constraint of Perfect reconstruction filter banks, the frequency band making wave filter is variable, and has linear phase character while having orthogonal property.

The technical scheme is that and be achieved in that:

One. know-why

The decomposition part of laplacian pyramid structure is as it is shown in figure 1, when H and G meets H G=I relation, this passage in LP structure constitutes the subsystem of a Perfect Reconstruction, can realize the Perfect Reconstruction of antithetical phrase spacing wave.According to this relation of H G=I, resolution filter h and the composite filter g condition needing to meet can be derived, can be expressed as follows,

$\left\{\begin{array}{c}h\left(n\right)=g(N-1-n)\\ \⟨g\left(n\right),g(n-\mathrm{mM})\⟩=\mathrm{\δ}\left(m\right)\end{array}\right.$

Wherein, N is the length of wave filter, and M is decimation factor, and the bandwidth of wave filter isM is shift count, and value is arbitrary integer, and δ (m) is dirac sequence,

Above-mentioned condition enable to bandwidth varying wave filter there is orthogonal and linear phase characteristic simultaneously.

Two. technical scheme

According to above-mentioned principle, technical scheme is as follows:

Technical scheme 1: a kind of low pass filter design method of bandwidth varying linear phase based on Laplace structure, including to low pass composite filter g_LDesign and to low pass resolution filter h_LDesign:

(1) low pass composite filter g is set_LLength N_L, cut-off frequecy of passbandStopband initial frequencyBandwidth isWherein M is decimation factor, M and N_LIt is integer, M >=2, N_LFor the integral multiple of M, $0<{\mathrm{\&omega;}}_{p_L}<\frac{\mathrm{\&pi;}}{M}<{\mathrm{\&omega;}}_{s_L}<\mathrm{\&pi;};$

(2) according to above-mentioned parameter, call the firpm function in MatLab and produce the lowpass prototype filter p of linear phase_LLeast square estimation p_L(n), n=0,1,2 ... N_L-1；

(3) lowpass prototype filter p is determined_LHalf coefficient q_L(k):

Work as N_LDuring for odd number,

Work as N_LDuring for even number,

(4) by q_LK (), as the initial value of majorized function fmincon, is optimized according to the following formula that optimizes:

$\begin{array}{c}\underset{g_{L0}}{\min }{\mathrm{\&phi;}}_L=\mathrm{\&alpha;}{\&Integral;}_0^{{\mathrm{\&omega;}}_{\mathrm{PL}}}|1-G_L\left(e^{\mathrm{j\&omega;}}\right){|}^2\mathrm{d\&omega;}+(1-\mathrm{\&alpha;}){\&Integral;}_{{\mathrm{\&omega;}}_{s_L}}^{\mathrm{\&pi;}}|G_L\left(e^{\mathrm{j\&omega;}}\right){|}^2\mathrm{d\&omega;}\\ s.t.\&lang;g_{L0}\left(n\right),g_{L0}(n-\mathrm{mM})\&rang;=\mathrm{\&delta;}\left(m\right)\end{array}$

Wherein, g_L0N () is low pass filter impulse response free variable, work as N_LDuring for even number, g_L0(n)=[q_L(k), fliplr (q_L(k))], work as N_LDuring for odd number, g_L0(n)=[q_L(k), fliplr (q_L(k-1)], fliplr overturns function, G about sequence in MatLab_L(e^jw) it is g_L0The frequency response of (n), α is weight, and value is 0 < α < 1, m is shift count, and value is arbitrary integer, and δ (m) is dirac sequence,

When meeting constraints, adjust weight α, as object function φ_LLow-pass reconstruction filters g when value reaches minimum, after being optimized_LHalf coefficient Q_L(k)；

(5) the half coefficient Q according to the low pass composite filter after optimizing_LK (), draws low pass composite filter g_LLeast square estimation g_L(n):

Work as N_LDuring for odd number, g_L(n)=[Q_L(k), fliplr (Q_L(k-1))]；

Work as N_LDuring for even number, g_L(n)=[Q_L(k), fliplr (Q_L(k))]；

(6) the Least square estimation g according to low pass composite filter_LThe Least square estimation h of (n) and low pass resolution filter_LN () meets time domain inverse relation: h_L(n)=g_L(N_L-1-n), by the Least square estimation g of low pass composite filter_LN () can obtain the Least square estimation h of low pass resolution filter_L(n)。

Technical scheme 2: a kind of Design of Bandpass method of bandwidth varying linear phase based on Laplace structure, including to bandpass composite filter g_BDesign and to bandpass resolution filter h_BDesign:

1) bandpass composite filter g is set_BLength N_B, cut-off frequecy of passband is respectivelyStopband initial frequency is respectivelyDecimation factor is M, and bandwidth isWherein M and N_BIt is integer, M >=2, N_BIt is the integral multiple of M,

2) according to above-mentioned parameter, call the firpm function in MatLab and produce the bandpass ptototype filter p of linear phase_BLeast square estimation p_B(n), n=0,1,2 ... N_B-1；

3) bandpass ptototype filter p is determined_BHalf coefficient q_B(k):

Work as N_BDuring for odd number,

Work as N_BDuring for even number,

4) by q_BK (), as the initial value of majorized function fmincon, is optimized according to the following formula that optimizes:

$\begin{array}{c}\underset{g_{B0}}{\min }{\mathrm{\&phi;}}_B=(1-\mathrm{\&alpha;}){\&Integral;}_0^{{\mathrm{\&omega;}}_{p_{B1}}}|G_B\left(e^{\mathrm{j\&omega;}}\right){|}^2\mathrm{d\&omega;}+\mathrm{\&alpha;}{\&Integral;}_{{\mathrm{\&omega;}}_{s_{B1}}}^{{\mathrm{\&omega;}}_{p_{B2}}}|1-G_B\left(e^{\mathrm{j\&omega;}}\right){|}^2\mathrm{d\&omega;}+(1-\mathrm{\&alpha;}){\&Integral;}_{{\mathrm{\&omega;}}_{s_{B2}}}^{\mathrm{\&pi;}}|G_B\left(e^{\mathrm{j\&omega;}}\right){|}^2\mathrm{d\&omega;}\\ s.t.\&lang;g_{B0}\left(n\right),g_{B0}(n-\mathrm{mM})\&rang;=\mathrm{\&delta;}\left(m\right)\end{array}$

Wherein, g_B0N () is band filter impulse response free variable, work as N_BDuring for even number, g_B0(n)=[q_B(k), fliplr (q_B(k))], work as N_BDuring for odd number, g_B0(n)=[q_B(k), fliplr (q_B(k-1)], fliplr overturns function, G about sequence in MatLab_B(e^jw) it is g_B0Frequency response, α is weight, and value is 0 < α < 1, m is shift count, and value is arbitrary integer, and δ (m) is dirac sequence,

When meeting constraints, adjust weight α, as object function φ_BValue reaches minimum, the half coefficient Q of the low pass composite filter after being optimized_B(k)；

5) the half coefficient Q according to the bandpass composite filter after optimizing_BK (), draws bandpass composite filter g_BLeast square estimation g_B(n):

Work as N_BDuring for odd number, g_B(n)=[Q_B(k), fliplr (Q_B(k-1))]；

Work as N_BDuring for even number, g_B(n)=[Q_B(k), fliplr (Q_B(k))]；

6) the Least square estimation g according to bandpass composite filter_BThe Least square estimation h of (n) and bandpass resolution filter_BN () meets time domain inverse relation: h_B(n)=g_B(N_B-1-n), by the Least square estimation g of bandpass composite filter_BN () can obtain the Least square estimation h of bandpass resolution filter_B(n)。

Technical scheme 3: the high pass filter method for designing of a kind of bandwidth varying linear phase based on Laplace structure, including to high pass composite filter g_HDesign and to high pass resolution filter h_HDesign:

(A) high pass composite filter g is set_HLength N_H, cut-off frequecy of passband isStopband initial frequency isBandwidth isWherein M is decimation factor, M and N_HIt is integer, M >=2, N_HFor the integral multiple of M,

(B) according to above-mentioned parameter, call the firpm function in MatLab and produce the high pass ptototype filter p of linear phase_HLeast square estimation p_H(n), n=0,1,2 ... N_H-1；

(C) high pass ptototype filter p is determined_HHalf coefficient q_H(k):

Work as N_HDuring for odd number,

Work as N_HDuring for even number,

(D) by q_HAs the initial value of majorized function fmincon, it is optimized according to the following formula that optimizes:

$\begin{array}{c}\underset{g_{H0}}{\min }{\mathrm{\&phi;}}_H=(1-\mathrm{\&alpha;}){\&Integral;}_0^{{\mathrm{\&omega;}}_{s_H}}{\left|G_H\left(e^{\mathrm{j\&omega;}}\right)\right|}^2\mathrm{d\&omega;}+\mathrm{\&alpha;}{\&Integral;}_{{\mathrm{\&omega;}}_{p_H}}^{\mathrm{\&pi;}}{|1-G_H\left(e^{\mathrm{j\&omega;}}\right)|}^2\mathrm{d\&omega;}\\ s.t.<g_{H0}\left(n\right),g_{H0}(n-\mathrm{mM})>=\mathrm{\&delta;}\left(m\right)\end{array}$

Wherein, g_H0N () is high pass filter impulse response free variable: work as N_HDuring for even number, g_H0(n)=[q_H(k) ,-fliplr (q_H(k))], work as N_HFor time strange, g_H0(n)=[q_H(k), fliplr (q_H(k-1)], fliplr overturns function, G about sequence in MatLab_H(e^jw) it is g_H0The frequency response of (n), α is weight, and value is 0 < α < 1, m is shift count, and value is arbitrary integer, and δ (m) is dirac sequence,

When meeting constraints, adjust weight α, as object function φ_HWhen value reaches minimum, the half coefficient Q of the high pass composite filter after being optimized_H(k)；

(E) the half coefficient Q according to the high pass composite filter after optimizing_HK (), draws high pass composite filter g_HLeast square estimation g_H(n):

Work as N_HDuring for odd number, g_H(n)=[Q_H(k), fliplr (Q_H(k-1))],

Work as N_HDuring for even number, g_H(n)=[Q_H(k) ,-fliplr (Q_H(k))]；

(F) the Least square estimation g according to high pass composite filter_HThe Least square estimation h of (n) and high pass resolution filter_HN () meets time domain inverse relation: h_H(n)=g_H(N_H-1-n), by the Least square estimation g of high pass composite filter_HN () can obtain the Least square estimation h of high pass resolution filter_H(n)。

The present invention compared with prior art has the advantage that

First, the wave filter of present invention design can have linear phase and orthogonal property simultaneously.

The western Daubechies wave filter of existing many shellfishes is except Ha Er Haar wave filter, and portion only has orthogonal property, and does not have linear phase characteristic；Though biorthogonal Biorthogonal wave filter tool linear phase characteristic, but do not possess orthogonal property.Other all wave filter designed by Perfect reconstruction filter banks, except Haar wave filter, portion only has one of them of orthogonal property and linear phase characteristic, and can not have the two characteristic simultaneously.The present invention breaks away from the design constraint of Perfect reconstruction filter banks, directly the wave filter in design laplacian pyramid structure LP.When resolution filter h and reconfigurable filter g satisfies condition:Time, wave filter just can have linear phase and orthogonal property simultaneously；

Second, the adaptive-bandwidth of the present invention.

Existing Ha Er Haar filter bandwidht isThe filter bandwidht of the present invention can beNamely spectral limit is

3rd, design process of the present invention is simple.

Existing Ha Er Haar wave filter is that composite filter g and resolution filter h portion are designed, the wave filter of the present invention is only designed with to composite filter g, then according to Least square estimation h (n) of resolution filter h and composite filter g relation h (the n)=g (N-1-n) of Least square estimation g (n), resolution filter h can be tried to achieve, greatly reduce the complexity of design.

Accompanying drawing explanation

Fig. 1 is existing Laplce structure LP decomposition unit separation structure figure；

Fig. 2 is existing multichannel Perfect reconstruction filter banks structure chart；

Fig. 3 is the design flow diagram of the present invention；

Fig. 4 is the present invention based on the bandwidth of pulling structure isLength is the performance simulation figure of the low pass filter of 24；

Fig. 5 is the present invention based on the bandwidth of pulling structure isLength is the performance simulation figure of the band filter of 60；

Fig. 6 is the present invention based on the bandwidth of pulling structure isLength is the performance simulation figure of the high pass filter of 60.

Detailed description of the invention

Below in conjunction with accompanying drawing, technical scheme and effect are described in further detail.

The wave filter of present invention design is based on laplacian pyramid structure, as it is shown in figure 1, its wave filter includes resolution filter and composite filter.Input signal x sequentially passes through resolution filter h, after upper and lower decimation factor M and composite filter g, obtains the approximation signal p of input signal x, and the residual error between input signal and its approximation signal p is d.

Example 1: designing the bandwidth based on Laplace structure isLength is the low pass filter of 24.

This example includes low pass composite filter g_LDesign and to low pass resolution filter h_LDesign two parts.

One, design low-pass reconstruction filters g_L

Step 1, sets low pass composite filter g_LParameter.

If length N_L=24, decimation factor M=4, cut-off frequecy of passbandStopband initial frequencyWherein r_LIt is that low pass intermediate zone regulates parameter, by changing r_LSize, it is possible to the convenient width adjusting intermediate zone, this example takes r_L=0.4.

Step 2, it is determined that lowpass prototype filter p_L。

Lowpass prototype filter p_LA substantially low-pass impulse response sequence p_LN (), this example obtains linear phase low pass ptototype filter p by the firpm function called in MatLab_LLeast square estimation p_L(n), n=0,1 ... N_L-1。

Described firpm function is a function designing limited long impulse response sequence FIR filter, has introduction in " Digital Signal Processing " book, and the filter section that firpm function is designed has linear phase characteristic.

Step 3, it is determined that lowpass prototype filter p_LThe Least square estimation q of half_L(k)。

Due to lowpass prototype filter p_LThere is linear phase characteristic, i.e. its Least square estimation p_LN () is centrosymmetric, in order to reduce computation complexity, therefore only to p_LN the Least square estimation of the half of () is optimized and just can obtain low pass composite filter g_LLeast square estimation g_L(n)。

Take the Least square estimation p of lowpass prototype filter_LThe half Least square estimation q of (n)_LK (), wherein works as N_LWhen being even number,Work as N_LWhen being odd number,Such as, if N_L=4, p_L(n)=[p_L(0)p_L(1)p_L(2)p_L(3)], p is taken_LN the Least square estimation of the half of () is:This example takes N_L=24, the therefore Least square estimation p of lowpass prototype filter_LThe half Least square estimation of (n)

Step 4, sets the parameter of majorized function fmincon.

If the initial value of majorized function fmincon is q_LK (), constraints is < g_L0(n), g_L0(n-mM) >=δ (m), object function isWherein: g_L0N () is low pass filter impulse response free variable, m is shift count, and value is arbitrary integer, g_L0(n-mM) it is g_L0N m the shift sequence of (), α is weight, 0 < α < 1, G_L(e^jω) it is g_L0The frequency response of (n), δ (m) is dirac sequence, $\mathrm{\&delta;}\left(m\right)=\left\{\begin{array}{cc}1,& m=0\\ 0,& m\&NotEqual;0\end{array}\right.;$

Described majorized function fmincon is majorized function conventional in a kind of MatLab；

g_L0N the value of () is according to N_LParity and determine:

Work as N_LWhen being even number, g_L0(n)=[q_L(k), fliplr (q_L(k))],

Work as N_LWhen being odd number, g_L0(n)=[q_L(k), fliplr (q_L(k-1))],

Fliplr overturns function about sequence in MatLab, for instance set q_L(k)=[q_L(0)q_L(1)q_L(2)q_L(3)], N_L=4, by sequence q_LK the upset of () left and right, is represented by fliplr (q_L(k))=[q_L(3)q_L(2)q_L(1)q_L(0)]；

Constraints < g_L0(n), g_L0(n-mM) >=δ (m) shows g_L0(n) and g_L0(n-mM) sequence is orthogonal, namely $<g_{L0}\left(n\right),g_{L0}(n-\mathrm{mM})>=\left\{\begin{array}{c}1,m=0\\ 0,m\&NotEqual;0\end{array}\right.;$

Object function is the passband to wave filter and stopband optimizes simultaneously, and limited degree to the passband of wave filter and stopband when weight α determines to optimize, α shows more greatly passband restriction more big, and namely passband is more smooth, and stop band ripple is more big.

Step 5, the half impulse response q of the impulse response of lowpass prototype filter_LK () determines low pass composite filter g_LHalf coefficient Q_L(k)。

Impulse response p to lowpass prototype filter_LThe half impulse response q of (n)_LK (), is optimized by majorized function fmincon, adjust the size of weight α value, and this example takes α=0.8, when meeting constraints, as object function φ_LTime minimum, its function return value is low pass composite filter g_LHalf coefficient Q_L(k)。

Step 6, it is determined that low pass composite filter g_LLeast square estimation g_L(n)。

Due to low pass composite filter g_LIt is linear-phase filter, its Least square estimation g_LN () is centrosymmetric, therefore can by low pass composite filter g_LHalf coefficient Q_LK () determines low pass composite filter g_LLeast square estimation g_L(n):

Work as N_LWhen being even number, g_L(n)=[Q_L(k), fliplr (Q_L(k))],

Work as N_LWhen being odd number, g_L(n)=[Q_L(k), fliplr (Q_L(k-1))]。

Two, design low pass resolution filter h_L

Step 7, it is determined that the Least square estimation h of low pass resolution filter_L(n)。

Set the impulse response g of low pass composite filter_LThe Least square estimation h of (n) and low pass resolution filter_LN () meets time domain inverse relation: h_L(n)=g_L(N_L-1-n), the Least square estimation h of low pass resolution filter can be obtained according to this time domain inverse relation_L(n), i.e. low pass resolution filter h_L。

Such as set g_L(n)=[g_L(0)g_L(1)g_L(2)g_L(3)], N_L=4, g_L(n) and h_LN () meets time domain inverse relation, then h_L(n)=g_L(N_L-n-1)=[g_L(3)g_L(2)g_L(1)g_L(0)]。

By above step, it is possible to obtain based on the bandwidth of Laplace structure beLength is the low pass filter of 24, and frequency spectrum supports and ranges forThe performance of this low pass filter is as shown in Figure 4.Wherein:

Fig. 4 (a) is low pass composite filter g_LTime-domain pulse response,

Fig. 4 (b) is low pass composite filter g_LFrequency domain response,

Fig. 4 (c) is low pass resolution filter h_LTime-domain pulse response,

Fig. 4 (d) is low pass resolution filter h_LFrequency domain response.

Can be seen that from Fig. 4 (a) and 4 (c), this low pass filter belongs to linear-phase filter, demonstrates the Least square estimation g of the low pass composite filter of this example design_LN () meets < g_L(n),g_L(n-mM) >=δ (m), therefore this low pass composite filter has orthogonality, again due to the impulse response g of low pass composite filter_LThe Least square estimation h of (n) and low pass resolution filter_LN () meets time domain inverse relation: h_L(n)=g_L(N_L-1-n), the therefore impulse response h of low pass resolution filter_LN () also meets < h_L(n),h_L(n-mM) >=δ (m), namely this low pass resolution filter has orthogonality.Can be seen that from Fig. 4 (b) and 4 (d), this low pass filter frequency spectrum supports scope and isNamely bandwidth is

Example 2 designs the bandwidth based on Laplace structureLength is the band filter of 60.

This example includes bandpass composite filter g_BDesign and to bandpass resolution filter h_BDesign two parts.

1, design bandpass composite filter g_B

Step one, sets bandpass composite filter g_BParameter.

If length N_B=60, decimation factor M=3, cut-off frequecy of passband is respectively Stopband initial frequency is respectivelyWherein r_BIt is that bandpass intermediate zone regulates parameter, by changing r_BSize, it is possible to the convenient width adjusting intermediate zone, this reality r_B=0.3.

Step 2, it is determined that bandpass ptototype filter p_B。

Bandpass ptototype filter p_BA substantially bandpass Least square estimation p_BN (), this example obtains the logical ptototype filter p of linear phase, band by the firpm function called in MatLab_BLeast square estimation p_B(n), n=0,1 ... N_B-1。

Step 3, it is determined that bandpass ptototype filter p_BThe Least square estimation q of half_B(k)。

Due to bandpass ptototype filter p_BThere is linear phase characteristic, i.e. its Least square estimation p_BN () center is symmetrical, in order to reduce computation complexity, therefore only to p_BN the Least square estimation of () half is optimized and just can obtain bandpass composite filter g_BLeast square estimation g_B(n)。

Take the Least square estimation p of bandpass ptototype filter_BThe half Least square estimation q of (n)_BK (), wherein works as N_BWhen being even number,Work as N_BWhen being odd number,Such as, if N_B=4, p_B(n)=[p_B(0)p_B(1)p_B(2)p_B(3)], p is taken_BN the Least square estimation of the half of () is:This example takes N_B=60, the therefore Least square estimation p of bandpass ptototype filter_BThe half Least square estimation of (n)

Step 4, sets the parameter of majorized function fmincon.

If the initial value of majorized function fmincon is q_BK (), constraints is < g_B0(n), g_B0(n-mM) >=δ (m), object function isWherein: g_B0N () is band filter impulse response free variable, m is shift count, and value is arbitrary integer, g_B0(n-mM) it is g_B0N m the shift sequence of (), α is weight, 0 < α < 1, G_B(e^jω) it is g_B0The frequency response of (n), δ (m) is dirac sequence,

g_B0N the value of () is according to N_BParity and determine:

Work as N_BWhen being even number, g_B0(n)=[q_B(k), fliplr (q_B(k))],

Work as N_BWhen being odd number, g_B0(n)=[q_B(k), fliplr (q_B(k-1))],

Fliplr overturns function, constraints < g about sequence in MatLab_B0(n), g_B0(n-mM) >=δ (m) shows g_B0(n) and g_B0(n-mM) sequence is orthogonal, namely

Object function is the passband to wave filter and stopband optimizes simultaneously, and limited degree to the passband of wave filter and stopband when weight α determines to optimize, α shows more greatly passband restriction more big, and namely passband is more smooth, and stop band ripple is more big.

Step 5, the half impulse response q of the impulse response of bandpass ptototype filter_BK () determines low pass composite filter g_BHalf coefficient Q_L(k)。

Optimized by majorized function fmincon, the impulse response p to bandpass ptototype filter_BThe half impulse response q of (n)_BK (), adjusts the size of weight α value, this example takes α=0.7, when meeting constraints, as object function φ_BTime minimum, its function return value is bandpass composite filter g_BHalf coefficient Q_B(k)。

Step 6, it is determined that bandpass composite filter g_BLeast square estimation g_B(n)。

Due to bandpass composite filter g_BIt is linear-phase filter, its Least square estimation g_BN () is centrosymmetric, therefore can by bandpass composite filter g_BHalf coefficient Q_BK () determines bandpass composite filter g_BLeast square estimation g_B(n):

Work as N_BWhen being even number, g_B(n)=[Q_B(k), fliplr (Q_B(k))],

Work as N_BWhen being odd number, g_L(n)=[Q_B(k), fliplr (Q_B(k-1))]。

2, design bandpass resolution filter h_B

Step 7, it is determined that the Least square estimation h of bandpass resolution filter_B(n)。

According to the impulse response g setting bandpass composite filter_BThe Least square estimation h of (n) and bandpass resolution filter_BN () meets time domain inverse relation: h_B(n)=g_B(N_B-1-n), the Least square estimation h of resolution filter can be drawn_B(n), i.e. bandpass resolution filter h_B。

Such as set g_B(n)=[g_B(0)g_B(1)g_B(2)g_B(3)], N_B=4, g_B(n) and h_BN () meets time domain inverse relation, then h_B(n)=g_B(N_B-n-1)=[g_B(3)g_B(2)g_B(1)g_B(0)]。

By above step, it is possible to obtain based on the bandwidth of Laplace structure beLength is the band filter of 60, and frequency spectrum supports and ranges forThe performance of this band filter is as shown in Figure 4.Wherein:

Fig. 4 (a) is bandpass composite filter g_BTime-domain pulse response,

Fig. 4 (b) is bandpass composite filter g_BFrequency domain response,

Fig. 4 (c) is bandpass resolution filter h_BTime-domain pulse response,

Fig. 4 (d) is bandpass resolution filter h_BFrequency domain response.

Can be seen that from Fig. 4 (a) and 4 (c), this band filter belongs to linear-phase filter, demonstrates the Least square estimation g of the bandpass composite filter of this example design_BN () meets < g_B(n), g_B(n-mM) >=δ (m), therefore this bandpass composite filter has orthogonality, again because of the impulse response g of bandpass composite filter_BThe Least square estimation h of (n) and bandpass resolution filter_BN () meets time domain inverse relation: h_B(n)=g_B(N_B-1-n), the therefore impulse response h of bandpass resolution filter_BN () also meets < h_B(n), h_B(n-mM) >=δ (m), namely this bandpass resolution filter has orthogonality.Can be seen that from Fig. 4 (b) and 4 (d), this band filter frequency spectrum supports scope and isNamely bandwidth is

Example 3. designs the bandwidth based on Laplace structureLength is the high pass filter of 60.

This example includes high pass composite filter g_HDesign and to high pass resolution filter h_HDesign two parts.

1), design high pass composite filter g_H

Step A, sets high pass composite filter g_HParameter.

If length N_H=60, decimation factor M=3, cut-off frequecy of passbandStopband initial frequencyWherein r_HIt is that high pass intermediate zone regulates parameter, by changing r_HSize, it is possible to the convenient width adjusting intermediate zone, this example takes r_H=0.2.

Step B, it is determined that high pass ptototype filter p_H。

High pass ptototype filter p_HA substantially high pass Least square estimation p_HN (), this example obtains linear phase high pass ptototype filter p by the firpm function called in MatLab_HLeast square estimation p_H(n), n=0,1 ..., N_H-1。

Step C, it is determined that high pass ptototype filter p_HThe Least square estimation q of half_H(k)。

Due to high pass ptototype filter p_HThere is linear phase characteristic, i.e. its Least square estimation p_HN () is centrosymmetric, in order to reduce computation complexity, therefore only to p_HN the Least square estimation of the half of () is optimized and just can obtain high pass composite filter g_HLeast square estimation g_H(n)。

Take the Least square estimation p of high pass ptototype filter_HThe half Least square estimation q of (n)_HK (), wherein works as N_HWhen being even number,Work as N_HWhen being odd number,Such as, if N_H=4, p_H(n)=[p_H(0)p_H(1)p_H(2)p_H(3)], p is taken_HN the Least square estimation of the half of () is:This example takes N_H=24, the therefore Least square estimation p of high pass ptototype filter_HThe half Least square estimation of (n)

Step D, sets the parameter of majorized function fmincon.

If the initial value of majorized function fmincon is q_HK (), constraints is < g_H0(n), g_H0(n-mM) >=δ (m), object function isWherein: g_H0N () is high pass filter impulse response free variable, m is shift count, and value is arbitrary integer, g_H0(n-mM) it is g_H0N m the shift sequence of (), α is weight, 0 < α < 1, G_H(e^jω) it is g_H0The frequency response of (n), δ (m) is dirac sequence,

g_H0N the value of () is according to N_HParity and determine:

Work as N_HWhen being even number, g_H0(n)=[q_H(k) ,-fliplr (q_H(k))],

Work as N_HWhen being odd number, g_H0(n)=[q_H(k), fliplr (q_H(k-1))],

Fliplr overturns function about sequence in MatLab, for instance set q_H(k)=[q_H(0)q_H(1)q_H(2)q_H(3)], N_H=4, by sequence q_HK the upset of () left and right, is represented by fliplr (q_H(k))=[q_H(3)q_H(2)q_H(1)q_H(0)]；

Constraints < g_H0(n), g_H0(n-mM) >=δ (m) shows g_H0(n) and g_H0(n-mM) sequence is orthogonal, namely $<g_{H0}\left(n\right),g_{H0}(n-\mathrm{mM})>=\left\{\begin{array}{c}1,m=0\\ 0,m\&NotEqual;0\end{array}\right.;$

Object function is the passband to wave filter and stopband optimizes simultaneously, and limited degree to the passband of wave filter and stopband when weight α determines to optimize, α shows more greatly passband restriction more big, and namely passband is more smooth, and stop band ripple is more big.

Step E, the half impulse response q according to the impulse response of high pass ptototype filter_HK () determines high pass composite filter g_HHalf coefficient Q_H(k)。

Impulse response p to high pass ptototype filter_HThe half impulse response q of (n)_HK (), is optimized by majorized function fmincon, adjust the size of weight α value, and this example takes α=0.8, when meeting constraints, as object function φ_HTime minimum, its function return value is high pass composite filter g_HHalf coefficient Q_H(k)。

Step F, it is determined that high pass composite filter g_HLeast square estimation g_H(n)。

Due to high pass composite filter g_HIt is linear-phase filter, its Least square estimation g_HN () is centrosymmetric, therefore can by high pass composite filter g_HHalf coefficient Q_HK () can determine that high pass composite filter g_HLeast square estimation g_H(n):

Work as N_HWhen being even number, g_H(n)=[Q_H(k) ,-fliplr (Q_H(k))],

Work as N_HWhen being odd number, g_H(n)=[Q_H(k), fliplr (Q_H(k-1))]。

Step G, the impulse response g according to high pass composite filter_HN () draws high pass resolution filter h_H。

Owing to setting the impulse response g of high pass composite filter_HThe Least square estimation h of (n) and high pass resolution filter_HN () meets time domain inverse relation: h_H(n)=g_H(N_H-1-n), therefore the Least square estimation h of high pass resolution filter can be drawn by this time domain inverse relation_H(n), i.e. high pass resolution filter h_H。

Such as set g_H(n)=[g_H(0)g_H(1)g_H(2)g_H(3)], N_H=4, g_H(n) and h_HN () meets time domain inverse relation, then h_H(n)=g_H(N_H-n-1)=[g_H(3)g_H(2)g_H(1)g_H(0)]。

By above step, it is possible to obtain based on the bandwidth of Laplace structure beLength is the high pass filter of 60, and frequency spectrum supports and ranges forThe performance of this low pass filter is as shown in Figure 5.Wherein:

Fig. 5 (a) is high pass composite filter g_HTime-domain pulse response,

Fig. 5 (b) is high pass composite filter g_HFrequency domain response,

Fig. 5 (c) is high pass resolution filter h_HTime-domain pulse response,

Fig. 5 (d) is high pass resolution filter h_HFrequency domain response.

Can be seen that from Fig. 5 (a) and (c), this high pass filter belongs to linear-phase filter, demonstrates the Least square estimation g of the high pass composite filter of this example design_HN () meets < g_H(n), g_H(n-mM) >=δ (m), therefore this high pass composite filter has orthogonality, the impulse response g of high pass composite filter again_HThe Least square estimation h of (n) and high pass resolution filter_HN () meets time domain inverse relation: h_H(n)=g_H(N_H-1-n), the therefore impulse response h of high pass resolution filter_HN () also meets < h_H(n), h_H(n-mM) >=δ (m), namely this high pass resolution filter has orthogonality.Can be seen that from Fig. 5 (b) and 5 (d), this high pass filter frequency spectrum supports scope and isNamely bandwidth is

Claims (3)

1. based on a low pass filter design method for the bandwidth varying linear phase of Laplace structure, including to low pass composite filter g_LDesign and to low pass resolution filter h_LDesign:

(1) low pass composite filter g is set_LLength N_L, cut-off frequecy of passbandStopband initial frequencyBandwidth isWherein M is decimation factor, M and N_LIt is integer, M >=2, N_LFor the integral multiple of M, $0<{\mathrm{\&omega;}}_{p_L}<\frac{\mathrm{\&pi;}}{M}<{\mathrm{\&omega;}}_{s_L}<\mathrm{\&pi;};$

(2) according to above-mentioned parameter, call the firpm function in MatLab and produce the lowpass prototype filter p of linear phase_LLeast square estimation p_L(n), n=0,1,2 ... N_L-1；

(3) lowpass prototype filter p is determined_LHalf coefficient q_L(k):

Work as N_LDuring for odd number, $q_L\left(k\right)=p_L(0:\frac{N_L-1}{2}),k=0,1,2...\frac{N_L-1}{2},$

Work as N_LDuring for even number, $q_L\left(k\right)=p_L(0:\frac{N_L}{2}-1),k=0,1,2...\frac{N_L}{2}-1;$

(4) by q_LK (), as the initial value of majorized function fmincon, is optimized according to the following formula that optimizes:

$\underset{g_{L0}}{min}{\mathrm{\&phi;}}_L=\mathrm{\&alpha;}{\&Integral;}_0^{{\mathrm{\&omega;}}_{p_L}}|1-G_L\left(e^{j\mathrm{\&omega;}}\right){|}^{2}d\mathrm{\&omega;}+(1-\mathrm{\&alpha;}){{\&Integral;}_{{\mathrm{\&omega;}}_{s_L}}^{\mathrm{\&pi;}}\left|G_L\left(e^{j\mathrm{\&omega;}}\right)\right|}^{2}d\mathrm{\&omega;}$

s.t.〈g_L0(n),g_L0(n-mM) >=δ (m)

Wherein, g_L0N () is low pass filter impulse response free variable, work as N_LDuring for even number, g_L0(n)=[q_L(k),fliplr(q_L(k))], work as N_LDuring for odd number, g_L0(n)=[q_L(k),fliplr(q_L(k-1)], fliplr overturns function, G about sequence in MatLab_L(e^jw) it is g_L0The frequency response of (n), α is weight, and value is 0 < α < 1, m is shift count, and value is arbitrary integer, and δ (m) is dirac sequence, $\mathrm{\&delta;}\left(m\right)=\left\{\begin{array}{cc}1,& m=0\\ 0,& m\&NotEqual;0\end{array}\right.;$

When meeting constraints, adjust weight α, as object function φ_LLow pass composite filter g when value reaches minimum, after being optimized_LHalf coefficient Q_L(k)；

(5) the half coefficient Q according to the low pass composite filter after optimizing_LK (), draws low pass composite filter g_LLeast square estimation g_L(n):

Work as N_LDuring for odd number, g_L(n)=[Q_L(k),fliplr(Q_L(k-1))]；

Work as N_LDuring for even number, g_L(n)=[Q_L(k),fliplr(Q_L(k))]；

(6) the Least square estimation g according to low pass composite filter_LThe Least square estimation h of (n) and low pass resolution filter_LN () meets time domain inverse relation: h_L(n)=g_L(N_L-1-n), by the Least square estimation g of low pass composite filter_LN () can obtain the Least square estimation h of low pass resolution filter_L(n)。

2. based on a Design of Bandpass method for the bandwidth varying linear phase of Laplace structure, including to bandpass composite filter g_BDesign and to bandpass resolution filter h_BDesign:

1) bandpass composite filter g is set_BLength N_B, cut-off frequecy of passband is respectivelyStopband initial frequency is respectivelyDecimation factor is M, and bandwidth isWherein M and N_BIt is integer, M >=2, N_BIt is the integral multiple of M,

2) according to above-mentioned parameter, call the firpm function in MatLab and produce the bandpass ptototype filter p of linear phase_BLeast square estimation p_B(n), n=0,1,2 ... N_B-1；

3) bandpass ptototype filter p is determined_BHalf coefficient q_B(k):

Work as N_BDuring for odd number, $q_B\left(k\right)=p_B(0:\frac{N_B-1}{2}),k=0,1,2...\frac{N_B-1}{2},$

Work as N_BDuring for even number, $q_B\left(k\right)=p_B(0:\frac{N_B}{2}-1),k=0,1,2...\frac{N_B}{2}-1;$

4) by q_BK (), as the initial value of majorized function fmincon, is optimized according to the following formula that optimizes:

$\underset{g_{B0}}{min}{\mathrm{\&phi;}}_B=(1-\mathrm{\&alpha;}){\&Integral;}_0^{{\mathrm{\&omega;}}_{p_{B1}}}|G_B\left(e^{j\mathrm{\&omega;}}\right){|}^{2}d\mathrm{\&omega;}+\mathrm{\&alpha;}{\&Integral;}_{{\mathrm{\&omega;}}_{s_{B1}}}^{{\mathrm{\&omega;}}_{p_{B2}}}|1-G_B\left(e^{j\mathrm{\&omega;}}\right){|}^{2}d\mathrm{\&omega;}+(1-\mathrm{\&alpha;}){\&Integral;}_{{\mathrm{\&omega;}}_{s_{B2}}}^{\mathrm{\&pi;}}|G_B\left(e^{j\mathrm{\&omega;}}\right){|}^{2}d\mathrm{\&omega;}$

s.t.〈g_B0(n),g_B0(n-mM) >=δ (m)

Wherein, g_B0N () is band filter impulse response free variable, work as N_BDuring for even number, g_B0(n)=[q_B(k),fliplr(q_B(k))], work as N_BDuring for odd number, g_B0(n)=[q_B(k),fliplr(q_B(k-1)], fliplr overturns function, G about sequence in MatLab_B(e^jw) it is g_B0Frequency response, α is weight, and value is 0 < α < 1, m is shift count, and value is arbitrary integer, and δ (m) is dirac sequence, $\mathrm{\&delta;}\left(m\right)=\left\{\begin{array}{cc}1,& m=0\\ 0,& m\&NotEqual;0\end{array}\right.;$

When meeting constraints, adjust weight α, as object function φ_BValue reaches minimum, the half coefficient Q of the low pass composite filter after being optimized_B(k)；

5) the half coefficient Q according to the bandpass composite filter after optimizing_BK (), draws bandpass composite filter g_BLeast square estimation g_B(n):

Work as N_BDuring for odd number, g_B(n)=[Q_B(k),fliplr(Q_B(k-1))]；

Work as N_BDuring for even number, g_B(n)=[Q_B(k),fliplr(Q_B(k))]；

6) the Least square estimation g according to bandpass composite filter_BThe Least square estimation h of (n) and bandpass resolution filter_BN () meets time domain inverse relation: h_B(n)=g_B(N_B-1-n), by the Least square estimation g of bandpass composite filter_BN () can obtain the Least square estimation h of bandpass resolution filter_B(n)。

3. based on a high pass filter method for designing for the bandwidth varying linear phase of Laplace structure, including to high pass composite filter g_HDesign and to high pass resolution filter h_HDesign:

(A) high pass composite filter g is set_HLength N_H, cut-off frequecy of passband isStopband initial frequency isBandwidth isWherein M is decimation factor, M and N_HIt is integer, M >=2, N_HFor the integral multiple of M, $0<{\mathrm{\&omega;}}_{s_H}<\mathrm{\&pi;}-\frac{\mathrm{\&pi;}}{M}<{\mathrm{\&omega;}}_{p_H}<\mathrm{\&pi;};$

(B) according to above-mentioned parameter, call the firpm function in MatLab and produce the high pass ptototype filter p of linear phase_HLeast square estimation p_H(n), n=0,1,2 ... N_H-1；

(C) high pass ptototype filter p is determined_HHalf coefficient q_H(k):

Work as N_HDuring for odd number, $q_H\left(k\right)=p_H(0:\frac{N_H-1}{2}),k=0,1,2...\frac{N_H-1}{2},$

Work as N_HDuring for even number, $q_H\left(k\right)=p_H(0:\frac{N_H}{2}-1),k=0,1,2...\frac{N_H}{2}-1;$

(D) by q_HAs the initial value of majorized function fmincon, it is optimized according to the following formula that optimizes:

$\underset{g_{H0}}{min}{\mathrm{\&phi;}}_H=(1-\mathrm{\&alpha;}){\&Integral;}_0^{{\mathrm{\&omega;}}_{s_H}}|G_H\left(e^{j\mathrm{\&omega;}}\right){|}^{2}d\mathrm{\&omega;}+\mathrm{\&alpha;}{\&Integral;}_{{\mathrm{\&omega;}}_{p_H}}^{\mathrm{\&pi;}}|1-G_H\left(e^{j\mathrm{\&omega;}}\right){|}^{2}d\mathrm{\&omega;}$

s.t.〈g_H0(n),g_H0(n-mM) >=δ (m)

Wherein, g_H0N () is high pass filter impulse response free variable: work as N_HDuring for even number, g_H0(n)=[q_H(k),-fliplr(q_H(k))], work as N_HFor time strange, g_H0(n)=[q_H(k),fliplr(q_H(k-1)], fliplr overturns function, G about sequence in MatLab_H(e^jw) it is g_H0The frequency response of (n), α is weight, and value is 0 < α < 1, m is shift count, and value is arbitrary integer, and δ (m) is dirac sequence, $\mathrm{\&delta;}\left(m\right)=\left\{\begin{array}{cc}1,& m=0\\ 0,& m\&NotEqual;0\end{array}\right.;$

When meeting constraints, adjust weight α, as object function φ_HWhen value reaches minimum, the half coefficient Q of the high pass composite filter after being optimized_H(k)；

(E) the half coefficient Q according to the high pass composite filter after optimizing_HK (), draws high pass composite filter g_HLeast square estimation g_H(n):

Work as N_HDuring for odd number, $g_H\left(n\right)=\&lsqb;Q_H\left(k\right),fliplr\left(Q_H\left(k-1\right)\right)\&rsqb;,$

Work as N_HDuring for even number, $g_H\left(n\right)=\&lsqb;Q_H\left(k\right),-fliplr\left(Q_H\left(k\right)\right)\&rsqb;;$

(F) the Least square estimation g according to high pass composite filter_HThe Least square estimation h of (n) and high pass resolution filter_HN () meets time domain inverse relation: h_H(n)=g_H(N_H-1-n), by the Least square estimation g of high pass composite filter_HN () can obtain the Least square estimation h of high pass resolution filter_H(n)。