CN102882491B

CN102882491B - A kind of sparse method for designing without frequency deviation linear phase fir notch filter

Info

Publication number: CN102882491B
Application number: CN201210405631.2A
Authority: CN
Inventors: 赵加祥; 徐微; 王洪杰
Original assignee: Nankai University
Current assignee: Nankai University
Priority date: 2012-10-23
Filing date: 2012-10-23
Publication date: 2016-04-13
Anticipated expiration: 2032-10-23
Also published as: CN102882491A

Abstract

The invention provides a design method of a sparse frequency-offset linear phase FIR notch filter for the first time. The method uses LASSO algorithm to determine the number and position of non-zero tap coefficients of the required filter, and then applies the least square algorithm to calculate the values of these coefficients. The sparsity of this filter can reduce the number of adders and multipliers used in its realization, thereby improving its operation speed, reducing operation errors and reducing energy consumption. Simulation results show that under the requirements of the same design index, the linear phase FIR notch filter designed by the present invention is more than 15% less than the similar filter designed by the best design method at home and abroad, and its non-zero tap coefficient number is less than 15%. And there is no frequency offset problem.

Description

A Design Method of Sparse No Frequency Offset Linear Phase FIR Notch Filter

技术领域technical field

本发明属于数字信号处理技术领域，提供了一种稀疏、高效、无频偏的线性相位FIR(有限脉冲响应)陷波滤波器的设计方法。The invention belongs to the technical field of digital signal processing, and provides a design method of a sparse, high-efficiency, and frequency-offset-free linear phase FIR (finite impulse response) notch filter.

背景技术Background technique

随着信息技术的迅猛发展及集成电路制作的改进，数字信号处理技术因其灵活，快速，精确等优点，广泛应用于各个领域，它主要用于对语音信息处理，电话信道上的数据传输，图像处理，生物制造工程，地震学，核爆炸探测等等。在数字信号处理过程中，数字滤波技术是其重要的组成部分，而从脉冲响应角度数字信号滤波器一般可分为有限脉冲响应(FIR)滤波器和(无限脉冲响应)IIR滤波器。而从频率响应角度可分为以下几种：低通滤波器，高通滤波器，带通滤波器，带阻滤波器，陷波滤波器。陷波滤波器用于消除输入信号中某些特定频率并对剩余频率相对没有影响。在信号处理、通信和生物医学工程等领域有广泛应用。稀疏的线性相位FIR(有限脉冲响应)陷波滤波器是系数具有稀疏特性(非零抽头系数的数目小于滤波器阶数)的陷波滤波器。稀疏的滤波器其实现所用的加法和乘法器数目远少于与其滤波效果相当的同类滤波器，因此，稀疏的滤波器具有运算速度高、运算误差小和能耗低等优点。With the rapid development of information technology and the improvement of integrated circuit production, digital signal processing technology is widely used in various fields due to its advantages of flexibility, speed and accuracy. It is mainly used for voice information processing, data transmission on telephone channels, Image processing, biofabrication engineering, seismology, nuclear explosion detection, and more. In the process of digital signal processing, digital filtering technology is an important part, and from the perspective of impulse response, digital signal filters can generally be divided into finite impulse response (FIR) filters and (infinite impulse response) IIR filters. From the perspective of frequency response, it can be divided into the following types: low-pass filter, high-pass filter, band-pass filter, band-stop filter, and notch filter. Notch filters are used to remove certain frequencies from an input signal and have relatively no effect on the remaining frequencies. It is widely used in the fields of signal processing, communication and biomedical engineering. A sparse linear-phase FIR (Finite Impulse Response) notch filter is a notch filter whose coefficients have sparse properties (the number of non-zero tap coefficients is smaller than the filter order). The number of adders and multipliers used in the implementation of sparse filters is far less than that of similar filters with comparable filtering effects. Therefore, sparse filters have the advantages of high computing speed, small computing errors, and low energy consumption.

目前对于FIR数字滤波器的设计方法主要分为窗函数法，频率抽样设计法，最优化/计算机辅助设计法，而对于已经提出的陷波滤波器的设计方法中，P.Zahradnik和M.Vlcek的快速分析理想等纹波设计法可算是经典，通过分析方法，估计滤波器阶数，以及可设计非零系数较少，通带纹波较小的滤波器，之后的改进算法，精确等纹波设计法用来解决前者频率偏移的问题，但却以放宽带宽为代价。但P.Zahradnik和M.Vlcek的经典方法的优点不可忽视，也是国内外学者争相对比的模板。At present, the design methods for FIR digital filters are mainly divided into window function method, frequency sampling design method, optimization/computer-aided design method, and for the design method of notch filter that has been proposed, P.Zahradnik and M.Vlcek The ideal equal-ripple design method can be regarded as a classic. Through the analysis method, the order of the filter can be estimated, and a filter with less non-zero coefficients and smaller pass-band ripple can be designed. The wave design method is used to solve the former frequency offset problem, but at the cost of widening the bandwidth. However, the advantages of the classic methods of P.Zahradnik and M.Vlcek cannot be ignored, and it is also a template for domestic and foreign scholars to compare.

发明内容Contents of the invention

本发明目的是设计实现小纹波，低抽头数的陷波滤波器，同时克服上述快速分析理想等纹波设计法的频偏问题，并提供一种全新的设计方法——可设计稀疏，高效，无频偏的线性相位FIR陷波滤波器的方法。The purpose of the present invention is to design and realize a notch filter with small ripple and low tap number, overcome the frequency offset problem of the above-mentioned rapid analysis ideal equiripple design method at the same time, and provide a new design method-can design sparse, high-efficiency , a method for frequency-offset-free linear-phase FIR notch filters.

本发明提供的稀疏无频偏线性相位FIR陷波滤波器的设计方法具体步骤如下：The specific steps of the design method of the sparse frequency-offset linear phase FIR notch filter provided by the present invention are as follows:

(下面以I型线性相位FIR滤波器为例，本发明同样适用于其他II,III,IV型线性相位FIR滤波器)：(below taking I type linear phase FIR filter as example, the present invention is equally applicable to other II, III, IV type linear phase FIR filter):

(1)根据设计要求，包括通带纹波δ、阻带带宽陷波频率ω_s和衰减深度A_notch，确定线性相位FIR陷波滤波器的初始阶数N和相应的离散化理想频率响应；H_d(e^jω)表示理想线性相位FIR陷波滤波器频率响应，其中ω∈[0,π]，j表示虚数单位，离散化理想频率响应向量y代表H_d(e^jω)当ω依次取L+1个离散点{0,ω₀,2ω₀,3ω₀,L,lω₀L,Lω₀}时的取值，其中表示离散化后的频率间隔，0≤l≤L，L为正整数，y表示为：(1) According to design requirements, including passband ripple δ, stopband bandwidth The notch frequency ω _s and the attenuation depth A _notch determine the initial order N of the linear phase FIR notch filter and the corresponding discretized ideal frequency response; H _d (e ^jω ) represents the frequency response of the ideal linear phase FIR notch filter , where ω∈[0,π], j represents the imaginary unit, and the discretized ideal frequency response vector y represents H _d (e ^jω ) when ω takes L+1 discrete points {0,ω ₀ ,2ω ₀ ,3ω ₀ ,L,lω ₀ L,Lω ₀ }, where Indicates the discretized frequency interval, 0≤l≤L, L is a positive integer, and y is expressed as:

$y the y = = {[[{H h}_{d d} (({e e}^{j j 00})),, {H h}_{d d} (({e e}^{{jω jω}_{00}})),, ... ... {H h}_{d d} (({e e}^{{jlω jlω}_{00}})) ... ...,, {H h}_{d d} (({e e}^{{jLω jLω}_{00}}))]]}^{T T} - - - - - - ((11))$

线性相位FIR陷波滤波器抽头系数向量β表示为The linear phase FIR notch filter tap coefficient vector β is expressed as

β＝[β₀,β₁,...,β_m,...β_M]^T＝[h_M,2h_M-1,…,2h_m…,2h₁,2h₀]^T(2)β=[β ₀ ,β ₁ ,...,β _m ,...β _M ] ^T ＝[h _M ,2h _M-1 ,...,2h _m ...,2h ₁ ,2h ₀ ] ^T (2)

其中M＝N/2，h_m代表FIR陷波滤波器的第m个抽头系数,0≤m≤M。将稀疏的线性相位FIR陷波滤波器设计问题转化为如下的数学优化问题：Wherein M=N/2, h _m represents the mth tap coefficient of the FIR notch filter, 0≤m≤M. Transform the sparse linear phase FIR notch filter design problem into the following mathematical optimization problem:

$\hat{β} = \underset{β}{\arg \min} {| | β | |}_{0},$ 且满足 $\frac{{| | Xβ - y | |}_{2}}{{| | y | |}_{2}} < ϵ - - - (3)$ $\hat{β} = \underset{β}{\arg \min} {| | β | |}_{0},$ and satisfied $\frac{{| | Xβ - the y | |}_{2}}{{| | the y | |}_{2}} < ϵ - - - (3)$

其中||·||₀表示0-范数运算，||β||₀即表示抽头系数向量中非零抽头的个数，||·||₂表示2-范数运算；ε是给定误差；与滤波器设计相关的范德蒙矩阵X表示为Where ||·|| ₀ represents the 0-norm operation, ||β|| ₀ represents the number of non-zero taps in the tap coefficient vector, ||·|| ₂ represents the 2-norm operation; ε is the given error; the Vandermonde matrix X associated with the filter design is expressed as

其中x_m表示X中的列向量。where x _m represents the column vector in X.

(2)设定非零系数对应的列向量集合的初始值为空集，即X₀＝{}，FIR陷波滤波器抽头系数向量β的初始估计值为全0向量，对公式(4)中构造的范德蒙矩阵X进行压缩选择，选择与理想滤波器系数变量β相关性比较大的变量组成集合X_k(k＝0,1,...,M)。(2) Set the initial value of the column vector set corresponding to the non-zero coefficient to be an empty set, that is, X ₀ ={}, the initial estimated value of the FIR notch filter tap coefficient vector β is a vector of all 0s, compresses and selects the Vandermonde matrix X constructed in formula (4), and selects variables that have a relatively large correlation with the ideal filter coefficient variable β to form a set X _k (k=0,1,...,M ).

公式(4)中的矩阵X可以表示为X＝(x₀,x₁,...,x_m,...,x_M)，其中x_m(0≤m≤M)表示X的列向量。滤波器系数β的估计值可以表示为其中 ${\hat{β}}_{m} = {2 h}_{M - m}, (0 \leq m \leq M - 1) .$ 估计值参量为：The matrix X in formula (4) can be expressed as X=(x ₀ ,x ₁ ,...,x _m ,...,x _M ), where x _m (0≤m≤M) represents the column vector of X . The estimated value of the filter coefficient β can be expressed as in ${\hat{β}}_{m} = {2 h}_{m - m}, (0 \leq m \leq m - 1) .$ Estimated value parameter for:

$\overset{^^}{μ μ} = = X x \overset{^^}{β β} = = {Σ Σ}_{j j = = 11}^{m m} {x x}_{j j} {\overset{^^}{β β}}_{j j} - - - - - - ((88))$

当前变量与残差的相关系数记为 Current variables and residuals The correlation coefficient is recorded as

$c c ((\overset{^^}{μ μ})) = = {X x}^{T T} ((y the y - - \overset{^^}{μ μ})) - - - - - - ((99))$

先假设从开始，则此时残差就是y。相应此时相关系数为：找出与当前残差y相关系数绝对值最大的那个变量x_j1,设并将这个变量加入到集合X₁＝[x_j1]，集合设为X_k＝[...s_j1x_j1...,j1∈Ψ](s_j1是所选择的变量x_j1与当前残差的相关系数的符号)。Assume first from At the beginning, the residual at this time is y. Correspondingly, the correlation coefficient at this time is: Find the variable x _j1 with the largest absolute value of the correlation coefficient with the current residual y, set And add this variable to the set X ₁ =[x _j1 ], set the set as X _k =[...s _j1 x _j1 ...,j1∈Ψ] (s _j1 is the selected variable x _j1 and the current residual The sign of the poor correlation coefficient).

(3)沿着第一个变量x_j1方向前进，直到找到第二个向量x_j2，使它与当前残差的相关系数与原来的x_j1与当前残差的相关系数相等。下标加入Ψ，并将x_j2加入集合。此时估计值为其中保证与x_j,x_i有相同的相关性，即 (3) Advance along the direction of the first variable x _j1 until the second vector x _j2 is found, so that the correlation coefficient between it and the current residual is equal to the correlation coefficient between the original x _j1 and the current residual. The subscript joins Ψ, and x _j2 joins the set. The estimated value at this time is in ensure have the same correlation with x _j , x _i , namely

在第k(2≤k≤M)次迭代中，附录二中证明的由LASSO算法确定的当前最小角方向为：In the k-th (2≤k≤M) iteration, the current minimum angular direction determined by the LASSO algorithm proved in Appendix II is:

u_k＝u(k)＝X_kw_k(10)u _k ＝u(k)＝X _k w _k (10)

其中：in:

${w w}_{k k} = = {A A}_{k k} {G G}_{k k}^{- - 11} {I I}_{k k} = = (({w w}_{k k 11,,} {w w}_{k k 22},, ... ...,, {w w}_{k k j j},, .... ....)) - - - - - - ((1111))$

I_A＝(1,1,...,1)^T(12)I _A = (1,1,...,1) ^T (12)

${G G}_{k k} = = {X x}_{k k}^{T T} {X x}_{k k},, {A A}_{k k} = = {(({I I}_{k k}^{T T} {G G}_{k k}^{- - 11} {I I}_{k k}))}^{- - \frac{11}{22}} - - - - - - ((1313))$

沿着最小角方向在X中继续寻找新的列向量，直到找到与当前残差相关系数绝对值与原来集合中的变量与当前残差相关系数绝对值相等的变量出现，加入集合X_k，确定当前每次前进的步长，Continue to search for new column vectors in X along the direction of the minimum angle until a variable with the absolute value of the correlation coefficient of the current residual error and the variable in the original set is equal to the absolute value of the current residual correlation coefficient appears, and join the set X _k to determine The current step size of each advance,

$\overset{^^}{γ γ} = = {min min}_{j j &Element; &Element; Ψ Ψ}^{+ +} {{\frac{c c - - {\overset{^^}{c c}}_{j j}}{{A A}_{k k} - - {a a}_{j j}},, \frac{c c + + {\overset{^^}{c c}}_{j j}}{{A A}_{k k} + + {a a}_{j j}}}} - - - - - - ((1414))$

而and

$\overset{~ ~}{γ γ} = = {min min}_{{γ γ}_{j j} > > 00} {{{γ γ}_{j j}}},, {γ γ}_{j j} = = - - {\overset{^^}{β β}}_{j j} / / {\overset{^^}{d d}}_{j j} - - - - - - ((1515))$

其中A_k＝X^Tμ＝(A_k1,A_k2,...,A_km)， where A _k =X ^T μ = (A _k1 ,A _k2 ,...,A _km ),

进行判断： Make a judgment:

回到(3)继续寻找下一个变量，当集合X_k包含了X中的所有列向量(即(x₀,x₁,...,x_m,...,x_M)这些列向量)时，集合X_k的构造完成。Go back to (3) and continue to look for the next variable, when the set X _k contains all column vectors in X (ie (x ₀ , x ₁ ,...,x _m ,...,x _M ) these column vectors) When , the construction of the set X _k is completed.

(4)利用当前的集合X_k，以及利用最小二乘法估计FIR陷波滤波器抽头系数向量β。根据预设稀疏度门限值，在所有集合X_k(k＝0,1,...,M)中确定适当的集合X_i，集合X_i＝[x_j1...x_ji]0≤i≤M中的各个向量的数目和位置即对应着滤波器非零抽头系数的数目和位置，X_i代表对应的陷波滤波器非零抽头系数有i个。则FIR陷波滤波器抽头系数向量：(4) Estimate the FIR notch filter tap coefficient vector β by using the current set X _k and the least square method. According to the preset sparsity threshold, determine the appropriate set Xi in all sets X _k ( _k =0,1,...,M ₎ , set Xi =[x _j1 ... x _ji ]0≤ The number and position of each vector in i≤M corresponds to the number and position of the non-zero tap coefficients of the filter, and X _i represents that there are i non-zero tap coefficients of the corresponding notch filter. Then the FIR notch filter tap coefficient vector:

$\overset{^^}{β β} = = {(({X x}_{i i} {X x}_{i i}^{T T}))}^{- - 11} {X x}_{i i}^{T T} y the y - - - - - - ((55))$

(5)对估计值进行转化，从而得到滤波器系数：(5) For the estimated value Transform to get the filter coefficients:

$\overset{^^}{β β} = = {[[{\overset{^^}{β β}}_{00},, {\overset{^^}{β β}}_{11},, ... ...,, {\overset{^^}{β β}}_{m m},, ... ... {\overset{^^}{β β}}_{M m}]]}^{T T} = = {(({h h}_{M m},, 22 {h h}_{M m - - 11},, ... ...,, 22 {h h}_{00}))}^{T T} - - - - - - ((66))$

(6)对于给定误差ε，判断如下不等式是否成立：(6) For a given error ε, judge whether the following inequality is true:

$\frac{| | | | X x \overset{^^}{β β} - - y the y | | {| |}_{22}}{| | | | y the y | | {| |}_{22}} < < ϵ ϵ - - - - - - ((77))$

若上式成立，则即为本发明所设计滤波器的抽头系数的具体数值，那么得到的FIR陷波滤波器的频率响应记为G(e^jω)：If the above formula holds, then Be the specific numerical value of the tap coefficient of the designed filter of the present invention, so the frequency response of the FIR notch filter that obtains is denoted as G(e ^jω ):

$G G (({e e}^{j j w w})) = = {h h}_{00} + + 22 {Σ Σ}_{m m = = 11}^{M m} {h h}_{m m} {e e}^{- - j j m m w w} - - - - - - ((1616))$

否则回到步骤(3)，继续增加所需集合，进而确定稀疏的线性相位FIR陷波滤波器所需的滤波器非零抽头系数的数目、位置以及具体的数值，迭代过程在计算得到的FIR滤波器满足设计要求时停止。Otherwise, go back to step (3), continue to increase the required set, and then determine the number, position and specific value of the non-zero tap coefficients of the filter required by the sparse linear phase FIR notch filter. The iterative process is obtained in the calculated FIR The filter stops when it meets the design requirements.

本发明的优点和积极效果：Advantage and positive effect of the present invention:

1、本发明首次提供了一种稀疏、无偏的线性相位FIR陷波滤波器设计方法；2、本发明可设计低非零抽头数的陷波滤波器，滤波器的稀疏性可使其实现所用的加法器乘法器数目减少，从而能提高其运算速度、减小运算误差和降低能耗，进而降低生产成本；3、仿真结果表明，在相同设计指标的要求下，本发明的非零抽头系数的数目比国内外最佳的同类滤波器的数目少15％以上，并且没有出现频偏问题。1, the present invention provides a kind of sparse, unbiased linear phase FIR notch filter design method for the first time; 2, the present invention can design the notch filter of low non-zero tap number, the sparsity of filter can make it realize The used adder multiplier number reduces, thereby can improve its operation speed, reduce operation error and reduce energy consumption, and then reduce production cost; 3, simulation result shows, under the requirement of same design index, non-zero tap of the present invention The number of coefficients is more than 15% less than that of the best filters of the same kind at home and abroad, and there is no frequency deviation problem.

附图说明Description of drawings

图1是：实现本发明的基于LASSO和最小二乘算法的线性相位FIR陷波滤波器设计方法流程图。Fig. 1 is: realize the flow chart of the linear phase FIR notch filter design method based on LASSO and least square algorithm of the present invention.

图2是：线性相位FIR陷波滤波器设计参数解释图。Figure 2 is an explanatory diagram of design parameters of a linear phase FIR notch filter.

图3是：利用本发明进行陷波滤波器设计的实施例一，可见通过本发明的设计，解决了频偏问题。同时抽头系数更少。FIG. 3 is the first embodiment of notch filter design using the present invention. It can be seen that the frequency offset problem is solved through the design of the present invention. At the same time, there are fewer tap coefficients.

具体实施方式detailed description

实施例1：Example 1:

第1、稀疏的线性相位FIR陷波滤波器阶数的估计；1. Estimation of sparse linear phase FIR notch filter order;

第2、稀疏的线性相位FIR陷波滤波器设计方法，包括单位脉冲响应的非零抽头系数数目、位置以及具系数数值的确定。2. A sparse linear phase FIR notch filter design method, including the determination of the number, position and value of non-zero tap coefficients of the unit impulse response.

为了验证该滤波器设计方法的有效性，对该方法进行了计算机模拟仿真。In order to verify the validity of the filter design method, a computer simulation is carried out on the method.

设计要求：利用文献[12]中所给出的设计指标，陷波频率为0.84π，阻带带宽0.061π，通带纹波-0.95dB.我们用LASSO和最小二乘结合算法来设计滤波器。Design requirements: Using the design indicators given in literature [12], the notch frequency is 0.84π, the stopband bandwidth is 0.061π, and the passband ripple is -0.95dB. We use the combination of LASSO and least squares algorithm to design the filter .

步骤一：根据FIR陷波滤波器的设计参数，估算陷波滤波器H_d(e^jω)的初始阶数N，通过计算可以得初始阶数应为：N＝76，即对于此种设计要求，本发明能够设计的滤波器阶数为76阶。Step 1: According to the design parameters of the FIR notch filter, estimate the initial order N of the notch filter H _d (e ^jω ), and the initial order N can be obtained by calculation: N=76, that is, for this design requirement , the filter order that the present invention can design is 76 orders.

步骤二：根据FIR陷波滤波器的设计参数要求，理想陷波滤波器的频率响应H_d(e^jω)如下：Step 2: According to the design parameter requirements of the FIR notch filter, the frequency response H _d (e ^jω ) of the ideal notch filter is as follows:

${H h}_{d d} (({e e}^{j j ω ω})) = = \{\begin{matrix} 11;; & 00 \leq \leq ω ω < < {ω ω}_{s the s} - - Δ Δ F f \\ 00;; & {ω ω}_{s the s} - - Δ Δ F f \leq \leq ω ω \leq \leq {ω ω}_{s the s} + + Δ Δ F f \\ 11;; & {ω ω}_{s the s} + + Δ Δ F f < < ω ω \leq \leq π π \end{matrix}$

其中陷波频率ω_s＝0.84π，阻带带宽通带纹波δ＝-0.95dB。Where the notch frequency ω _s =0.84π, the stopband bandwidth Passband ripple δ=-0.95dB.

步骤三：使用LASSO算法和最小二乘算法求解步骤一中构造的L-2范数误差最小的问题，进而得到稀疏的FIR陷波滤波器的系数。Step 3: Use the LASSO algorithm and the least squares algorithm to solve the problem of the smallest error of the L-2 norm constructed in step 1, and then obtain the coefficients of the sparse FIR notch filter.

利用本发明算法对变量集X进行变量压缩选择。对于76阶滤波器，本算法进行选择后可得到所有抽头系数的所有情况，即变量集X_k(k＝0,1,...,76)。The algorithm of the present invention is used to perform variable compression selection on the variable set X. For the 76th-order filter, after selection by this algorithm, all cases of all tap coefficients can be obtained, that is, the variable set X _k (k=0, 1, . . . , 76).

根据设定的稀疏门限值，在得到的变量集中，选择适当的变量集X_j，利用最小二乘算法进行估计逼近，得到的滤波器系数为：According to the set sparse threshold value, in the obtained variable set, select the appropriate variable set X _j , and use the least square algorithm to estimate and approximate, and the obtained filter coefficient is:

$\overset{^^}{β β} = = {[[{\overset{^^}{β β}}_{00},, {\overset{^^}{β β}}_{11},, ... ...,, {\overset{^^}{β β}}_{M m}]]}^{T T} = = {(({h h}_{M m},, 22 {h h}_{M m - - 11},, ... ...,, 22 {h h}_{00}))}^{T T} = = {(({X x}_{j j} {X x}_{j j}^{T T}))}^{- - 11} {X x}_{j j}^{T T} y the y$

再计算误差，若误差小于给定误差：Then calculate the error, if the error is less than the given error:

$\frac{| | | | X x \overset{^^}{β β} - - y the y | | {| |}_{22}}{| | | | y the y | | {| |}_{22}} < < ϵ ϵ$

若上式成立，那么相应的FIR滤波器的频率响应G(e^jω)：If the above formula is established, then the frequency response G(e ^jω ) of the corresponding FIR filter:

$G G (({e e}^{j j w w})) = = {h h}_{00} + + 22 {Σ Σ}_{n no = = 11}^{M m} {h h}_{n no} {e e}^{- - j j n no w w}$

否则若误差大于给定误差，则继续从所选的所有变量集X_k(k＝0,1,...,76)中选择适当的，进行最小二乘估计，当得到满足设计要求的FIR滤波器时停止迭代。最终得到的线性相位FIR陷波滤波器是非零系数为63的滤波器。与快速分析理想等纹波设计法相比节省了18.2％，同时解决了频偏问题。Otherwise, if the error is greater than the given error, continue to select the appropriate one from all the selected variable sets X _k (k=0,1,...,76), and perform the least squares estimation, when the FIR that meets the design requirements is obtained Stop iterations when filtering. The resulting linear-phase FIR notch filter is a filter with nonzero coefficients of 63. Compared with the fast analysis ideal equiripple design method, it saves 18.2%, and solves the problem of frequency deviation at the same time.

在表-1中，我们分别比较了两种算法得到的FIR陷波滤波器的阶数、非零抽头权重的数量、阻带带宽(在δ＝-0.95dB处)、通带纹波和陷波频率ω_s处的衰减等几项关键指标。In Table-1, we compared the order of the FIR notch filter obtained by the two algorithms, the number of non-zero tap weights, the stopband bandwidth (at δ = -0.95dB), the passband ripple and the notch Several key indicators such as the attenuation at the wave frequency ω _s .

表-1Table 1

在表-2中，我们列出了利用本发明中所涉及的算法所设计了线性相位FIR陷波滤波器的时域脉冲响应系数。In Table-2, we list the time-domain impulse response coefficients of the linear-phase FIR notch filter designed using the algorithm involved in the present invention.

表-2Table 2

利用本发明算法设计例子(1)得到的陷波滤波器的脉冲响应The impulse response of the notch filter that utilizes the algorithm design example (1) of the present invention to obtain

陷波滤波器的理想频响H_d(e^jω)可以表示为：The ideal frequency response H _d (e ^jω ) of the notch filter can be expressed as:

附录一Appendix I

式(1)理想陷波滤波器频响函数离散化的证明陷波滤波器的理想频响H_d(e^jω)可以表示为：Equation (1) Proof of Discretization of Frequency Response Function of Ideal Notch Filter The ideal frequency response H _d (e ^jω ) of notch filter can be expressed as:

其中j表示虚数单位，ω表示频率， where j represents the imaginary unit, ω represents the frequency,

I型N阶线性相位FIR陷波滤波器的频响函数H(e^jω)表示为：The frequency response function H(e ^jω ) of type I N-order linear phase FIR notch filter is expressed as:

$H h (({e e}^{j j ω ω})) = = {Σ Σ}_{n no = = 00}^{N N} {h h}_{n no} {e e}^{- - j j ω ω n no} = = {e e}^{- - j j M m ω ω} {H h}_{z z e e r r o o} (({e e}^{j j ω ω}))$

其中h_n(0≤n≤N)是滤波器H(e^jω)的时域冲激响应，且满足h_n＝h_N-n(0≤n≤N，N是正偶数)。H_zero(e^jω)是滤波器H(e^jω)的零相位频响，其定义如下：Where h _n (0≤n≤N) is the time-domain impulse response of the filter H(e ^jω ), and h _n =h _Nn (0≤n≤N, N is a positive even number). H _zero (e ^jω ) is the zero-phase frequency response of filter H(e ^jω ), which is defined as follows:

$\begin{matrix} {H h}_{z z e e r r o o} (({e e}^{j j ω ω})) = = {h h}_{M m} + + 22 {Σ Σ}_{m m = = 11}^{M m} {h h}_{M m - - m m} c c o o s the s ((m m ω ω)) \\ = = ((11,, cos cos ((ω ω)),, ... ...,, cos cos ((M m ω ω)))) {(({h h}_{M m},, 22 {h h}_{M m - - 11},, ... ...,, 22 {h h}_{00}))}^{T T} \end{matrix}$

其中M＝N/2。那么I型线性相位FIR陷波滤波器的设计问题可以转化为如下的L-2范数误差最小问题：where M=N/2. Then the design problem of type I linear phase FIR notch filter can be transformed into the following problem of minimum L-2 norm error:

$\frac{{| | | | {H h}_{zero zero} (({e e}^{jω jω})) - - {H h}_{f f} (({e e}^{jω jω})) | | | |}_{22}}{{| | | | {H h}_{d d} (({e e}^{jω jω})) | | | |}_{22}} < < ϵ ϵ,, &ForAll; &ForAll; ω ω &Element; &Element; [[00,, π π]]$

对频率ω进行均匀采样，则上述问题的离散表达形式为：The frequency ω is uniformly sampled, then the discrete expression of the above problem is:

$\frac{{| | | | Xβ Xβ - - y the y | | | |}_{22}}{{| | | | y the y | | | |}_{22}} < < ϵ ϵ$

其中：in:

β＝[h_M,2h_M-1,…,2h_m…,2h₁,2h₀]^T β=[h _M ,2h _M-1 ,…,2h _m …,2h ₁ ,2h ₀ ] ^T

$y the y = = {[[{H h}_{d d} (({e e}^{j j 00})),, ... ...,, {H h}_{d d} (({e e}^{{jω jω}_{00}})),, ... ...,, {H h}_{d d} (({e e}^{{jLω jLω}_{00}}))]]}^{T T}$

其中L与总的采样点数有关(总的采样点数为L+1)，0≤l≤L，0≤m≤M。Where L is related to the total number of sampling points (the total number of sampling points is L+1), 0≤l≤L, 0≤m≤M.

附录二Appendix II

式(9)LASSO算法确定最小角方向的证明Equation (9) LASSO algorithm to determine the proof of the minimum angle direction

利用LASSO算法确定最小角方向的具体实现步骤如下：The specific implementation steps of using the LASSO algorithm to determine the direction of the minimum angle are as follows:

矩阵X可以表示为X＝(x₀,x₁,...,x_m,...,x_M)，其中x_m(0≤m≤M)表示X的列向量。滤波器系数β的估计值可以表示为其中 ${\hat{β}}_{0} = h_{M}, {\hat{β}}_{m} = {2 h}_{M - m},$ $(0 < m \leq M - 1) .$ 定义参量为The matrix X can be expressed as X=(x ₀ , x ₁ , . . . , x _m , . . . , x _M ), where x _m (0≤m≤M) represents a column vector of X. The estimated value of the filter coefficient β can be expressed as in ${\hat{β}}_{0} = h_{m}, {\hat{β}}_{m} = {2 h}_{m - m},$ $(0 < m \leq m - 1) .$ define parameters for

$\overset{^^}{μ μ} = = X x \overset{^^}{β β} = = {Σ Σ}_{j j = = 11}^{m m} {x x}_{j j} {\overset{^^}{β β}}_{j j}$

当前变量与残差之间的相关系数记为 The correlation coefficient between the current variable and the residual is recorded as

$c c ((\overset{^^}{μ μ})) = = {X x}^{T T} ((y the y - - \overset{^^}{μ μ}))$

先假设从开始，则此时残差就是y。相应此时相关系数为：找出与当前残差y相关系数绝对值最大的那个变量x_j1,设并将这个变量加入到活动变量集中，活动变量集设为X_k＝[...s_j1x_j1...,j1∈Ψ](s_j1是所选择的变量x_j1与当前残差的相关系数的符号)，继续寻找相关性相同的变量加入到活动变量集X_k＝[...s_j1x_j1...,j1∈Ψ]中。Assume first from At the beginning, the residual at this time is y. Correspondingly, the correlation coefficient at this time is: Find the variable x _j1 with the largest absolute value of the correlation coefficient with the current residual y, set And add this variable to the active variable set, the active variable set is set to X _k =[...s _j1 x _j1 ...,j1∈Ψ] (s _j1 is the correlation between the selected variable x _j1 and the current residual The sign of the coefficient), and continue to look for variables with the same correlation to add to the active variable set X _k =[...s _j1 x _j1 ...,j1∈Ψ].

设Assume

I_A＝(1,1,...,1)^T I _A ＝(1,1,...,1) ^T

${G G}_{k k} = = {X x}_{k k}^{T T} {X x}_{k k},, {A A}_{k k} = = {(({I I}_{k k}^{T T} {G G}_{k k}^{- - 11} {I I}_{k k}))}^{- - \frac{11}{22}}$

${W W}_{k k} = = {A A}_{k k} {G G}_{k k}^{- - 11} {I I}_{k k} = = (({w w}_{k k 11} {w w}_{k k 22},, ... ...,, {w w}_{k k j j},, .... ....))$

那么当前的最小角方向为：Then the current minimum angular direction is:

u_k＝u(k)＝X_kw_k u _k ＝u(k)＝X _k w _k

因而A_kI_k中的元素都相同，所以这样的选择可以保证当前已选入集合的变量与最小角方向之间的相关系数均相等，即相当于已选入变量集的变量的“角平分线”，保证了LARS的要求。because The elements in A _k I _k are all the same, so this selection can ensure that the correlation coefficients between the variables currently selected into the set and the direction of the minimum angle are equal, which is equivalent to the "angle bisection" of the variables that have been selected into the variable set line", which guarantees the requirements of LARS.

参考文献：references:

[1].刘树棠.黄建国.离散时间信号处理.西安交通大学出版社.2001.9.[1]. Liu Shutang. Huang Jianguo. Discrete Time Signal Processing. Xi'an Jiaotong University Press. 2001.9.

[2].王兵锋.FIR数字滤波器设计与仿真：[硕士论文].哈尔滨理工大学.2009.3.[2]. Wang Bingfeng. FIR Digital Filter Design and Simulation: [Master Thesis]. Harbin University of Science and Technology. 2009.3.

[3].SuhashChandraDuttaRoyandBalbirKumaretal,‘notchfilterdesignreview’,ElectronicsandEnergeticsVOL.14,No.3,December2001,295-327.[3]. Suhash Chandra Dutta Roy and Balbir Kumaretal, 'notch filter design review', Electronics and Energetics VOL.14, No.3, December 2001, 295-327.

[4].Wei,D.；,"Non-convexoptimizationforthedesignofsparsefirfilters,"StatisticalSignalProcessing,2009.SSP'09.IEEE/SP15thWorkshopon,vol.,no.,pp.117-120,Aug.312009-Sept.32009[4].Wei,D.;,"Non-convex optimization for the design of sparse filters,"Statistical Signal Processing,2009.SSP'09.IEEE/SP15thWorkshopon,vol.,no.,pp.117-120,Aug.312009-Sept.32009

[5].Deshpande,R.；Kumar,B.；Jain,S.B.；,"Highlynarrowrejectionbandwidthfiniteimpulseresponsenotchfiltersforcommunication,"Communications,IET,vol.4,no.18,pp.2208-2216,December172010[5]. Deshpande, R.; Kumar, B.; Jain, S.B.;, "Highlynarrow rejectionbandwidthfiniteimpulseresponsenotchfiltersforcommunication," Communications, IET, vol.4, no.18, pp.2208-2216, December 172010

[6].Zahradnik,P.；Vlcek,M.；,"FastanalyticaldesignalgorithmsforFIRnotchfilters,"CircuitsandSystemsI:RegularPapers,IEEETransactionson,vol.51,no.3,pp.608-623,March2004[6]. Zahradnik, P.; Vlcek, M.;, "Fast analytical design algorithms for FIR notch filters," Circuits and Systems I: Regular Papers, IEEE Transactionson, vol.51, no.3, pp.608-623, March 2004

[7].S.C.DuttaRoy,S.B.JainandB.Kumar,'DesignofdigitalFIRnotchfilters'.IEEProc.-Vis.ImageSignalProcess,Vot.141,NO.5,October1994.[7]. S.C.DuttaRoy, S.B.JainandB.Kumar, 'DesignofdigitalFIRnotchfilters'. IEEProc.-Vis.ImageSignalProcess, Vot.141, NO.5, October1994.

[8].Chien-ChengTsengandSoo-ChangPei.'SparseFIRnotchfilterdesignanditsapplication'.IEE1997,ElectronicsLettersOnlineNO:19970766.[8]. Chien-Cheng Tseng and Soo-Chang Pei. 'Sparse FIR notch filter design and its application'. IEE1997, Electronics Letters Online NO: 19970766.

[9].Zahradnik,P.；Vlcek,M.；,"AnanalyticalprocedureforcriticalfrequencytuningofFIRfilters,"CircuitsandSystemsII:ExpressBriefs,IEEETransactionson,vol.53,no.1,pp.72-76,Jan.2006[9]. Zahradnik, P.; Vlcek, M.;, "Ananalytical procedure for critical frequency tuning of FIR filters," Circuits and Systems II: Express Briefs, IEEE Transactionson, vol.53, no.1, pp.72-76, Jan.2006

[10].KIchige,MIwaki,andRIshii.AccurateestimationofminimumfilterlengthforoptimumFIRdigitalfilters[J],CircuitsandSystemsII:AnalogandDigitalSignalProcessing,2000,47(10):1008-1016.[10]. KIchige, MIwaki, and RIshii. Accurate estimation of minimum filter length for optimal FIR digital filters [J], Circuits and Systems II: Analog and Digital Signal Processing, 2000, 47(10): 1008-1016.

[11].Chien-ChengTsengandSoo-ChangPei.'SparseFIRnotchfilterdesignanditsapplication'.IEE1997,ElectronicsLettersOnlineNO:19970766.[11]. Chien-Cheng Tseng and Soo-Chang Pei. 'Sparse FIR notch filter design and its application'. IEE1997, Electronics Letters Online NO: 19970766.

[12].陈彬.基于遗传算法的数字滤波器的实现：[硕士论文].南京林业大学,2009.[12]. Chen Bin. Realization of Digital Filter Based on Genetic Algorithm: [Master Thesis]. Nanjing Forestry University, 2009.

[13].刘欢.基于遗传算法的可重构FIR滤波器研究与设计[硕士论文].电子科技大学,2009.[13]. Liu Huan. Research and Design of Reconfigurable FIR Filter Based on Genetic Algorithm [Master Thesis]. University of Electronic Science and Technology of China, 2009.

[14].Zarahdnik,P.；Simak,B.；Vlcek,M.；,"FastDesignofPreciseEquirippleNotchFIRFilters,"Networking,InternationalConferenceonSystemsandInternationalConferenceonMobileCommunicationsandLearningTechnologies,2006.ICN/ICONS/MCL2006.InternationalConferenceon,vol.,no.,pp.154,23-29April2006.[14]. Zarahdnik, P.; Simak, B.; Vlcek, M.;, "Fast Design of Precise Equiripple NotchFIRFilters," Networking, International Conference on Systems and International Conference on Mobile Communications and Learning Technologies, 2006.ICN/ICONS/MCL2006.International Conference on, vol., no., pp.154, 23- 29 April 2006.

[15].EfronandTrevorHastieetal,‘Leastangleregression’,TheAnnalsofStatistics,2004,32:407-499[15]. Efron and Trevor Hastiie et al, 'Leastangle regression', The Annals of Statistics, 2004, 32:407-499

[16].张勇进.带弹性网的稀疏主成分分析：[硕士论文].华中科技大学.2006.[16]. Zhang Yongjin. Sparse Principal Component Analysis with Elastic Net: [Master Thesis]. Huazhong University of Science and Technology. 2006.

[17].柯郑林.Lasso及其相关方法在多元线性回归模型中的应用[博士论文].北京交通大学,2011[17]. Ke Zhenglin. Application of Lasso and related methods in multiple linear regression model [PhD dissertation]. Beijing Jiaotong University, 2011

[18].刘志平，石林英.最小二乘原理及其MATLAB实现.中国西部科技.2008.第七卷第17期总第142期.[18]. Liu Zhiping, Shi Linying. The Principle of Least Squares and Its Implementation in MATLAB. Science and Technology of Western China. 2008. Volume VII, Issue 17, Issue 142.

Claims

1. a kind of design method of sparse frequency-offset linear phase FIR notch wave filter, it is characterized in that the method can design the notch wave filter of low non-zero tap number, make it realize that the used adder multiplier number reduces, thereby Can improve its operation speed, reduce operation error and reduce energy consumption, the specific steps of the method include:

1. According to the design requirements, passband ripple δ, stopband bandwidth The notch frequency ω _s and the attenuation depth A _notch determine the initial order N of the linear phase FIR notch filter and the corresponding discretized ideal frequency response H _d (e ^jω ), where ω∈[0,π], j represents The imaginary number unit, the discretized ideal frequency response vector y represents H _d (e ^jω ) when ω takes L+1 discrete points {0,ω ₀ ,2ω ₀ ,3ω ₀ ,…,lω ₀ …,Lω ₀ } in turn value, where Indicates the discretized frequency interval, 0≤l≤L, L is a positive integer, and y is expressed as:

y the y = = {[[{H h}_{d d} (({e e}^{j j 00})),, {H h}_{d d} (({e e}^{{jω jω}_{00}})),, ... ... {H h}_{d d} (({e e}^{{jlω jlω}_{00}})) ... ...,, {H h}_{d d} (({e e}^{{jLω jLω}_{00}}))]]}^{T T} - - - - - - ((11))

The linear phase FIR notch filter tap coefficient vector β is expressed as

β=[β ₀ ,β ₁ ,...,β _m ,...β _M ] ^T ＝[h _M ,2h _M-1 ,...,2h _m ...,2h ₁ ,2h ₀ ] ^T (2)

Among them, M=N/2, h _m represents the mth tap coefficient of the FIR notch filter, 0≤m≤M; the sparse linear phase FIR notch filter design problem is transformed into the following mathematical optimization problem:

\hat{β} = \underset{β}{\arg m i no} | | β | |_{0},

and satisfied

\frac{| | x β - the y | |_{2}}{| | the y | |_{2}} < ϵ - - - (3)

Where ||·|| ₀ represents the 0-norm operation, ||β|| ₀ represents the number of non-zero taps in the tap coefficient vector, ||·|| ₂ represents the 2-norm operation; ε is the given error; the Vandermonde matrix X associated with the filter design is expressed as

where x _m represents the column vector in X;

2. Set the initial value of the column vector set corresponding to the non-zero coefficient to be an empty set, that is, X ₀ ={}, the initial estimated value of the FIR notch filter tap coefficient vector β is a vector of all 0s, the estimated value of the FIR notch filter frequency response vector μ Calculate the correlation coefficient Among them, 0≤m≤M, x _m is the column vector in the matrix X, find out the difference with the current residual The variable with the largest absolute value of the correlation coefficient and add this variable to the collection

x_{1} = [x_{j_{1}}];

3. In the kth iteration, use the LASSO algorithm and the vectors in the set X _k-1 to determine the minimum angle direction, and find a new column vector of X in X along this direction Will join the collection where 2≤k≤M;

Fourth, use the current set X _k and use the least squares method to estimate the FIR notch filter tap coefficient vector:

\overset{^^}{β β} = = {(({X x}_{k k} {X x}_{k k}^{T T}))}^{- - 11} {X x}_{k k}^{T T} y the y - - - - - - ((55))

5. Estimated value Transform to get the filter coefficients:

\overset{^^}{β β} = = {[[{\overset{^^}{β β}}_{00},, {\overset{^^}{β β}}_{11},, ... ...,, {\overset{^^}{β β}}_{m m},, ... ... {\overset{^^}{β β}}_{M m}]]}^{T T} = = {(({h h}_{M m},, 22 {h h}_{M m - - 11},, ... ...,, 22 {h h}_{00}))}^{T T} - - - - - - ((66))

6. For a given error ε, judge whether the following inequality is true,

\frac{| | | | X x \overset{^^}{β β} - - y the y | | {| |}_{22}}{| | | | y the y | | {| |}_{22}} < < ϵ ϵ - - - - - - ((77))

If the above formula holds, then That is, the specific value of the tap coefficient of the filter designed for this method, otherwise return to step 3, continue to increase the required set, and then determine the number of non-zero tap coefficients of the filter required by the sparse linear phase FIR notch filter , location, and specific values.