CN107294511A - 一种低复杂度的可变分数时延滤波方法及滤波器 - Google Patents

一种低复杂度的可变分数时延滤波方法及滤波器 Download PDF

Info

Publication number
CN107294511A
CN107294511A CN201710385496.2A CN201710385496A CN107294511A CN 107294511 A CN107294511 A CN 107294511A CN 201710385496 A CN201710385496 A CN 201710385496A CN 107294511 A CN107294511 A CN 107294511A
Authority
CN
China
Prior art keywords
mrow
mtd
filter
mtr
msub
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201710385496.2A
Other languages
English (en)
Other versions
CN107294511B (zh
Inventor
黄翔东
张博
马欣
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tianjin University
Original Assignee
Tianjin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tianjin University filed Critical Tianjin University
Priority to CN201710385496.2A priority Critical patent/CN107294511B/zh
Publication of CN107294511A publication Critical patent/CN107294511A/zh
Application granted granted Critical
Publication of CN107294511B publication Critical patent/CN107294511B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • H—ELECTRICITY
    • H03—ELECTRONIC CIRCUITRY
    • H03H—IMPEDANCE NETWORKS, e.g. RESONANT CIRCUITS; RESONATORS
    • H03H17/00—Networks using digital techniques
    • H03H17/0009—Time-delay networks
    • H03H17/0018—Realizing a fractional delay

Landscapes

  • Complex Calculations (AREA)

Abstract

本发明公开了一种低复杂度的可变分数时延滤波方法及滤波器，包括：根据频率采样向量长度、截止频率确定滤波器过渡带和边界整数；通过归一化的卷积窗、频率采样向量长度、边界整数获取全相位滤波器；对全相位滤波器进行三次样条插值，获取样条系数；将分数时延值的不同取值p∈[0,0.5]或p∈[‑0.5,0)、以及样条系数带入泰勒解析式中求得可变分数时延滤波器法罗结构的子滤波器系数。本方法将全相位滤波器的设计、三次样条插值与泰勒系数展开进行了有机的综合，实现了可变分数时延滤波器Farrow结构抽头系数的快速配置。

Description

一种低复杂度的可变分数时延滤波方法及滤波器
技术领域
本发明涉及数字信号处理技术领域，尤其涉及一种低复杂度的可变分数时延滤波方法及滤波器。
背景技术
可变分数时延(Variable fractional delay，VFD)数字滤波器广泛应用于采样率转换[1]、时延估计[2]、梳状滤波器设计[3]、离散时间信号插值[4]、图像插值[5]与雷达波束形成[6]等场合。文献[7]提出一种用于VFD滤波器设计的Farrow(法罗)结构，该结构由多个直接型有限冲击响应(Finite Impulse Response，FIR)子滤波器组成，只需用不同幂次方的时延参数对这些FIR滤波器做线性组合即可灵活调节滤波器的群延时，因而20多年来，Farrow结构成为VFD滤波器的主流框架[8-14]。相应地，VFD滤波器的设计，实际上就是对Farrow结构多个子滤波器的系数进行配置的问题。实际应用中，总是期望该系数配置以高效、解析的方式实现，然而现有的两种主流方法(极小极大方法[9-11]与加权最小二乘法(Weighted Least Squares，WLS)[12-14])很难满足这个需求。
极小极大方法(Minimax)通过最小化滤波器的当前传输特性与期望传输特性之间的最大幅度误差来获得滤波器系数。文献[9]指出，该方法要求对所有子滤波器的系数分别优化，并且每次优化都涉及多次迭代更新，因此运算复杂度很高；此外，每个子滤波器的期望特性都源于理想滤波器的泰勒展开对应项的传输特性，文献[13]指出，由于理想滤波器的冲激响应衰减非常缓慢，这必然导致VFD滤波器设计需耗费大量的子滤波器，从而进一步增加了运算复杂度。因此，相较于优化设计方法，解析设计法更受欢迎。
与极小极大方法不同，WLS方法通过求解线性方程获得滤波器系数[14](该线性方程由滤波器传输特性与理想滤波器传输特性的加权误差的最小二乘准则导出)，因此该方法是一种解析设计方法。然而，该线性方程的求解包含矩阵求逆，故其运算复杂度会随着滤波器阶数的增加而急速增加。
发明内容
本发明提供了一种低复杂度的可变分数时延滤波方法及滤波器，本发明将全相位滤波器的设计、三次样条插值与泰勒系数展开进行了有机的综合，实现了可变分数时延滤波器Farrow结构抽头系数的快速配置，详见下文描述：
一种低复杂度的可变分数时延滤波方法，所述方法包括以下步骤：
根据频率采样向量长度、截止频率确定滤波器过渡带和边界整数；
通过归一化的卷积窗、频率采样向量长度、边界整数获取全相位滤波器；
对全相位滤波器进行三次样条插值，获取样条系数；
将分数时延值的不同取值p∈[0,0.5]或p∈[-0.5,0)、以及样条系数带入泰勒解析式中求得可变分数时延滤波器法罗结构的子滤波器系数。
其中，所述通过归一化的卷积窗、频率采样向量长度、边界整数获取全相位滤波器具体为：
其中，g(n)为全相位滤波器；wc(n)为归一化的卷积窗；N为频率采样向量长度；K为边界整数。
其中，所述边界整数滤波器过渡带Δω＝2π/N；ωc为截止频率。
进一步地，所述对全相位滤波器进行三次样条插值，获取样条系数具体为：
其中，bn,0、bn,1、bn,2、bn,3分别为样条系数；g(n)为全相位滤波器；g(n+1)为全相位滤波器g(n)以1为单位的时域延拓；g″n与g″n+1分别为在n和n+1两个位置算出数值的二阶导数。
进一步地，当p∈[0,0.5]时，子滤波器系数为：
其中，a(n,m)为子滤波器系数；bn,i为样条系数。
进一步地，当p∈[-0.5,0)时，子滤波器系数为：
其中，a(n,m)为子滤波器系数；bn,i为样条系数。
一种低复杂度的可变分数时延滤波器，所述滤波器包括：外部RAM、数字信号处理器、以及输出驱动及显示电路，
所述外部RAM用于接收滤波器分数时延、截止频率、频率采样向量长度、即归一化的卷积窗；
所述数字信号处理器经过处理得到全相位滤波器及其传输曲线；利用三次样条插值及泰勒系数展开快速计算得到可变分数时延滤波器法罗结构的子滤波器系数及其传输曲线；所述输出驱动及显示电路用于输出显示。
进一步地，所述数字信号处理器为DSP。
本发明提供的技术方案的有益效果是：
1、本方法将全相位滤波器的设计、三次样条插值与泰勒系数展开进行了有机的综合，实现了可变分数时延滤波器Farrow结构抽头系数的快速配置；
2、理论分析和仿真实验都可以证明，本发明设计出的可变分数时延滤波器具有高灵活性(对群延时和截止频率)、高精度与低复杂度。
附图说明
图1为可变分数时延滤波器的Farrow结构的示意图；
图2为全相位滤波器的传输曲线(N＝10，K＝3)的示意图；
图3为三种冲激响应关系的示意图；
图4为VFD滤波器的传输曲线的示意图；
图5为VFD滤波器的群延时曲线的示意图；
图6为两种设计在不同滤波器阶数时的执行时间的示意图；
图7为本发明的硬件实施图；
图8为DSP内部程序流图；
表1为三步解析设计流程；
表2为两方法的群延时误差(p＝0.3,N＝18)。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。
为了解决现有技术中存在的问题，有必要提出一种低复杂度的可变分数时延滤波方法及滤波器。本发明实施例将全相位滤波器解析设计、三次样条插值和泰勒级数展开做有机结合，推导出Farrow结构子滤波器的解析配置公式，实现了低复杂度VFD滤波器的高效、解析设计。
实施例1
101：根据频率采样向量长度、截止频率确定滤波器过渡带和边界整数；通过归一化的卷积窗、频率采样向量长度、边界整数获取全相位滤波器；
102：对全相位滤波器进行三次样条插值，获取样条系数；
103：将分数时延值的不同取值(p∈[0,0.5]或p∈[-0.5,0))、以及样条系数bn,m带入泰勒解析式中求得VFD滤波器Farrow结构的子滤波器系数。
综上所述，本发明实施例通过上述步骤101-步骤103实现了可变分数时延滤波器Farrow结构抽头系数的快速配置，具有高灵活性(对群延时和截止频率)、高精度与低复杂度。
实施例2
下面结合具体的计算公式、实例、表1对实施例1中的方案进行进一步地介绍，详见下文描述：
201：给定工程需求的分数时延p、频率采样向量长度N、截止频率ωc和一个归一化的卷积窗wc(n)，确定滤波器的滤波器过渡带Δω＝2π/N和边界整数
202：根据解析式
计算得到全相位滤波器g(n),-N+1≤n≤N-1。
203：对全相位滤波器g(n)进行三次样条插值，根据解析式计算得到样条系数bn,m(-N+1≤n≤N-1,0≤m≤3)；
其中，g(n+1)为上一步推导出的全相位滤波器解析式；g″n+1与g″n分别为在n和n+1两个位置算出的g(t)的数值二阶导数。
204：根据p值的不同情况(p∈[0,0.5]或p∈[-0.5,0))，将bn,m代入下式求得VFD滤波器Farrow结构的子滤波器系数a(n,m)，(-N+2≤n≤N-2,0≤m≤3)，即：
一、Farrow滤波结构
对于以Farrow结构作为框架的VFD滤波器，其传输函数可表示为[7]
其中，分数p(p∈[-0.5,0.5])是用于调节群延时的参数，hn(p)为VFD滤波器的系数，Z为滤波器的Z变换，即通过滤波器系数计算得到滤波器传输函数。
而hn(p)可进一步表示为关于p的多项式：
其中，M为Farrow结构中子滤波器的个数。联立式(1)、式(2)，有
其中，Gm(z)为每个子滤波器的传输函数。
因此，H(z,p)可由图1所示的Farrow结构实现。相应地，VFD滤波器的设计，即转化为图1中子滤波器系数a(n,m)的合理配置问题。
二、VFD滤波器的解析设计
本发明实施例提出的VFD滤波器设计包含3个部分(即全相位滤波器设计、三次样条插值和泰勒级数展开)，这3部分设计都可解析实现，故本发明设计的滤波器的计算复杂度低。
1、全相位滤波器的解析设计
如文献[15,16]中指出，给定频率采样向量H和一个长度为2N-1的窗wc(n)，其中H为：
wc(n)由一个长度为N的常用窗f(n)(归一化因子为与一个长度为N的反转矩形窗(RN(-n))卷积得到，即：
wc(n)＝f(n)*RN(-n),n＝-N+1,...,N-1 (5)
需将H定义域延拓后的IDFT结果h(n)与wc(n)相乘，便可得到全相位滤波器系数g(n)的解析表达式为：
文献[18]中证明，滤波器g(n)的传输特性G(jω)可通过对频率向量H做内插得到，并且其内插函数恰好是归一化卷积窗wc(n)的傅里叶频谱Wc(jω)，因此G(jω)等于：
其中，H(k)为频率采样向量；ω为滤波器的数字角频率。
众所周知，任意常用窗(例如：汉宁窗、布莱克曼窗)的傅里叶频谱F(jω)可以表示如下：
其中，Fg(ω)为常用窗的幅频响应。
同样地，矩形窗RN(n)的傅里叶频谱RN(jω)可以表示为：
结合式(5)、(8)与(9)，可以得到
其中，RN(-jω)为反转矩形窗(RN(-n))的傅里叶频谱；j0根据公式(10)计算得到的，相位相消后为0。
公式(10)表明Wc(jω)在整个频域是零相位的，因此式(7)中的传输特性G(jω)也是零相位的，即全相位滤波器g(n),n＝-N+1,...,N-1的群延时等于0。
此外，文献[18]证明了卷积窗频谱Wc(jω)具有理想的频率采样特性，传输曲线G(jω)恰好通过H中的N个频率采样点，即：
G(jkΔω)＝H(k),k＝0,...,N-1 (11)
其中，G(jkΔω)即为传输曲线G(jω)，其中ω＝kΔω。
设定N＝10，K＝3，公式(6)中的wc(n)为汉宁单窗，将这些参数代入式(6)即可得到全相位滤波器g(n)，其幅度曲线如图2所示。
从图2中可看出，由于G(jω)通过所有的频率设置点，很好的逼近了期望的传输曲线Gd(jω)，故其在过渡带ω∈[K-1,K]Δω上近似为直线，从而3dB截至频率ωc近似满足：
反之，若给定参数N与截止频率ωc，则参数K可以计算如下：
这里[·]表示四舍五入为整数。
需要强调，因为定义域n∈[-N+1,N-1]，所以式(6)中全相位滤波器g(n)的群延时恒等于0。为获得分数时延，需要将离散定义域n∈[-N+1,N-1]延拓为连续定义域t∈[-N+1,N-1]，对g(n)施加差值措施即可实现这个扩展。
2、基于三次样条插值的定义域延拓
出于降低计算复杂度和节省硬件资源的考虑，将离散冲激响应g(n)转换为连续函数g(t)的插值函数不需要具有很高的阶数，因此本发明实施例引入三次样条插值(M＝3)对全相位滤波器系数g(n)做定义域延拓。
对于2N-1个全相位滤波器系数g(-N+1),...,g(N-1)，存在2N-2个单位区间t∈[n,n+1],n＝-N+1,...,N-2。故将取值离散的g(n)插值延拓为取值连续的函数g(t)后，该连续函数可表示为2N-2个分段函数gn(t)的逼近，即：
其中，gn(t)是如下所示的三次多项式：
根据三次样条插值的定义[20]，式(15)中的4个参数可分别计算如下：
其中，g″n与g″n+1分别为在n和n+1两个位置算出的数值二阶导数。
3、Farrow结构的子滤波器系数配置
为得到图1所示的Farrow结构子滤波器系数a(n,m)，需对g(t)进行泰勒级数展开。具体而言，除去t＝-N+1和t＝N-1这两个边界点，泰勒展开的固定点可选为剩下的2N-3个全相位滤波器离散取值点t＝n,n∈[-N+2,N-2]。不失一般性，令t＝n+p，分为两种情况进行讨论。
1)第一情况：0≤p≤0.5；
该情况下，t落在区间[n,n+1)内，因此g(t)可以被gn(t)的泰勒系数展开式逼近：
结合式(2)与(17)，可推导出Farrow结构的子滤波器系数a(n,m)具有如下形式：
根据式(15)，式(18)中的各阶导数可表示为：
进一步结合式(18)与(19)，可归纳得到子滤波器系数的统一公式为：
2)第二情况：-0.5≤p＜0；
该情况下，t落在区间[n-1,n)内，因此g(t)可以被gn-1(t)的泰勒系数展开式逼近：
经过类似于第一情况(0≤p≤0.5)的推导，该第二情况的子滤波器系数的统一公式被归纳为：
三、三步解析设计
总结以上3部分，本发明实施例提出的VFD滤波器可归纳为表1中如下3个步骤：
表1三步解析设计流程
从表1看出，本发明实施例仅需4个简单解析式即可配置出子滤波器系数，所有步骤都不涉及任何系数优化迭代和矩阵求逆操作，因此运算复杂度非常低。
四、性能分析
1)VFD滤波器与全相位滤波器的关系
众所周知，三次样条插值能够保证在任何相邻区间内曲线连续。故在区间t∈[n-1,n)和t∈[n,n+1)的连接点t＝n处，插值多项式gn(t)、gn-1(t)应满足：
结合式(23)与(17)、(21)可知，截掉两个边界系数g(-N+1),g(N-1)，全相位滤波器g(n)即对应图1中Farrow结构的子滤波器G0(z)。
用表示本发明实施例设计出的VFD滤波器，从前述的定义域拓展可推知：g(n)与都可通过对连续冲激响应g(t)做等间隔采样得到，它们的采样点分别是t＝n与t＝n+p(如图3所示)，由此可得：
因此，VFD滤波器近似为全相位滤波器g(n)分数时延后的版本。
2)截止频率与群延时的可变性
根据傅里叶变换的性质，从式(24)可推知g(n)的传输特性G(jω)满足：
因此，理论上是G(jω)相移后的版本，两者的幅度传输特性是一致的。既然g(n)可通过式(6)灵活配置，所以的截止频率也可由输入参数ωc被灵活调节。
从式(25)中可推知，在通带内，G(jω)的相位特性与应满足：
因此，它们的群延时与τ(ω)(为与的倒数)满足：
全相位滤波器的群延时τ(ω)等于0，因此可得到：
故的群延时可通过参数p灵活调节。
3)群延时误差分析
基于上述设计分析，本发明实施例的群延时误差主要由两方面因素造成：
1)冲激响应的截断误差；2)三次样条插值的逼近误差。
对于截断误差：在式(20)和(22)中，整数n∈[-N+2,N-2]，这保证了泰勒系数展开的定义域为t＝n+p∈[-N+1,N-1]。然而，当分数时延p给定时，第m个子滤波器的系数a(n,m)与定点t＝n存在一对一的关系。
因此本发明实施例中每个子滤波器的冲激响应都必须经历截断的过程，相较于全相位滤波器的2N-1个抽头系数g(-N+1),…,g(N-1)(通过(6)计算得到)，本发明实施例中的每个子滤波器实际包含2N-3个抽头系数，后面将会证明(实验四)，这个截断误差可忽略不计。
对于逼近误差：因为三次样条插值的阶数为3，所以它很难对一个尖锐的连续函数进行插值操作。连续曲线g(t)会随着截止频率ωc的增大变得尖锐，当ωc高于某个特定频率时，在式(14)中，g(t)与某些分段函数gn(t)之间的误差变得越来越大，这不可避免地降低了群延时的精度。
综上所述，本发明实施例通过上述设计，实现了可变分数时延滤波器Farrow结构抽头系数的快速配置，具有高灵活性(对群延时和截止频率)、高精度与低复杂度。
实施例3
下面结合具体的实例、图2-图6、表2对实施例1、2中的方案进行可行性验证，详见下文描述：
实验一、幅频响应仿真
令输入参数p＝0.3，N＝10，ωc＝0.5π(由式(13)可算得K＝3)，wc(n)为汉宁单窗卷积窗。图4给出了本发明实施例设计的VFD滤波器的幅度传输曲线对比图2和图4可看出，与全相位滤波器的幅频曲线|G(jω)|几乎是重叠的(仅轻微的偏离频率采样点，该小偏离由群延时误差分析中提及的边界截断引起)，这验证了两者幅度响应一致的结论。幅度响应的一致性，保证了群延时的理论特性。
实验二、相频响应仿真
令输入参数p＝0.3，N＝10，ωc＝0.6π，wc(n)为汉宁单窗卷积窗。图5给出了VFD滤波器的群延时曲线与理想群延时曲线τ(ω)＝-pω。可以看出，实际群延时曲线与理想群延时曲线在通带几乎完全重合。从而验证了式(28)的正确性。
实验三、本方法与解析WLS设计的比较
本节对两种解析设计(本方法与WLS设计[14])的群延时误差进行比较。给定相同的时延参数p＝0.3和频率采样整数N＝18，与WLS方法不同(其每个子滤波器耗费2N-1＝35个抽头系数)，本方法仅耗费2N-3＝33个抽头系数；另外，由于采用了3次样条插值(M＝3)，本方法子滤波器个数恒为4，而WLS设计则需考虑不同的M值情况进行测试(文献[14]指出：除去m＝0路外，第(2m-1)路和第2m路的子滤波器系数必须满足特定的解析关系，故M只能为偶数)。其群延时测试误差如表2所示。
表2两方法的群延时误差(p＝0.3,N＝18)
从表2中可总结出如下规律：
1、各种截止频率情况下，本发明实施例的群延时误差均很小(处于10-2级)，相比较而言，本方法更适合于具有低、中截止频率(如ωc＝0.3π，ωc＝0.5π)的VFD滤波器设计。
从表2可看出，随着ωc值升高，本发明实施例的群延时误差增大(从0.0005增大到0.0212)。这是因为，众所周知，低通滤波器的截止频率ωc值越大，其理想的连续冲激响应g(t)波形越陡峭，使得本发明实施例所用的三次样条插值函数对g(t)的逼近误差增大，从而群延时误差也略高于WLS方法的群延时误差。
2、WLS方法仅仅在截止频率为高频值(如ωc＝0.9π)时具有较小的群延时误差，仅适合于截止频率为高频值的VFD滤波器设计。
从表2可看出，当截止频率为低、中频(ωc＝0.3π，ωc＝0.5π)时，WLS方法的群延时误差比本发明实施例高出成百倍，可以说几乎完全不能适用。这是因为，文献[14]指出，WLS方法内含的子滤波器系数约束关系迫使其第1个子滤波器(m＝0)的冲击响应必须是单位响应，即：
a(n,0)＝δ(n) (29)
从而第1路子滤波器G0(z)具有全通的传输特性，也就是说，其截止频率正好等于数字角频率的最高频值π。
进一步，注意式(3)中，子滤波器Gm(z)被pm加权，而分数时延参数p∈[-0.5,0.5]，所以随着m增大，Gm(z)对全局传输函数H(z,p)的影响变得越来越小，远远小于G0(z)对H(z,p)的影响。
综上，G0(z)在全局传输特性H(z,p)中占主导地位，既然G0(z)的截止频率固定为最高值π，迫使整个Farrow结构的截止频率不能偏离π太远，因此，解析的WLS方法仅适用于设计具有高截止频率的VFD滤波器。与此相比，如前所述，本发明实施例的第1路子滤波器G0(z)为全相位滤波器，而不是单位脉冲响应滤波器，故完全不受此限制。
实验四、算法运行时间比较
设定p＝0.3，ωc＝0.7π，WLS设计中的M＝4，将本发明实施例与WLS设计在不同滤波器阶数(N＝8,9,…,31)情况下的算法运行时间进行比较，结果如图6所示。
从图6中可看出，WLS方法耗费的时间随着滤波器阶数变大明显增加(从8.03s升至48.16s)，而本发明实施例的运行时间曲线直观上几乎保持不变，大致都等于4.8s。这是因为，表1所列出的本发明的设计步骤仅仅包含4个解析式的简单计算，完全没有涉及矩阵求逆操作；而WLS方法包含矩阵求逆运算，且求逆矩阵的尺寸随N值增大而变大，故其计算复杂度远远高于本发明。
实施例4
本发明实施例是与实施例1和2中方法实施例对应的装置部分，即，本发明实施例公开了一种低复杂度的可变分数时延滤滤波器，详见下文描述：
在图7中，首先将所需的滤波器分数时延p、截止频率ωc、频率采样向量长度N及归一化的卷积窗wc(n)存入外部RAM中，再将它们实时输入到DSP中，经过DSP内部核心算法，得到全相位滤波器及其传输曲线；利用三次样条插值及泰勒系数展开快速计算得到可变分数时延滤波器Farrow结构的子滤波器系数及其传输曲线，最后由输出驱动及显示电路将其实时显示出来。
其中，图7的DSP(Digital Signal Processor，数字信号处理器)为核心器件，在计算滤波器系数的过程中，完成如下主要功能：
1)调用内部核心算法，对全相位滤波器的解析公式进行构建，计算出所需的滤波器系数，对滤波器系数进行傅里叶变换，得到滤波器传输曲线；
2)控制滤波器参数输入时间，并根据需要实时调整所需要的参数值；
3)将滤波器设计结果实时输出至驱动和显示模块。
DSP器件的内部程序流程如图8所示。
本发明实施例将三次样条插值及泰勒系数展开算法植入DSP器件内，基于此完成高精度、低复杂度、高效的VFD滤波器设计。
图8流程分为如下几个步骤：
(1)首先根据具体需要的分数时延p、滤波器截止频率ωc、频率采样向量长度N及归一化的卷积窗wc(n)计算滤波器设计所需的参数K；
(2)CPU主控器从I/O端口读取滤波器参数，进入内部RAM；
(3)根据推导出的3步设计法进行VFD滤波器设计是DSP算法最核心的部分，运行该算法后，即可得到目标滤波器系数及其滤波器传输曲线；
(4)判断本发明是否满足实际需求，若不满足，程序返回，重新根据要求设定滤波器参数；
(5)直至设计结果符合实际要求，然后通过DSP的输出总线输出至外部显示驱动设备，将滤波器设计结果进行数码显示。
需指出，由于采用了DSP实现，使得整个滤波器设计变得更为灵活快捷，可根据滤波器设计过程中的实际需要，灵活变换滤波器参数，使之最终符合工程需要。
综上所述，本发明实施例通过上述对滤波器的设计，实现了可变分数时延滤波器Farrow结构抽头系数的快速配置，具有高灵活性(对群延时和截止频率)、高精度与低复杂度。
参考文献
[1]陈彩莲，于宏毅，沈彩耀，等.采样率转换中Farrow滤波器实现结构研究[J].信息工程大学学报,2009,10(3):329-332.
[2]Olsson M，Johansson H，Lowenborg P.Time-delay estimation usingFarrow-based fractional-delay FIR filters:Filter approximation vs.estimationerrors[C]//Signal Processing Conference,2006,European.IEEE,2006:1-5.
[3]Pei S C，Tseng C C.A comb filter design using fractional-sampledelay[J].IEEE Trans.on Circuits&Systems II Analog&Digital Signal Processing,1998,45(5):649-653.
[4]Liu G S，Wei C H.A new variable fractional sample delay filter withnonlinear interpolation[J].1992,39(2):123-126.
[5]Deng T B.HIGH-RESOLUTION IMAGE INTERPOLATION USING LAGRANGE-TYPEVARIABLE FRACTIONAL DELAY FILTER[C]//ITC-CSCC:International TechnicalConference on Circuits Systems,Computers and Communic-ations.2006:541-544.
[6]韦文，李宁，汤俊，等.基于分数时延的宽带自适应波束形成[J].清华大学学报(自然科学版),2011,(07):988-992.
[7]Farrow C W.A continuously variable digital delay element[C]//Circuits and Systems,1988.,IEEE International Symposium on.IEEE,1988:2641-2645.
[8]Dam H H.Variable Fractional Delay FIR Filter Design with aBicriteria and Coefficient Relationship[J].IEEE Trans.on Circuits&Systems IIExpress Briefs,2014,61(1):36-40.
[9]Johansson H，Hermanowicz E.Two-Rate Based Low-Complexity VariableFractional Del-ay FIR Filter Structures[J].IEEE Trans.on Circuits&Systems IRegular Papers,2013,60(1):136-149.
[10]Johansson H，Lowenborg P.On the design of adjustable fractionaldelay FIR filters[J].IEEE Trans.on Circuits&Systems II Analog&Digital SignalProcessing,2003,50(4):164-169.
[11]Deng T B.Minimax Design of Low Complexity Allpass VariableFractional Delay Digital Filters[J].Circuits&Systems I Regular Papers IEEETrans.on,2010,57(8):2075-2086.
[12]Hai H D，Cantoni A，Teo K L，et al.Variable Digital Filter WithLeast-Square Criterion and Peak Gain Constraints[J].2007,54(1):24-28.
[13]Deng T B，Lian Y.Weighted-Least-Squares Design of VariableFractional-Delay FIR Filters Using Coefficient Symmetry[J].Signal ProcessingIEEE Trans.on,2006,54(8):3023-3038.
[14]Huang Y D，Pei S C，Shyu J J.WLS Design of Variable Fractional-Delay FIR Filters Using Coefficient Relationship[J].IEEE Trans.on Circuitsand Systems II:Express Briefs,2009,56(3):220-224.
[15]黄翔东,王兆华.基于双相移组合全相位法的FIR陷波器设计[J].系统工程与电子技术,2008,30(1):14-18.
[16]黄翔东,王兆华,吕卫.陷波频率点可精确控制的高效FRM陷波器设计[J].系统工程与电子技术,2009,31(10):2320-2322.
[17]黄翔东,韩溢文,闫子阳,马欣.基于全相位滤波的互素谱分析的高效设计[J].系统工程与电子技术,2017,39(1):23-33.
[18]Huang X，Jing S，Wang Z，et al.Closed-Form FIR Filter Design Basedon Convolution Window Spectrum Interpolation[J].IEEE Trans.on SignalProcessing,2016,64(5):1173-1186.
[19]Huang X，Wang Y，Yan Z，et al.Closed-Form FIR Filter Design withAccurately Controllable Cut-Off Frequency[J].Circuits Systems&SignalProcessing,2016:1-21.
[20]Vaseghi S V.Advanced digital signal processing and noisereduction[M].John Wiley,2000.
本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。
本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种低复杂度的可变分数时延滤波方法，其特征在于，所述方法包括以下步骤：
根据频率采样向量长度、截止频率确定滤波器过渡带和边界整数；
通过归一化的卷积窗、频率采样向量长度、边界整数获取全相位滤波器；
对全相位滤波器进行三次样条插值，获取样条系数；
将分数时延值的不同取值p∈[0,0.5]或p∈[-0.5,0)、以及样条系数带入泰勒解析式中求得可变分数时延滤波器法罗结构的子滤波器系数。
2.根据权利要求1所述的一种低复杂度的可变分数时延滤波方法，其特征在于，所述通过归一化的卷积窗、频率采样向量长度、边界整数获取全相位滤波器具体为：
<mrow> <mi>g</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <msub> <mi>w</mi> <mi>c</mi> </msub> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> <mi>N</mi> </mfrac> <mfrac> <mrow> <mi>sin</mi> <mo>&amp;lsqb;</mo> <mi>n</mi> <mrow> <mo>(</mo> <mn>2</mn> <mi>K</mi> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> <mi>&amp;pi;</mi> <mo>/</mo> <mi>N</mi> <mo>&amp;rsqb;</mo> </mrow> <mrow> <mi>sin</mi> <mrow> <mo>(</mo> <mi>&amp;pi;</mi> <mi>n</mi> <mo>/</mo> <mi>N</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>n</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mo>-</mo> <mi>N</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mo>-</mo> <mn>1</mn> <mo>&amp;rsqb;</mo> <mo>&amp;cup;</mo> <mo>&amp;lsqb;</mo> <mn>1</mn> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mn>1</mn> <mo>&amp;rsqb;</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mn>2</mn> <mi>K</mi> <mo>-</mo> <mn>1</mn> </mrow> <mi>N</mi> </mfrac> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>
其中，g(n)为全相位滤波器；wc(n)为归一化的卷积窗；N为频率采样向量长度；K为边界整数。
3.根据权利要求1或2所述的一种低复杂度的可变分数时延滤波方法，其特征在于，所述边界整数滤波器过渡带Δω＝2π/N；ωc为截止频率。
4.根据权利要求1所述的一种低复杂度的可变分数时延滤波方法，其特征在于，所述对全相位滤波器进行三次样条插值，获取样条系数具体为：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>,</mo> <mn>0</mn> </mrow> </msub> <mo>=</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>,</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mi>g</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <msubsup> <mi>g</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> <mrow> <mo>&amp;prime;</mo> <mo>&amp;prime;</mo> </mrow> </msubsup> <mo>+</mo> <mn>2</mn> <msubsup> <mi>g</mi> <mi>n</mi> <mrow> <mo>&amp;prime;</mo> <mo>&amp;prime;</mo> </mrow> </msubsup> </mrow> <mn>6</mn> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>,</mo> <mn>2</mn> </mrow> </msub> <mo>=</mo> <mfrac> <msubsup> <mi>g</mi> <mi>n</mi> <mrow> <mo>&amp;prime;</mo> <mo>&amp;prime;</mo> </mrow> </msubsup> <mn>2</mn> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>,</mo> <mn>3</mn> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <msubsup> <mi>g</mi> <mrow> <mi>n</mi> <mo>+</mo> <mn>1</mn> </mrow> <mrow> <mo>&amp;prime;</mo> <mo>&amp;prime;</mo> </mrow> </msubsup> <mo>-</mo> <msubsup> <mi>g</mi> <mi>n</mi> <mrow> <mo>&amp;prime;</mo> <mo>&amp;prime;</mo> </mrow> </msubsup> </mrow> <mn>6</mn> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，bn,0、bn,1、bn,2、bn,3分别为样条系数；g(n)为全相位滤波器；g(n+1)为全相位滤波器g(n)以1为单位的时域延拓；g″n与g″n+1分别为在n和n+1两个位置算出数值的二阶导数。
5.根据权利要求1所述的一种低复杂度的可变分数时延滤波方法，其特征在于，当p∈[0,0.5]时，子滤波器系数为：
<mrow> <mi>a</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>,</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <msup> <mi>n</mi> <mi>i</mi> </msup> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>,</mo> <mi>i</mi> </mrow> </msub> <munderover> <mo>&amp;Pi;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>3</mn> <mo>-</mo> <mi>i</mi> </mrow> <mrow> <mi>m</mi> <mo>+</mo> <mn>2</mn> <mo>-</mo> <mi>i</mi> </mrow> </munderover> <mrow> <mo>(</mo> <mn>3</mn> <mo>-</mo> <mi>j</mi> <mo>)</mo> </mrow> <msup> <mi>n</mi> <mrow> <mi>i</mi> <mo>-</mo> <mi>m</mi> </mrow> </msup> </mrow> <mrow> <mi>m</mi> <mo>!</mo> </mrow> </mfrac> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>3</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>n</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mo>-</mo> <mi>N</mi> <mo>+</mo> <mn>2</mn> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mn>2</mn> <mo>&amp;rsqb;</mo> </mrow>
其中，a(n,m)为子滤波器系数；bn,i为样条系数。
6.根据权利要求1所述的一种低复杂度的可变分数时延滤波方法，其特征在于，当p∈[-0.5,0)时，子滤波器系数为：
<mrow> <mi>a</mi> <mrow> <mo>(</mo> <mi>n</mi> <mo>,</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>i</mi> </mrow> </msub> <msup> <mi>n</mi> <mi>i</mi> </msup> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>m</mi> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>0</mn> </mrow> <mn>3</mn> </munderover> <msub> <mi>b</mi> <mrow> <mi>n</mi> <mo>-</mo> <mn>1</mn> <mo>,</mo> <mi>i</mi> </mrow> </msub> <munderover> <mo>&amp;Pi;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>3</mn> <mo>-</mo> <mi>i</mi> </mrow> <mrow> <mi>m</mi> <mo>+</mo> <mn>2</mn> <mo>-</mo> <mi>i</mi> </mrow> </munderover> <mrow> <mo>(</mo> <mn>3</mn> <mo>-</mo> <mi>j</mi> <mo>)</mo> </mrow> <msup> <mi>n</mi> <mrow> <mi>i</mi> <mo>-</mo> <mi>m</mi> </mrow> </msup> </mrow> <mrow> <mi>m</mi> <mo>!</mo> </mrow> </mfrac> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>3</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>,</mo> <mi>n</mi> <mo>&amp;Element;</mo> <mo>&amp;lsqb;</mo> <mo>-</mo> <mi>N</mi> <mo>+</mo> <mn>2</mn> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mn>2</mn> <mo>&amp;rsqb;</mo> </mrow>
其中，a(n,m)为子滤波器系数；bn,i为样条系数。
7.一种用于权利要求1所述的低复杂度的可变分数时延滤波方法的滤波器，其特征在于，所述滤波器包括：外部RAM、数字信号处理器、以及输出驱动及显示电路，
所述外部RAM用于接收滤波器分数时延、截止频率、频率采样向量长度、即归一化的卷积窗；
所述数字信号处理器经过处理得到全相位滤波器及其传输曲线；利用三次样条插值及泰勒系数展开快速计算得到可变分数时延滤波器法罗结构的子滤波器系数及其传输曲线；所述输出驱动及显示电路用于输出显示。
8.根据权利要求7所述的滤波器，其特征在于，所述数字信号处理器为DSP。
CN201710385496.2A 2017-05-26 2017-05-26 一种低复杂度的可变分数时延滤波方法及滤波器 Expired - Fee Related CN107294511B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710385496.2A CN107294511B (zh) 2017-05-26 2017-05-26 一种低复杂度的可变分数时延滤波方法及滤波器

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710385496.2A CN107294511B (zh) 2017-05-26 2017-05-26 一种低复杂度的可变分数时延滤波方法及滤波器

Publications (2)

Publication Number Publication Date
CN107294511A true CN107294511A (zh) 2017-10-24
CN107294511B CN107294511B (zh) 2020-09-08

Family

ID=60094582

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710385496.2A Expired - Fee Related CN107294511B (zh) 2017-05-26 2017-05-26 一种低复杂度的可变分数时延滤波方法及滤波器

Country Status (1)

Country Link
CN (1) CN107294511B (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110504948A (zh) * 2019-08-19 2019-11-26 中国人民解放军陆军工程大学 信号处理方法及等间隔分数时延滤波器
CN111367196A (zh) * 2020-03-05 2020-07-03 上海机电工程研究所 W波段宽带可变分数延时方法及系统
CN112583379A (zh) * 2020-12-04 2021-03-30 北京理工大学 可变小数延时滤波器及信道延时模拟装置
CN113037248A (zh) * 2020-12-04 2021-06-25 北京理工大学 一种基于分段频域优化的小数延时滤波器设计方法
CN114006668A (zh) * 2021-10-29 2022-02-01 中国人民解放军国防科技大学 卫星信道免系数更新的高精确时延滤波方法和装置
CN114189229A (zh) * 2021-11-25 2022-03-15 北京理工大学 一种基于自适应分段算法的小数延时滤波方法
CN114221637A (zh) * 2021-11-25 2022-03-22 北京理工大学 一种基于系数关系结构的分段小数延时滤波方法
CN116187064A (zh) * 2023-02-14 2023-05-30 中国科学院国家空间科学中心 一种连续信号时间序列二阶导数的数值仿真方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN100433779C (zh) * 2006-04-28 2008-11-12 天津大学 全相位离散余弦变换内插核函数及用于图像分辨率的改变
US7339503B1 (en) * 2006-09-29 2008-03-04 Silicon Laboratories Inc. Adaptive asynchronous sample rate conversion
CN102832620A (zh) * 2012-08-31 2012-12-19 天津理工大学 基于加窗全相位fft的apf谐波检测与控制系统及方法
CN103607180B (zh) * 2013-11-27 2016-06-29 东南大学 光电互感器中多数字源采样率转换的快速滤波方法
CN105676171B (zh) * 2016-01-04 2018-06-22 国家无线电监测中心 单通道双基站超短波信号空间定位方法

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110504948B (zh) * 2019-08-19 2022-10-25 中国人民解放军陆军工程大学 信号处理方法及等间隔分数时延滤波器
CN110504948A (zh) * 2019-08-19 2019-11-26 中国人民解放军陆军工程大学 信号处理方法及等间隔分数时延滤波器
CN111367196A (zh) * 2020-03-05 2020-07-03 上海机电工程研究所 W波段宽带可变分数延时方法及系统
CN111367196B (zh) * 2020-03-05 2023-08-18 上海机电工程研究所 W波段宽带可变分数延时方法及系统
CN113037248B (zh) * 2020-12-04 2022-11-11 北京理工大学 一种基于分段频域优化的小数延时滤波器设计方法
CN112583379B (zh) * 2020-12-04 2022-05-20 北京理工大学 可变小数延时滤波器及信道延时模拟装置
CN113037248A (zh) * 2020-12-04 2021-06-25 北京理工大学 一种基于分段频域优化的小数延时滤波器设计方法
CN112583379A (zh) * 2020-12-04 2021-03-30 北京理工大学 可变小数延时滤波器及信道延时模拟装置
CN114006668A (zh) * 2021-10-29 2022-02-01 中国人民解放军国防科技大学 卫星信道免系数更新的高精确时延滤波方法和装置
CN114006668B (zh) * 2021-10-29 2024-02-20 中国人民解放军国防科技大学 卫星信道免系数更新的高精确时延滤波方法和装置
CN114189229A (zh) * 2021-11-25 2022-03-15 北京理工大学 一种基于自适应分段算法的小数延时滤波方法
CN114221637A (zh) * 2021-11-25 2022-03-22 北京理工大学 一种基于系数关系结构的分段小数延时滤波方法
CN116187064A (zh) * 2023-02-14 2023-05-30 中国科学院国家空间科学中心 一种连续信号时间序列二阶导数的数值仿真方法
CN116187064B (zh) * 2023-02-14 2024-03-12 中国科学院国家空间科学中心 一种连续信号时间序列二阶导数的数值仿真方法

Also Published As

Publication number Publication date
CN107294511B (zh) 2020-09-08

Similar Documents

Publication Publication Date Title
CN107294511A (zh) 一种低复杂度的可变分数时延滤波方法及滤波器
CN104283527A (zh) 一种边界频带可快速配置的高效滤波器方法及其装置
CN104579243B (zh) 一种dft调制滤波器组的快速设计方法
Lataire et al. Frequency-domain weighted non-linear least-squares estimation of continuous-time, time-varying systems
Huang et al. Closed-form design of variable fractional-delay FIR filters with low or middle cutoff frequencies
CN105300437B (zh) 一种vlbi基带信号小数时延仿真方法
CN108011615A (zh) 一种信号处理的方法和装置
CN106982045B (zh) 一种基于socp优化的eir-cmfb结构的设计方法
Pei et al. Design of variable comb filter using FIR variable fractional delay element
Goel et al. Analysis of Dirichlet, Generalized Hamming and Triangular window functions in the linear canonical transform domain
Liu et al. Optimal design of frequency-response masking filters with reduced group delays
Koshita et al. High-performance variable band-pass/band-stop state-space digital filters using Gramian-preserving frequency transformation
Dhabu et al. A new time-domain approach for the design of variable FIR filters using the spectral parameter approximation technique
Koshita et al. Variable state-space digital filters using series approximations
Doroslovački On nontrivial analytic signals with positive instantaneous frequency
Petrović et al. New procedure for harmonics estimation based on Hilbert transformation
US7287050B2 (en) Parametric recursive digital filter
Rahmani et al. Improved state-space all-digital filters via series approximations
EP2827439B1 (en) A method for equalizing the distortion caused by losses in couplings in a microwave filter and a filter produced with said method
Milić et al. Design of first–order differentiator utilising FIR and IIR sub–filters
Zhang et al. Design of a novel CEM‐FRM‐based low‐complexity channelised radar transmitter
Litovski et al. Synthesis of Analog Gm-C Hilbert Transformer and Its Implementation for Band-Pass Filter Design
Galván-Sánchez et al. An asynchronous MATE-multirate method for the modeling of electric power systems
CN103178524A (zh) 一种用于有源电力滤波器的数字滤波器
Rasmussen New method of FIR comb filtering

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200908