CN101190136B - 实时产生滤波器系数的方法和装置 - Google Patents

Abstract

本发明公开了一种实时产生滤波器系数的方法和装置。所述方法包括：基于当前的系数索引查询变换的窗函数表从而得到变换的窗函数值；产生当前的截止角频率；基于当前的系数索引、滤波器阶数和当前的截止角频率产生查表地址并基于该查表地址查询正弦表从而获得正弦值；将变换的窗函数值与正弦值相乘从而获得滤波器系数。所述装置包括存储有变换的窗函数表的第一存储器、存储有正弦表的第二存储器、查表地址产生模块和第一乘法器。本发明易于实现，消耗的硬件资源小，并且高度灵活。尤其适用于阶数很高的有限冲激响应滤波器的硬件实现。

Description

实时产生滤波器系数的方法和装置

技术领域

本发明涉及一种实时产生滤波器系数的方法和装置，并且更具体地说，涉及一种用于超声成像系统的实时产生滤波器系数的方法和装置。

背景技术

在超声成像系统中，前端发射波形进入人体器官后其幅度以及相位被调制，回波必须被解调后才能得到调制波的信息，解调是一个频移以及低通滤波的过程。超声波进入人体后，在不同的深度，随频率衰减的特性不同，越深的地方，高频成分衰减越多，因此回波的中心频率随深度逐渐减小，回波带宽亦随深度逐渐减小，因此滤波器的截止角频率亦随时间逐渐减小。滤波器的截止角频率决定滤波器的系数，因此滤波器的系数亦随时间相应调整。

图1给出典型的超声成像系统的信号处理流程。超声回波信号经过波束合成后形成射频回波信号。射频回波经过频移后得到基带成分以及倍频成分，经过低通滤波后仅得到基带成分，对回波数据解交织处理，得到同一位置点上不同时间的扫描数据，然后经过壁滤波，壁滤波是一个高通滤波器，去掉较强的低频组织运动成分，保留较弱的高频血流运动成分。对壁滤波后的数据进行自相关运算，得到血流的速度矢量在坐标轴上的分量以及血流能量信息，随后对这些血流参数进一步处理得到血流的速度、能量和方差信息，经过扫描变换后，将这些参数信息显示出来。

从原理上来说，解调的过程比较经典且简单，频移加滤波，只需要知道滤波器的系数就可以完成滤波的过程。得到滤波器系数的方法有很多，一是事先准备好系数，由系统软件将系数写入到可编程器件的片内随机存取存储器或者片外存储设备，当用到某一系数的时候再从片内随机存取存储器或片外存储设备中取出来，但是当滤波的阶数都是上百阶左右且都用较高的位数来量化时，如果事先将系数存储到片内随机存取存储器，对于随时间变化的滤波器来说存储量是无法容忍的。例如：假设滤波器的阶数为100，16比特量化，存储20组系数，总的存储量为：100*16*20＝32K比特。如果存储到片外存储设备，存储量可以接受，但却要花费额外的代价来添加片外存储设备，因此对于滤波器来说并不可取。

在专利6,142,942(U.S.Pat.No.6,142,942)中，Clark介绍了一种自适应滤波器系数产生的方法。该方法不需要存储大量的系数，系数是实时产生的。利用回波数据的统计特性，在不同的时间段产生系数，然后对数据进行滤波运算，将滤波后的结果作进一步处理。

如上所述，在现有的技术解决方案中，有两种产生滤波器系数的方法，一是所有系数预先产生，放在片内或者片外存储设备中，当用到的时候再从存储设备中取出，对于滤波器阶数很小而且存储组数少的时候，该方法可行。但在超声成像系统中，滤波器阶数很高，通常有上百阶，而且要存储几十组不同的系数，会占用大量的存储资源。对于自适应滤波器来说，不需要存储系数，系数由不同时间段数据的统计特性实时产生。但该方法同样要对数据作大量的统计信号处理，算法复杂繁琐。对于有限资源的可编程器件来说同样不可取。

发明内容

本发明的目的是根据变化的截止角频率，实时计算产生每一时间段的滤波器系数。由于只要知道滤波器的截止角频率，就可以得到滤波器的系数，一个截止角频率值对应着一组系数，存储截止角频率消耗的片内随机存取存储器比存储系数消耗的片内随机存取存储器要少得多。

根据本发明的第一方面，提供一种实时产生滤波器系数的方法，所述方法包括：第一查表步骤，用于基于当前的系数索引对存储有变换的窗函数表的第一存储器查表从而获得第一中间值；查表地址产生步骤，用于基于当前的系数索引、滤波器阶数和当前的截止角频率产生查表地址；第二查表步骤，用于基于查表地址对存储有正弦表的第二存储器查表从而获得第二中间值；以及第一乘法步骤，用于将第一中间值与第二中间值相乘从而获得滤波器系数。

根据本发明第一方面的实时产生滤波器系数的方法，其中查表地址产生步骤包括：索引变换步骤，用于对当前的系数索引进行索引变换从而获得第三中间值；截止角频率产生步骤，用于产生当前的截止角频率；以及第二乘法步骤，用于将第三中间值与当前的截止角频率相乘从而获得查表地址。

根据本发明第一方面的实时产生滤波器系数的方法，其中索引变换步骤中的索引变换为：

T_n＝2n+1-N，

其中T_n为第三中间值，n为当前的系数索引，n＝0，1，....N-1，N为滤波器阶数。

根据本发明第一方面的实时产生滤波器系数的方法，其中截止角频率产生步骤包括：第三查表步骤，用于基于计数器产生的当前的时间段序号对存储有截止角频率变化量表的第三存储器查表从而获得当前的截止角频率变化量；以及加法步骤，用于将存储在寄存器中的上一个截止角频率与当前的截止角频率变化量相加从而获得当前的截止角频率。

根据本发明第一方面的实时产生滤波器系数的方法，其中第一乘法步骤中的滤波器系数的计算公式如下：

$h\left(n\right)=\sin [(2n+1-N)*k*\mathrm{pi}/M]\·\frac{(0.54-0.46\cos (2*\mathrm{pi}*n/(N-1)))}{\mathrm{pi}*(n-(N-1)/2)}$

其中，h(n)为滤波器系数，N为滤波器阶数，n为当前的系数索引，n＝0，1，....N-1，k为当前的截止角频率，pi为圆周率，M为2的整次幂。

根据根据本发明的第二方面，提供一种实时产生滤波器系数的装置，所述装置包括：第一存储器，用于存储变换的窗函数表，以便用当前的系数索引对第一存储器查表从而获得第一中间值；查表地址产生模块，用于基于当前的系数索引、滤波器阶数和当前的截止角频率产生查表地址；第二存储器，用于存储正弦表，以便用查表地址对第二存储器查表从而获得第二中间值；以及第一乘法器，用于将第一中间值与第二中间值相乘从而获得滤波器系数。

根据本发明第二方面的实时产生滤波器系数的装置，其中查表地址产生模块包括：索引变换单元，用于对当前的系数索引进行索引变换从而获得第三中间值；截止角频率产生单元，用于产生当前的截止角频率；以及第二乘法器，用于将第三中间值与当前的截止角频率相乘从而获得查表地址。

根据本发明第二方面的实时产生滤波器系数的装置，其中索引变换单元进行的索引变换为：

T_n＝2n+1-N，

其中T_n为第三中间值，n为当前的系数索引，n＝0，1，....，N-1，N为滤波器阶数。

根据本发明第二方面的实时产生滤波器系数的装置，其中截止角频率产生单元包括：计数器，用于产生当前的时间段序号；第三存储器，用于存储截止角频率变化量表，以便用当前的时间段序号对第三存储器查表从而获得当前的截止角频率变化量；寄存器，用于存储上一个截止角频率；以及加法器，用于将上一个截止角频率与当前的截止角频率变化量相加从而获得当前的截止角频率。

根据本发明第二方面的实时产生滤波器系数的装置，其中计数器为均匀计数器。

根据本发明第二方面的实时产生滤波器系数的装置，其中第一存储器和第三存储器为随机存取存储器，第二存储器为只读存储器。

根据本发明第二方面的实时产生滤波器系数的装置，其中第一乘法器计算滤波器系数的公式如下：

$h\left(n\right)=\sin [(2n+1-N)*k*\mathrm{pi}/M]\·\frac{(0.54-0.46\cos (2*\mathrm{pi}*n/(N-1)))}{\mathrm{pi}*(n-(N-1)/2)}$

其中，h(n)为滤波器系数，N为滤波器阶数，n为当前的系数索引，n＝0，1，....，N-1，k为当前的截止角频率，pi为圆周率，M为2的整次幂。

本发明实现简单，消耗的可编程硬件资源极小，并且高度灵活。通过调整正弦表以及变换的窗函数表里面存储的值，可以配置为各种不同选通类型的滤波器，也可以配置为不同窗函数的滤波器，还可以配置为截止角频率任意变化的滤波器。尤其适用于阶数很高的有限冲激响应滤波器的硬件实现。

附图说明

通过结合以下附图，并参考以下对具体实施方式的详细说明，可以对本发明有更透彻的理解：

图1为典型的超声成像系统的信号处理流程图；

图2为根据本发明的实时产生滤波器系数的流程图；

图3为写变换的窗函数表的时序图；

图4为读变换的窗函数表的时序图；

图5为一线时间等分序号图；

图6为根据本发明的产生当前的截止角频率的时序图；

图7为根据本发明的产生查表地址以及查表时序图；

图8为输出滤波器系数的时序图；

图9为根据本发明的实时产生滤波器系数的装置的结构框图；以及

图10为根据本发明的查表地址产生模块的结构框图。

具体实施方式

在本发明的第一实施例中，提供一种用于实时产生滤波器系数的方法。所述方法包括：基于当前的系数索引对存储有变换的窗函数表的第一存储器查表从而获得第一中间值；基于当前的系数索引、滤波器阶数和当前的截止角频率产生查表地址；基于查表地址对存储有正弦表的第二存储器查表从而获得第二中间值；以及将第一中间值与第二中间值相乘从而获得滤波器系数。下面将具体说明实时产生滤波器系数的方法。

首先说明滤波器的几个基本概念。设计滤波器的一个最常见的方法是理想滤波器加窗的方法，根据滤波器的截止角频率，得到理想滤波器的系数，这样得到的滤波器是非因果且不可实现，然后将理想滤波器的系数和窗函数相乘，移位后得到因果可实现的滤波器系数。在以下对本发明实施例的详细说明中，为便于解释和说明，假定所加的窗函数为汉明窗，并以实时产生低通滤波器系数为例。但是这并不意味着本发明仅限于此。

假设低通滤波器的阶数为N(N为偶数)，则基于理想滤波器加窗的方法可得低通滤波器的系数为：

$h\left(n^{\′}\right)=\frac{\sin \left(n^{\′}{\mathrm{\ω}}_c\right)}{n^{\′}{\mathrm{\ω}}_c}\·(0.54+0.46\cos (2*\mathrm{pi}*n^{\′}/N))---\left(1\right)$

其中 $n^{\′}=\mathrm{round}(-\frac{N}{2}+1),......,\mathrm{round}\left(\frac{N}{2}\right)$

令 $n=n^{\′}+\mathrm{floor}\left(\frac{N}{2}\right)-1,$ 移位后则

$h\left(n\right)=\sin \left[(n-\frac{N-1}{2}){\mathrm{\ω}}_c\right]\·\frac{(0.54-0.46\cos (2*\mathrm{pi}*n/(N-1)))}{\mathrm{pi}*(n-(N-1)/2)}---\left(2\right)$

其中ω_c为低通滤波器的截止角频率，n＝0，1，....，N-1，为低通滤波器的系数索引。于是低通滤波器的系数由两部分相乘得到，一部分仅与正弦值有关，另外一部分和汉明窗有关。

假如ω_c＝k*2pi/M，则k为ω_c对k*2pi/M量化后的值，k＝0～M-1，pi为圆周率，M取2的整次幂，如256，则上式(2)中乘号左边的项为：

$\sin \left[(n-\frac{N-1}{2}){\mathrm{\ω}}_c\right]=\sin [(2n+1-N)*k*\mathrm{pi}/M]---\left(3\right)$

该值可以通过查找一个长度为M的正弦表得到，该长度为M的正弦表存储由0到pi均匀分成M个角度的正弦值。式(2)乘号右边的项为变换的窗函数，可以预先计算得出并存储在随机存取存储器中，也可基于系数索引n查表得到对应的值。因此主要通过两次查表和一次乘法就可以得到低通滤波器的一个系数。

图2给出实时产生滤波器系数的流程图。基于当前的系数索引查询变换的窗函数表从而得到变换的窗函数值(步骤202)；产生当前的截止角频率(步骤204)；基于当前的系数索引、滤波器阶数和当前的截止角频率产生查表地址并基于该查表地址查询正弦表从而获得正弦值(步骤206)；将变换的窗函数值与正弦值相乘从而获得滤波器系数(步骤208)。以下对各个步骤进行详细说明。

对于N阶低通滤波器来说，其系数索引为n，n＝0，1，....，N-1。在步骤202中，用该索引n直接查询变换的窗函数表。从式(2)可以看出，该表不是标准的汉明窗，而是除以一个值，得到变换的窗函数表。将变换的窗函数表存储在双口随机存取存储器中，该双口随机存取存储器一端只写，另外一端只读，其深度为N。对于不同阶数的滤波器来说，其初始值不一样，系统在初始化的时候，由软件接口写入，每次滤波器的阶数需要变化时，表里面的内容被更新。其值量化的精度视要求而定，通常为16比特有符号数。在滤波器的设计中，为了得到更好的性能，可以改变所加窗的类型，比如布莱克曼窗或者汉宁窗等。图3为写变换的窗函数表的时序图，图4为读变换的窗函数表的时序图。

对于超声系统黑白信号处理来说，随着探测深度的增加，其中心频率向低频偏移，带宽降低，因此设计滤波器的时候，随着时间的变化，滤波器的截止角频率逐渐降低。在步骤204中，产生当前的截止角频率，具体过程如下。将截止角频率变化量量化后，存储于片内随机存取存储器中，位宽用3比特补码表示，既可递增，也可递减，还可以不变，视其内容而定。该随机存取存储器容量较小，由系统在初始化的时候写入。将一线的数据在时间上分为均匀的若干段，由计数器来产生均匀的时间段，用时间段序号查询截止角频率变化量表，用查询出来的值和存储在寄存器中的上一个截止角频率相加，即为当前时间段的截止角频率k。图5为一线时间等分序号图，图6为产生当前的截止角频率的时序图。

在步骤206中，基于系数索引n、滤波器阶数N和当前的截止角频率k产生查表地址并查询正弦表从而获得正弦值，具体过程如下。在查询正弦表之前，需要完成索引变换，即由n、N变为(2n+1-N)得到中间值T_n。然后计算T_n*k得到查表地址后直接查询正弦表。正弦表的深度和截止角频率量化的精度有关，比如将0至2pi区间的角度值256等分，则式(3)中的截止角频率k的取值范围为0至127，正弦表的深度为256，其值通常用16比特量化。考虑到正弦值的对称性，只需存储1/4象限的值，因此存储量又可以降为1/4，乘完后的结果取低8位，用低6位直接查询，高两位作为符号、旋转判断的依据。通常正弦值较为固定，无需用软件配置写入。可以指定初始值后做成只读存储器，仅仅当截止角频率的量化精度或正弦值的量化精度改变时，该只读存储器的内容需重新初始化。图7为产生查表地址以及查表时序图。

在步骤208中，将变换的窗函数值与正弦值相乘从而获得滤波器系数。在一些实施例中，对于同一个系数索引，查询变换的窗函数表后得出变换的窗函数值，查询正弦表得出正弦值，两者的延时量不同，得到正弦值的中间计算步骤多，因此延时量更多，需要将查询出的变换的窗函数值同步延时，在时序上保证经由同一个系数索引查询出的正弦值和变换的窗函数值对齐。查询的结果均为补码输出，然后相乘，相乘结果截取高位后作为系数输出。图8为输出滤波器系数的时序图。

本实施例提供的实时产生低通滤波器系数的方法可以容易地扩展为实时产生高通、带通或带阻滤波器系数的方法。

在本发明的第二实施例中，提供一种用于实时产生滤波器系数的装置，以下通过图9和图10进行详细说明。

图9为根据本发明的实时产生滤波器系数的装置的结构框图。如图9所示，该装置包括：第一存储器902，第二存储器904，查表地址产生模块906以及第一乘法器908。第一存储器902为随机存取存储器，用于存储变换的窗函数表。第二存储器904为只读存储器，用于存储正弦表。用当前的系数索引n对第一存储器902查表从而得到变换的窗函数值H_a。查表地址产生模块906用于根据当前的系数索引n、滤波器阶数N和当前的截止角频率k_D产生查表地址，下面将对其进行进一步说明。利用由查表地址产生模块906产生的查表地址对第二存储器904查表从而获得正弦值S_n。用第一乘法器908将变换的窗函数值H_a与正弦值S_n相乘从而获得滤波器系数。

图10为根据本发明的查表地址产生模块1006的结构框图。查表地址产生模块1006包括索引变换单元1010，用于对当前的系数索引n进行索引变换从而获得中间值T_n，即T_n＝2n+1-N；截止角频率产生单元1012，用于产生当前的截止角频率k_D；以及第二乘法器1014，用于将当前的截止角频率k_D与中间值T_n相乘从而获得查表地址。如图所示，该截止角频率产生单元1012包括计数器1016、第三存储器1018、寄存器1020和加法器1022。第三存储器1018为随机存取存储器，用于存储截止角频率变化量表。寄存器1020用于存储上一个截止角频率k_D-1。计数器1016为均匀计数器，用于产生当前的时间段序号，以便用该时间段序号对第三存储器1018查表从而获得当前的截止角频率变化量A_c。加法器1022将当前的截止角频率变化量A_c与上一个截止角频率k_D-1相加从而获得当前的截止角频率k_D。

在本实施例中，上一个截止角频率k_D-1、当前的截止角频率变化量A_c和当前的截止角频率k_D均为对2*pi/M量化后的值，其中pi为圆周率，M为2的整次幂，如256。

从上面可以看出，本发明的存储量非常小，包括截止角频率变化量表、1/4周期的正弦表和变换的窗函数表，用ALTERA的可编程器件实现时，6个M4K的存储单元即可。具体实施中仅用到两个乘法器，用可编程器件的专用硬件乘法器即可实现，无需逻辑单元实现乘法器的功能。具体实施中用逻辑单元来实现下列功能：均匀计数器、加法器和寄存器，仅用到200个逻辑单元，消耗的资源非常少。考虑到可以将滤波器设置为其它的类型，比如高通滤波器、带通滤波器和带阻滤波器，只需要更改正弦表的值即可。如果将滤波器配置为截止角频率不随时间变化的滤波器，仅需将截止角频率变化量表的内容置为0。如果需要其它类型的窗函数，比如巴特利特窗、矩形窗等等，只需要更改变换的窗函数表的值。对于需要并行输出多个滤波器系数的情况，可以用复制硬件，并行查询，也可以用高倍时钟复用系数索引后查询输出，然后将查询的结果串并变换。本发明适合于阶数很高比如上百阶以上、同时系数需要随时间变化的滤波器设计。

本发明不仅如上所述可应用于超声成像系统中，还可以应用于雷达信号处理系统和通信系统中。

以上通过特定的实施例对本发明进行了详细的描述，但本发明并不限于上述实施例。在不脱离本发明范围的前提下，可以对本发明进行各种修改和变更。本发明的范围由所附权利要求书限定。在本发明的说明书和权利要求书中所用的诸如“第一”、“第二”和“第三”等词语仅仅是为了便于描述，而不具有任何限制意义。

Claims (11)

1.一种用于超声成像系统的实时产生滤波器系数的方法，其特征在于，所述方法包括：

第一查表步骤，用于基于当前的系数索引对存储有变换的窗函数表的第一存储器查表从而获得第一中间值；

查表地址产生步骤，用于基于所述当前的系数索引、滤波器阶数和当前的截止角频率产生查表地址；

第二查表步骤，用于基于所述查表地址对存储有正弦表的第二存储器查表从而获得第二中间值；以及

第一乘法步骤，用于将所述第一中间值与所述第二中间值相乘从而获得所述滤波器系数，

其中，所述滤波器系数用于对所述超声成像系统中的回波数据进行滤波。

2.如权利要求1所述的方法，其特征在于，所述查表地址产生步骤包括：

索引变换步骤，用于对所述当前的系数索引进行索引变换从而获得第三中间值；

截止角频率产生步骤，用于产生所述当前的截止角频率；以及

第二乘法步骤，用于将所述第三中间值与所述当前的截止角频率相乘从而获得所述查表地址，

其特征在于，所述索引变换步骤中的索引变换为：

T_n＝2n+1-N，

其中T_n为所述第三中间值，n为所述当前的系数索引，n＝0，1，....，N-1，N为所述滤波器阶数。

3.如权利要求2所述的方法，其特征在于，所述截止角频率产生步骤包括：

第三查表步骤，用于基于计数器产生的当前的时间段序号对存储有截止角频率变化量表的第三存储器查表从而获得当前的截止角频率变化量；以及

加法步骤，用于将存储在寄存器中的上一个截止角频率与所述当前的截止角频率变化量相加从而获得所述当前的截止角频率。

4.如权利要求1到3中任一项所述的方法，其特征在于，所述第一乘法步骤中的滤波器系数的计算公式如下：

$h\left(n\right)=\sin [(2n+1-N)*k*\mathrm{pi}/M]\·\frac{(0.54-0.46\cos (2*\mathrm{pi}*n/(N-1)))}{\mathrm{pi}*(n-(N-1)/2)}$

其中，h(n)为所述滤波器系数，n为所述当前的系数索引，n＝0，1，....，N-1，N为所述滤波器阶数，k为所述当前的截止角频率，pi为圆周率，M为2的整次幂。

5.一种实时产生滤波器系数的装置，其特征在于，所述装置包括：

第一存储器，用于存储变换的窗函数表，以便用当前的系数索引对所述第一存储器查表从而获得第一中间值；

查表地址产生模块，用于基于所述当前的系数索引、滤波器阶数和当前的截止角频率产生查表地址；

第二存储器，用于存储正弦表，以便用所述查表地址对所述第二存储器查表从而获得第二中间值；以及

第一乘法器，用于将所述第一中间值与所述第二中间值相乘从而获得所述滤波器系数。

6.如权利要求5所述的装置，其特征在于，所述查表地址产生模块包括：

索引变换单元，用于对所述当前的系数索引进行索引变换从而获得第三中间值；

截止角频率产生单元，用于产生所述当前的截止角频率；以及

第二乘法器，用于将所述第三中间值与所述当前的截止角频率相乘从而获得所述查表地址，

其特征在于，所述索引变换单元进行的索引变换为：

T_n＝2n+1-N，

其中T_n为所述第三中间值，n为所述当前的系数索引，n＝0，1，....，N-1，N为所述滤波器阶数。

7.如权利要求6所述的装置，其特征在于，所述截止角频率产生单元包括：

计数器，用于产生当前的时间段序号；

第三存储器，用于存储截止角频率变化量表，以便用所述当前的时间段序号对所述第三存储器查表从而获得当前的截止角频率变化量；

寄存器，用于存储上一个截止角频率；以及

加法器，用于将所述上一个截止角频率与所述当前的截止角频率变化量相加从而获得所述当前的截止角频率。

8.如权利要求7所述的装置，其特征在于，所述计数器是均匀计数器。

9.如权利要求5到7中任一项所述的装置，其特征在于，所述第一存储器为随机存取存储器，所述第二存储器为只读存储器。

10.如权利要求7所述的装置，其特征在于，所述第三存储器为随机存取存储器。

11.如权利要求5到7中任一项所述的装置，其特征在于，所述第一乘法器计算滤波器系数的公式如下：

$h\left(n\right)=\sin [(2n+1-N)*k*\mathrm{pi}/M]\·\frac{(0.54-0.46\cos (2*\mathrm{pi}*n/(N-1)))}{\mathrm{pi}*(n-(N-1)/2)}$

其中，h(n)为所述滤波器系数，N为所述滤波器阶数，n为所述当前的系数索引，n＝0，1，....，N-1，k为所述当前的截止角频率，pi为圆周率，M为2的整次幂。

Images