基于GSPAP的子带自适应声反馈消除算法
梁维谦 1 , 郑方 2 , 陈朝阳 1 , 陈高鋆 3     
1. 集美大学 信息工程学院, 海上通信与智能电子系统福建省高等学校重点实验室, 厦门 361021;
2. 清华大学 信息技术研究院, 语音和语言技术中心, 北京 100084;
3. 厦门莱亚特医疗器械有限公司, 厦门 361009
摘要:声反馈抑制需求鲁棒、高效的自适应滤波器。该文提出一种基于Gauss-Seidel伪放射投影(Gauss-Seidel pseudo affine projection,GSPAP)的子带自适应声反馈消除算法。通过加权重叠相加滤波器组进行子带划分,子带上采用参考信号的自相关矩阵取代能量对滤波器的自适应步长进行归一化;在自相关矩阵的迭代公式中引入基于输入信号和参考信号能量最大值的自整定系数,增强算法鲁棒性;选用二阶GSPAP迭代法对自相关矩阵解算求逆,以平衡算法性能与复杂度。实验结果表明:在相同滤波器长度的条件下,该文方法获得11~22 dB的稳态增益增量,比时域归一化最小均方误差(normalized least mean square,NLMS)方法提升20%~55%。
关键词声反馈抑制    子带    稳态增益    高斯-赛德尔伪放射投影(GSPAP)    
GSPAP based sub-band adaptive feedback cancellation algorithm
LIANG Weiqian1, ZHENG Fang2, CHEN Chaoyang1, CHEN Gaojun3     
1. Key Laboratory of Maritime Communication and Intelligent Electronic Systems of Fujian Province, College of Information Engineering, Jimei University, Xiamen 361021, China;
2. Center for Speech and Language Technologies, Research Institute of Information Technology, Tsinghua University, Beijing 100084, China;
3. Xiamen LA and Associates Medical Equipment Co. LTD., Xiamen 361009, China
Abstract: Acoustic feedback cancellation needs robust and efficient adaptive filters. A Gauss-Seidel pseudo affine projection (GSPAP) sub-band adaptive filter based algorithm was developed using the auto-correlation matrix of the reference signal in the sub-band domain, instead of the power, to normalize the adaptation step size. A self-adjusting coefficient calculated from the maximum power of the input and output sub-band signals, was used to improve the robustness of the auto-correlation matrix. A two-order GSPAP algorithm was used to find the inverse of the auto-correlation matrix to balance the performance and complexity. Results from the added stable gain (ASG) experiment indicate that this algorithm gives 11 to 22 dB improved feedback cancellation for different types of hearing aids, which is 22% to 55% better than given by the time-domain normalized least mean square (NLMS) algorithm.
Key words: acoustic feedback cancellation     sub-band     stable gain     Gauss-Seidel pseudo affine projection (GSPAP)    

目前声反馈抑制的算法主要有2类[1]:1) 基于啸叫检测与抑制的算法,2) 基于自适应滤波的反馈消除算法。

基于啸叫检测与抑制的算法[1-2]是检测系统以及信号的特征,并以此为依据判断系统是否存在啸叫;若检测到啸叫,则对啸叫频率成分进行增益抑制或者对信号相位进行调整,从而破坏啸叫产生的条件,达到消除反馈的目的。用于啸叫检测的声学特征主要有子带能量及其导数、自相关函数的过零率、相位及其导数等[1-2]。由于此类方法要先检测到啸叫发生,而后再进行啸叫抑制处理,因此又被称为被动反馈管理。

基于自适应滤波的反馈消除算法[3-6]则是通过自适应滤波器实时估计声反馈信号,并将其在系统传递函数的分母项中减除。自适应滤波器估计信号越接近于真实声反馈信号,则系统所能提供的稳态增益越高[1]。此类方法根据处理域的不同可分为时域和频域方法。时域方法在时域构建自适应滤波器,处理延时较小,但由于输入声信号和受话器输出声信号之间存在强相关性,造成自适应滤波器的估计偏差问题[1]。因此需要引入白化滤波器去除两者之间的相关性,称作Filter-X方法[3-4]。白化滤波器可以是初始化固定系数或自适应系数的有限冲激响应(finite impulse response, FIR)或无限冲激响应(infinite impulse response, IIR)滤波器[1, 3],或是线性预测编码(linear predictive code, LPC)滤波器[4]。频域方法则是在频域构建自适应滤波器,时频变换相当于对信号进行了白化处理,因此无需额外的去相关步骤[5-6]。频域方法相比于时域方法计算复杂度较高,抑制声反馈的效果也较好,但系统稳定性弱于时域方法。自适应滤波器通常采用归一化最小均方误差(normalized least mean square, NLMS)准则进行迭代。在NLMS的框架下,人们还提出了一系列的约束自适应方法[1, 3, 7],例如将滤波器系数约束在预先测定的初值附近等。

针对NLMS在收敛速度和稳态误差方面的不足[7]以及频域方法的稳定性问题,本文提出一种基于Gauss-Seidel伪仿射投影(Gauss-Seidel pseudo affine projection, GSPAP)[8-9]的子带自适应声反馈消除方法。构建子带自适应滤波器:根据输入信号和输出信号的能量计算出自适应滤波器的自整定系数;由自整定系数稳定迭代更新输出信号的自相关矩阵;在二阶伪仿射投影条件下,推导出自相关矩阵求逆的快速迭代解算公式;由自相关矩阵的逆、输出信号和残差信号迭代更新自适应滤波器系数。本文方法有效提升了系统的最大稳态增益。

1 自适应声反馈消除问题

助听器等开放式电声系统的简化传递函数[1]

$ \begin{array}{*{20}{c}} {Y\left( {f,t} \right) \cong }\\ {S\left( {f,t} \right)\frac{{C\left( {f,t} \right) + G\left( {f,t} \right)\left[ {M\left( f \right)A\left( f \right)R\left( f \right)} \right]}}{{1 - G\left( {f,t} \right)\left[ {M\left( f \right)A\left( f \right)R\left( f \right)B\left( {f,t} \right)} \right]}}.} \end{array} $ (1)

其中:输入声信号S(f, t)经麦克风M(f)转化为模拟电信号;经模数转换得到数字信号,经数字信号处理器G(f, t)得到处理后的数字信号,再经数模转换得到模拟电信号;经功率放大器A(f)、受话器R(f),输出声信号Y(f, t),Y(f, t)通过耳道传输激励耳膜。由于麦克风和受话器之间不可能进行完全的声隔离,因此存在前馈声音通路C(f, t)和反馈声音通路B(f, t),推导得出如式(1) 所示的助听器系统的响应函数[1]。反馈通路B(f, t)的存在,使得系统响应函数出现了可能导致系统不稳定(即引起反馈啸叫)的分母项。为保证式(1) 的分母项不为零,自适应声反馈消除引入自适应滤波器W(f, t),将式(1) 的系统传递函数改写为

$ Y \cong S\frac{{C + G\left[ {W \cdot C + M \cdot A \cdot R} \right]}}{{1 - G\left[ {M \cdot A \cdot R \cdot B - W} \right]}}. $ (2)

其中,为简化公式省略了频率自变量f和时间自变量t

2 子带自适应声反馈消除算法

图 1给出了本文的子带自适应声反馈消除算法框图。输入时域信号s(t)经过加权重叠相加(weighted overlap add, WOLA)[10]分析滤波器组得到时频信号S(n, k);S(n, k)在减去自适应滤波器输出O(n, k)后得到残差信号E(n, k);E(n, k)再经过动态范围压缩、噪声抑制和输出端自动增益控制(automatic gain control-output, AGCo)等其他算法处理后得到输出时频信号Y(n, k);Y(n, k)经过WOLA合成滤波器组[10]得到时域的输出信号y(t),Y(n, k)同时作为自适应滤波器的输入,以最小化E(n, k)均方函数为目标来优化自适应滤波器系数W(n, k)。

图 1 子带自适应声反馈消除算法框图

定义第n帧第k个子带输出信号矢量Yn, k和子带滤波器系数矢量Wn, k分别为

$ {\mathit{\boldsymbol{Y}}_{n,k}} = {\left[ {Y\left( {n,k} \right)Y\left( {n - 1,k} \right) \cdots Y\left( {n - L + 1,k} \right)} \right]^{\rm{T}}}, $ (3)
$ {\mathit{\boldsymbol{W}}_{n,k}} = {\left[ {W\left( {n,k} \right)W\left( {n - 1,k} \right) \cdots W\left( {n - L + 1,k} \right)} \right]^{\rm{T}}}. $ (4)

其中:L表示子带滤波器的长度,(*)T表示转置运算。则残差信号为

$ E\left( {n,k} \right) = S\left( {n,k} \right) - W_{n - 1}^{\rm{H}},{}_k{Y_{n,k}}. $ (5)

其中,(*)H表示共轭转置运算。

WOLA是一种可通过快速Fourier变换实现的具有更高分辨率和更低延时的均匀的频谱分析方法[10]。实验中采用奇型堆叠的WOLA滤波器组,其分析窗长La=4 ms、输入块长度R=1 ms、Fourier变换窗长N=2 ms和合成窗长Ls=2 ms,采样频率fs=16 kHz。

3 基于GSPAP的子带自适应滤波器

本文将S(n, k)作为主信号,Y(n, k)作为参考信号,基于最小二乘法(least squares, LS)求解最优滤波器[1, 3]如下所示:

$ {\mathit{\boldsymbol{w}}_{{\rm{LS}}}} = \mathit{\boldsymbol{\dot R}}_y^{ - 1}{{\mathit{\boldsymbol{\hat c}}}_{sy}}. $ (6)

其中:$\boldsymbol{\hat R}_y^{-1} $表示参考信号的自相关矩阵的逆,${\boldsymbol{\hat c}_{sy}} $表示主信号与参考信号的互相关矢量,wLS表示LS准则下滤波器系数的最优解。下面基于GSPAP[8-9]推导自相关矩阵求逆的迭代公式以及滤波器系数的自适应更新公式,计算步骤如图 2所示。

图 2 基于GSPAP的子带自适应滤波器

1) 根据S(n, k)和Y(n, k)的能量计算出自相关矩阵迭代公式的自整定系数。

a) 计算S(n, k)和Y(n, k)的能量σS2(n, k)和σY2(n, k):

$ \begin{array}{*{20}{c}} {\sigma _S^2\left( {n,k} \right)=}\\ {\left\{ \begin{array}{l} {\left| {S\left( {n,k} \right)} \right|^2},当\;\sigma _S^2\left( {n - 1,k} \right) < {\left| {S\left( {n,k} \right)} \right|^2};\\ \left( {1 - {\alpha _\sigma }} \right){\left| {S\left( {n,k} \right)} \right|^2} + {\alpha _\sigma }\sigma _S^2\left( {n - 1,k} \right),其他. \end{array} \right.} \end{array} $ (7)
$ \begin{array}{*{20}{c}} {\sigma _Y^2\left( {n,k} \right)=}\\ {\left\{ \begin{array}{l} {\left| {Y\left( {n,k} \right)} \right|^2},当\;\sigma _Y^2\left( {n - 1,k} \right) < {\left| {Y\left( {n,k} \right)} \right|^2};\\ \left( {1 - {\alpha _\sigma }} \right){\left| {Y\left( {n,k} \right)} \right|^2} + {\alpha _\sigma }\sigma _Y^2\left( {n - 1,k} \right),其他. \end{array} \right.} \end{array} $ (8)

其中:迭代因子ασ=1-L-1L为子带自适应滤波器的长度。此处采用一阶冲击-释放IIR滤波器迭代计算子带主信号和参考信号的能量。滤波器的冲击时间参数为0,以保证自适应算法的稳定性;释放时间参数ασ的选择与滤波器长度L相关是由于只有最新的L帧信号影响到自相关矩阵的估计和滤波器系数的更新。

b) 计算自整定系数δ(n, k)。通过比较主信号和参考信号能量得到子带的最大能量Δ(n, k):

$ \Delta \left( {n,k} \right) = \max \left\{ {L\left( {M - 1} \right)\sigma _Y^2\left( {n,k} \right),L\sigma _S^2\left( {n,k} \right)} \right\}. $ (9)

其中,M是自适应滤波器伪仿射投影的阶数,M的取值为从1到正无穷大的正整数。理论上,M的取值越大,滤波器跟踪反馈通路的速度越快,精度越高,但相应的计算量也就越大。由Δ(n, k)计算自整定系数δ(n, k):

$ \delta \left( {n,k} \right) = \left\{ \begin{array}{l} \Delta \left( {n,k} \right),当\Delta \left( {n,k} \right) \ge \delta \left( {n - 1,k} \right);\\ \left( {1 - {\alpha _\delta }} \right)\Delta \left( {n,k} \right) + {\alpha _\delta }\delta \left( {n - 1,k} \right),其他. \end{array} \right. $ (10)

其中:αδ=1-(fs/R)-1fs为采样频率,R为输入块长度。此处亦采用一阶冲击-释放IIR滤波器进行迭代。实验表明,释放时间参数αδ的选择无需严格等于1-(fs/R)-1,只需0.5≤αδ < 1,以保证自整定系数不会下降得过快。

选择主信号或参考信号能量的最大值迭代计算自整定系数,用于修正自相关矩阵的迭代公式,如式(10) 所示,避免了可能出现的由于参考信号过小导致自相关矩阵的逆过大从而使得滤波器系数Wn, k更新过快而出现不稳定甚至发散的情况,增强了自适应滤波器的鲁棒性,同时也间接实现了对滤波器步长的控制。

2) 利用上述自整定系数修正Y(n, k)的自相关矩阵${\boldsymbol{\hat R}_{n, k}} $的迭代公式:

$ {{\mathit{\boldsymbol{\hat R}}}_{n,k}} = {{\mathit{\boldsymbol{\hat R}}}_{n - 1,k}} + \mathit{\boldsymbol{\xi }}_{n,k}^ * \mathit{\boldsymbol{\xi }}_{n,k}^T - \mathit{\boldsymbol{\xi }}_{n - L,k}^ * \mathit{\boldsymbol{\xi }}_{n - L,k}^T + \delta \left( {n,k} \right)\mathit{\boldsymbol{I}}. $ (11)

其中:ξn, k=Y(n, k)…Y(nM-1, k)TI是大小为M×M的单位矩阵,M为伪仿射投影阶数,(*)*表示共轭运算。实验表明M取值大于2时对自适应滤波器的收敛效果无明显改进,且会大大增加算法复杂度。因此对算法的性能与复杂度进行折中,取阶数M=2,

$ {{\mathit{\boldsymbol{\hat R}}}_{n,k}} = \left[ {\begin{array}{*{20}{c}} {{{\hat R}_{n,k}}\left( {0,0} \right)}&{{{\hat R}_{n,k}}\left( {0,1} \right)}\\ {{{\hat R}_{n,k}}\left( {1,0} \right)}&{{{\hat R}_{n,k}}\left( {1,1} \right)} \end{array}} \right]. $ (12)

3) 在伪仿射投影阶数M=2的条件下,为进一步简化算法复杂度,定义如下矩阵方程:

$ {{\mathit{\boldsymbol{\hat R}}}_{n,k}}{\mathit{\boldsymbol{P}}_{n,k}} = {\left[ {\begin{array}{*{20}{c}} 1&{{{\bf{0}}^{\rm{T}}}} \end{array}} \right]^{\rm{T}}}. $ (13)

其中:Pn, k=[Pn, k(0)Pn, k(1)… Pn, k(M-1)]T,是大小为M×1的待求解的自相关矩阵的一维逆矩阵;0T是大小为1×(M-1) 的全0矩阵。采用一维Gauss-Seidel迭代法解算逆矩阵Pn, k

$ \left\{ \begin{array}{l} {P_{n,k}}\left( 0 \right) = \frac{{1 - {R_{n,k}}\left( {0,1} \right){P_{n - 1,k}}\left( 1 \right)}}{{{R_{n,k}}\left( {0,0} \right)}},\\ {P_{n,k}}\left( 1 \right) = \frac{{ - {R_{n,k}}\left( {1,0} \right){P_{n,k}}\left( 0 \right)}}{{{R_{n,k}}\left( {1,1} \right)}}. \end{array} \right. $ (14)

4) 利用当前帧的逆矩阵Pn, k解算结果、Y(n, k)和E(n, k)迭代更新Wn, k

$ {\mathit{\boldsymbol{W}}_{n,k}} = {\mathit{\boldsymbol{W}}_{n - 1,k}} + \mu {\mathit{\boldsymbol{V}}_{n,k}}E\left( {n,k} \right). $ (15)

其中:Wn, k表示第n-1帧的滤波器系数矢量,μ为自适应滤波器的控制步长,Vn, k=ζn, kPn, kζn, k的表达式如下:

$ {\mathit{\boldsymbol{\zeta }}_{n,k}} = {\left[ {{Y_{n,k}}{Y_{n - 1,k}} \cdots {Y_{n - M + 1,k}}} \right]^{\rm{T}}}. $ (16)

其中,

$ {\mathit{\boldsymbol{Y}}_{n,k}} = {\left[ {Y\left( {n,k} \right)Y\left( {n - 1,k} \right) \cdots Y\left( {n - L + 1,k} \right)} \right]^{\rm{T}}}. $ (17)

当伪仿射投影阶数M=1时,本文方法退化为NLMS自适应滤波器。

4 实验结果与分析 4.1 方法与步骤

采用稳态增益增量(added stable gain,ASG)指标来检验本文方法的声反馈消除性能[1, 11]。最大稳态增益(maximum stable gain,MSG)定义为助听器即将发生啸叫的临界状态时助听器系统所达到的全频带整体增益。该参数表征了助听器系统所能提供的最大增益,决定了助听器适用的功率验配范围[11]。采用助听器专用测试设备Fonix8000[12],本文按照如下步骤测量ASG:1) 在1 kHz频点校准助听器使得电增益等同于声增益[11];2) 将助听器系统放入Fonix8000的消声箱,受话器接入2CC声耦合腔[11-12],设置输入声信号为60 dB SPL(sound pressure level)的复合噪声信号[12],关闭助听器的反馈消除算法,观察输出声信号的频响曲线;3) 调节助听器的整体电增益,同时观察声增益曲线,找到助听器即将发生啸叫的临界曲线(此时若再增加电增益就会在声增益曲线上出现声压级异常高的啸叫频点),并记录下此时助听器的电增益值作为MSG(单位为dB SPL);4) 开启反馈消除算法,按照步骤3) 的方法重新测量,得到有反馈消除条件下的MSG;5) 将步骤4) 得到的MSG减去步骤3) 得到的MSG,得到反馈消除算法的ASG值(单位为dB)。

声反馈消除的性能与麦克风和受话器之间的声学结构密切相关,因此对ITE(in the ear)、ITC(in the canal)、CIC(completely in the canal)、BTE(behind the ear)、RITE(receiver in the ear)和OE(open ear)等6种助听器声学结构[11]进行了本文方法与时域NLMS方法[1]的对比实验。

4.2 结果与分析

图 3给出了OE类型助听器在无反馈消除时按照节4.1所述步骤测量MSG的示例。曲线1在3.5和7 kHz附近分别出现约80和45 dB SPL 2个尖峰,表示助听器此时已发生啸叫,对助听器的整体电增益进行微调得到临界的曲线2,则曲线2对应的电增益值即为MSG。

图 3 OE类型助听器在无反馈消除时MSG测量示例

表 1给出了本文方法(伪仿射投影阶数M=2,滤波器长度L=1/2/3) 与时域NLMS方法(滤波器长度Lms=2 ms)的ASG实验对比结果。Lms表示以毫秒为单位的时域滤波器长度,它与子带滤波器长度L的函数关系为

表 1 本文方法与时域NLMS的ASG实验结果
实验条件 受话器到麦克风
的延时/ms
无反馈消除
MSG/dB SPL
时域NLMS
ASG/dB
本文方法ASG/dB
L=1 L=2 L=3
BTE 0.50 41 11 11 15 16
ITE 0.31 31 11 10 17 17
ITC 0.31 25 10 8 12 12
RITE 0.41 24 17 15 22 20
CIC 0.25 20 8 5 10 11
OE 0.60 15 12 10 16 18

$ {L_{{\rm{ms}}}} = L \times R \times 1000/{f_{\rm{s}}}. $ (18)

其中:R为WOLA滤波器组的输入块长度,实验中R=1 ms,fs=16 kHz。

影响啸叫产生的因素主要有助听器的增益、喇叭和麦克风之间的隔离度、喇叭和麦克风之间的距离。上述6种助听器按隔离度从大到小依次是BTE、ITE、ITC、CIC、RITE和OE,隔离度越小越易发生啸叫。采用随机相位多重正弦(random phase multisine, RPM)信号[13]测得这6种助听器的受话器到麦克风的延时,延时从大到小依次是OE、BTE、RITE、ITE、ITC和CIC,延时越小说明喇叭和麦克风之间的距离越近也就越容易发生啸叫。综合隔离度和距离因素,测得这6种助听器的MSG从高到低依次是:MSGBTE=41 dB SPL、MSGITE=31 dB SPL、MSGITC=25 dB SPL、MSGRITE=24 dB SPL、MSGCIC= 20 dB SPL和MSGOE=15 dB SPL。对应BTE、ITE、ITC、BITE、CIC和OE,本文方法(L=2) 与时域NLMS方法相比,ASG分别从11、11、10、17、8和12 dB提高到15、17、12、22、10和16 dB,提升20%~55%。本文方法随着滤波器长度L的增加,ASG也有所提升:L=2比L=1提升36%~100%,较为明显;L=3比L=2提升-9%~12.5%,提升相对较小,其中BTE类型还略有下降,然而计算复杂度却明显增加。

表 2给出了本文方法(L=2) 的伪仿射投影阶数的ASG实验结果。随着伪仿射投影阶数M的增加,ASG也有所提升:M=2比M=1提升9%~89%,较为明显;M=3比M=2提升-10%~6%,提升相对较小,其中ITC类型还略有下降,然而计算复杂度却显著增加。因此本文最终选择滤波器长度L=2,伪仿射投影阶数M=2。

表 2 本文方法伪仿射投影阶数M的ASG实验结果
实验条件 无反馈消除
MSG/dB SPL
ASG/dB
M=1 M=2 M=3
BTE 41 9 17 18
ITE 31 11 12 12
ITC 25 8 10 9
RITE 24 12 15 18
CIC 20 15 22 23
OE 15 10 16 16

5 结论

本文提出一种基于GSPAP的子带自适应声反馈消除方法。在子带上采用自相关矩阵取代NLMS方法中的参考信号能量对滤波器的自适应迭代项进行归一化,收敛速度更快且稳态误差也更小;在自相关矩阵的迭代公式中引入基于输入信号和输出信号能量最大值的自整定系数,避免了可能出现的由于参考信号过小导致自相关矩阵的逆过大从而使得滤波器系数更新过快而出现不稳定甚至发散的情况,增强了自适应滤波器的鲁棒性,也起到了变步长的作用;采用二阶Gauss-Seidel迭代法对自相关矩阵解算求逆,进一步降低算法复杂度。实验结果表明:在相同滤波器长度的条件下,本文方法比时域NLMS方法对ASG有明显提升(20%~55%)。综合考虑算法性能和复杂度,最终选择滤波器长度L=2,伪仿射投影阶数M=2。

参考文献
[1] Kates J M. Digital Hearing Aids[M]. San Diego, USA: Plural Publishing INC, 2008.
[2] Nakagawa C R C, Nordholm S, YAN Weiyang. Feedback cancellation with probe shaping compensation[J]. IEEE Signal Processing Letters, 2014, 21(3): 365–369. DOI:10.1109/LSP.2014.2301832
[3] Hellgren J. Analysis of feedback cancellation in hearing aids with Filter-X LMS and the direct method of closed loop identification[J]. IEEE Transactions on Speech and Audio Processing, 2002, 10(2): 119–131. DOI:10.1109/89.985549
[4] MA Guilin. Feedback Suppression in Digital Hearing Instruments [D]. Ballerup, Denmark: Technical University of Denmark, 2010. https://www.researchgate.net/publication/15114361_Digital_Feedback_Suppression_DFS_Clinical_Experiences_when_Fitting_a_DFS_Hearing_Instrument_on_Children
[5] Wyrsch S, Kaelin A. Adaptive feedback cancellation in subbands for hearing aids [C]//Proceedings of the International Conference on Acoustic, Speech and Signal Processing (ICASSP'99). Phoenix, AZ, USA: IEEE, 1999: 921-924.
[6] Gil-Cacho J M, Waterschoot V T, Moonen M, et al. Transform domain prediction error method for improved acoustic echo and feedback cancellation [C]//Proceedings of the 20th European Signal Processing Conference. Bucharest, Romania: EURASIP, 2012: 2422-2426.
[7] Vanderkooy J. Aspects of MLS measuring systems[J]. Journal of the Audio Engineering Society, 1994, 42(4): 219–231.
[8] Lee J, Park Y C, Youn D H. Robust pseudo affine projection algorithm with variable step-size[J]. Electronic Letters, 2008, 44(3): 250–252. DOI:10.1049/el:20083379
[9] Sheikhzadeh H, Whyte K R L, Brennan R L. Partial update subband implementation of complex pseudo-affine projection algorithm on oversampled filterbanks [C]//Proceedings of the International Conference on Acoustic, Speech and Signal Processing (ICASSP 2005). Philadelphia, PA, USA: IEEE, 2005, 4: 373-376.
[10] On Semiconductor. WOLA Filterbank Coprocessor: Introductory Concepts and Techniques[R/OL]. April, 2009. http://www.onsemi.com/pub/Collateral/AND8382-D.PDF.
[11] Valente M. Hearing Aids: Standards, Options, and Limitations[M]. 2nd ed, New York, USA: Thieme Medical Publishers, 2002.
[12] Frye Electronics INC. FONIX 8000 Hearing Aid Test System Operator's Manual Version 2.41[R/OL]. August, 2013. http://www.frye.com/wp/wp-content/uploads/2013b/manuals/8000v2.41.pdf.
[13] 厦门莱亚特医疗器械有限公司. 一种测量助听器反馈路径的方法: 中国, CN201410685349[P]. 2015. 01. 28. Xiamen LA and Associates Medical Equipment Co Ltd. A Method for Measuring the Feedback Path of Hearing Aids: China, CN201410685349[P]. 2015.01.28. (in Chinese)