参考链接：
语音自适应滤波器的设计
说说鸡尾酒会问题和程序实现
【PUDN】自适应变步长的龙格库塔法
参考文献：
[1]黄颖. 基于麦克风阵列手机消噪方案的应用与实现[D].上海交通大学,2012.
[2]王超. 频域卷积混合盲分离研究[D].上海大学,2008.
【论文】自适应谱线增强器在生物雷达中的应用

一、理论基础

在实际应用中，常常无法的带信号和噪声统计特性的鲜艳信息，在此情况下自适应滤波技术能够获得极佳的滤波性能，因此具有很好的应用价值，常用的自适应滤波器技术有：最小均方（LMS）自适应滤波器，递推最小二乘（RLS）滤波器，格型滤波器和无限冲激响应（IIR)滤波器，。这些自适应滤波技术的应用又包括：自适应噪声抵消，自适应谱线增强和陷波。
LMS自适应滤波器事使滤波器的输出信号于期望响应之间的误差的均方值为最小，因此称为最小均方（LMS）自适应滤波器。

1、基于LMS算法

构成自适应数字滤波器的基本部件是自适应线性组合器，设线性组合的M个输入为x(k-1)…,x(k-M),其输出有y(k)是这些输入加权后的线性组合，即
∑i=0MWix(k−i)\sum_{i=0}^{M}{W_{i}x(k-i)}i=0∑MWix(k−i)
定义权向量，W=[W1,W2,W3,...,Wm]W=[W_{1},W_{2},W_{3},...,W_{m}]W=[W1,W2,W3,...,Wm]X(k)=[x(k−1)...,x(k−M)]TX(k)=[x(k-1)...,x(k-M)]_{}^{T}X(k)=[x(k−1)...,x(k−M)]T令d(k)代表“所期望的响应”，并定义误差信号ε(k)=d(k)−y(k)=d(k)−∑i=0MWix(k−i)\varepsilon(k)=d(k)-y(k)=d(k)-\sum_{i=0}^{M}{W_{i}x(k-i)}ε(k)=d(k)−y(k)=d(k)−i=0∑MWix(k−i)写成向量形式为：ε(k)=d(k)−WTX(k)=d(k)−XTW(k)\varepsilon(k)=d(k)-W_{}^{T}X(k)=d(k)-X_{}^{T}W(k)ε(k)=d(k)−WTX(k)=d(k)−XTW(k)误差平方为：ε2(k)=d2(k)−2d(k)XT(k)W+WTX(k)XT(k)W\varepsilon^2(k)=d^2(k)-2d(k)X^T(k)W+W^TX(k)X^T(k)Wε2(k)=d2(k)−2d(k)XT(k)W+WTX(k)XT(k)W上式两边取数学期望为，得均方误差为
E(ε2(k))=E(d2(k))−2E(d(k)XT(k)W)+WTE(X(k)XT(k))WE(\varepsilon^2(k))=E(d^2(k))-2E(d(k)X^T(k)W)+W^TE(X(k)X^T(k))WE(ε2(k))=E(d2(k))−2E(d(k)XT(k)W)+WTE(X(k)XT(k))W定义互相关函数行向量RxdTR_{xd}^{T}RxdT:
RxdT=E(d(k)XT(k))——（1−9）R_{xd}^{T}=E(d(k)X^T(k))——（1-9）RxdT=E(d(k)XT(k))——（1−9）和自相关函数矩阵
RXX=E(X(k)XT(k))——（1−10）R_{XX}=E(X(k)X^T(k))——（1-10）RXX=E(X(k)XT(k))——（1−10）则均方误差可更改表述为
E(ε2(k))=E(d2(k))−2RxdTW+WTRXXW——（1−11）E(\varepsilon^2(k))=E(d^2(k))-2R_{xd}^{T}W+W^TR_{XX}W——（1-11）E(ε2(k))=E(d2(k))−2RxdTW+WTRXXW——（1−11）
利用式1-10求最佳权系数向量得精确解需要知道RXXR_{XX}RXX和RxdR_{xd}Rxd的先验统计知识，而且还需要进行矩阵求逆等运算。这种算法的根据是最优化方法中的最速下降法。根据最速下降法，“下一时刻”权系数向量W(k+1)应该等于“现时刻”权系数向量W（k）加上一个负均方误差梯度−▽(k)-\triangledown(k)−▽(k)的比例项，即W(k+1)=W(k)−μ▽(k)W(k+1)=W(k)-\mu\triangledown(k)W(k+1)=W(k)−μ▽(k)式中，μ\muμ是一个控制收敛速度于稳定性的常数，称之为收敛因子。
不难看出，LMS算法有两个关键，梯度▽(k)\triangledown(k)▽(k)的计算以及收敛因子μ\muμ的选择。

1.1▽(k)\triangledown(k)▽(k)的近似计算

精确计算梯度▽(k)\triangledown(k)▽(k)是十分苦难的，一种粗略的但是却十分有效的计算▽(k)\triangledown(k)▽(k)的近似方法是:直接取ε2(k)\varepsilon^2(k)ε2(k)作为均方误差E(ε2(k))E(\varepsilon^2(k))E(ε2(k))的估计值，即
▽^(k)=▽[ε2(k)]=2ε(k)▽(k)\hat{\triangledown}(k)=\triangledown[\varepsilon^2(k)]=2\varepsilon(k)\triangledown(k)▽^(k)=▽[ε2(k)]=2ε(k)▽(k)得到梯度估计值
▽^(k)=−2ε(k)X(k)\hat{\triangledown}(k)=-2\varepsilon(k)X(k)▽^(k)=−2ε(k)X(k)于是，Widrow-Hoff LMS算法最终为
W(k+1)=W(k)+2με(k)X(k)——（1−14）W(k+1)=W(k)+2\mu\varepsilon(k)X(k)——（1-14）W(k+1)=W(k)+2με(k)X(k)——（1−14）式子1-14的实现方框图如下图所示

下面分析梯度估计值▽^(k)\hat{\triangledown}(k)▽^(k)的无偏性。▽^(k)\hat{\triangledown}(k)▽^(k)的数学期望为
——（1-15）
在上面的推导过程中，利用了d(k)和ε(k)\varepsilon(k)ε(k)二者皆为标量的事实。在得到最后的结果时，利用了式子1-9和1-15表明，梯度估值▽^(k)\hat{\triangledown}(k)▽^(k)是无偏估计。

1.2μ\muμ的选择

对权系数向量更新公式1-14两边取数学期望，得
E(W(k+1))=E(W(lk))+2μE(ε(k)X(k))=(I−2μRXX)E(W(k))+2μRxdE(W(k+1))=E(W(lk))+2\mu E(\varepsilon(k)X(k))=(I-2\mu R_{XX})E(W(k))+2\mu R_{xd} E(W(k+1))=E(W(lk))+2μE(ε(k)X(k))=(I−2μRXX)E(W(k))+2μRxd
式子中，I为单位矩阵，Rxd=E(d(k)X(k))R_{xd}=E(d(k)X(k))Rxd=E(d(k)X(k))和RXX=E(X(k)XT(k))R_{XX}=E(X(k)X^T(k))RXX=E(X(k)XT(k))。
当时，k=0k=0k=0时
E(W(1))=(I−2μRXX)E(W(0))+2μRxdE(W(1))=(I-2\mu R_{XX})E(W(0))+2\mu R_{xd}E(W(1))=(I−2μRXX)E(W(0))+2μRxd
当时，k=1k=1k=1时,利用上式结果，则有
E(W(2))=(I−2μRXX)E(W(1))+2μRxd(I−2μRXX)E(W(0))+2μ∑i=01(I−2μRXX)iRxdE(W(2))=(I-2\mu R_{XX})E(W(1))+2\mu R_{xd}(I-2\mu R_{XX})E(W(0))+2\mu\sum_{i=0}^{1}(I-2\mu R_{XX})^iR_{xd}E(W(2))=(I−2μRXX)E(W(1))+2μRxd(I−2μRXX)E(W(0))+2μi=0∑1(I−2μRXX)iRxd
起始时，E(W(0))=W(0)E(W(0))=W(0)E(W(0))=W(0)
故重复以上迭代至k+1，则有
E(W(k+1))=(I−2μRXX)k+1W(0)+2μ∑i=01(I−2μRXX)iRxd——（1−18）E(W(k+1))=(I-2\mu R_{XX})^{k+1}W(0)+2\mu\sum_{i=0}^{1}(I-2\mu R_{XX})^iR_{xd}——（1-18）E(W(k+1))=(I−2μRXX)k+1W(0)+2μi=0∑1(I−2μRXX)iRxd——（1−18）由于RXXR_{XX}RXX是实值得对称阵，我们可以写出其特征值分解式
RXX=Q∑QT=Q∑Q−1——（1−19）R_{XX}=Q\sum{Q^T}=Q\sum{Q^{-1}}——（1-19）RXX=Q∑QT=Q∑Q−1——（1−19）这里，我们利用了正定阵Q的性质QT=Q−1Q^T=Q^{-1}QT=Q−1,且∑=diag(λ1,...λM)\sum=diag(\lambda_{1},...\lambda_{M})∑=diag(λ1,...λM)是对角阵，其对角元素λi\lambda_{i}λi是RXXR_{XX}RXX的特征值。将式1-18带入式1-19后得
E(W(k+1))=(I−2μQ∑Q−1)k+1W+2μ∑i=0k(I−2μQ∑Q−1)iRxd——（1−20）E(W(k+1))=(I-2\mu Q\sum{Q^{-1}})^{k+1}W+2\mu\sum_{i=0}^{k}(I-2\mu Q\sum{Q^{-1}})^iR_{xd}——（1-20）E(W(k+1))=(I−2μQ∑Q−1)k+1W+2μi=0∑k(I−2μQ∑Q−1)iRxd——（1−20）
E(W(k+1))=Q∑−1Q−1Rxd=RXX−1Rxd=Wopt——（1−20）E(W(k+1))=Q\sum^{-1}Q^{-1}R_{xd}=R_{XX}^{-1}R_{xd}=W_{opt}——（1-20）E(W(k+1))=Q∑−1Q−1Rxd=RXX−1Rxd=Wopt——（1−20）由此可见，当迭代次数无限增加时，权系数向量得数学期望可收敛至wiener解，其条件是对角阵(I−2μ∑)(I-2\mu\sum)(I−2μ∑)的所有对角元素均小于1，即0<μ<1λmax0<\mu<\frac{1}{\lambda_{max}}0<μ<λmax1
其中，λmax\lambda_{max}λmax是RXXR_{XX}RXX的最大特征值。μ\muμ称为收敛因子，它决定达到式1-20的速率。事实上，W(k)W(k)W(k)收敛于WoptW_{opt}Wopt由比值d=λmax/λmind=\lambda_{max}/\lambda_{min}d=λmax/λmin决定，该比值叫做谱动态范围。大的ddd值预示要花很长的时间才会收敛到最佳权值。克服这一困难的方法之一是产生正交数据。

验证函数直接使用LMS实现成matlab中的一个函数，在测试LMS算法中调用LMS函数即可实现以下图形。

LMS算法matlab实现：

function [yn,W,en] = LMS(xn,dn,M,mu,itr)
%LMS(least mean squre)最小均方差算法
%%%%%%%%%%%%%输入参数：%%%%%%%%%%%%%%%%%%%
%xn 输入信号的序列（列向量）
%dn 所期望的响应序列（列向量）
%M  滤波器的阶数（标量）
%mu 收敛因子（步长）（标量）要求大于0，小于xn的相关矩阵最大特征值的倒数1/lamuda[max]
%itr 迭代次数（标量）默认为xnde changdu ,M<itr<length(xn)
%%%%%%%%%%%%%%输出参数：%%%%%%%%%%%%%%%%%
%W 滤波器的权值矩阵（矩阵）大小为M*itr
%en 误差序列（itr*1）(列向量)
%yn 实际输出序列（列向量）
%%%%%%%%%%%%%%%参数个数必须时4个或者5个%%%%%%%%%%%%%
if nargin == 4itr = length(xn);
elseif nargin = 5if itr>length(xn) | itr<Merror('迭代次数过大或过小');end
else error('请检查输入参数的个数');
end
%%%%%%%%%%%%初始化参数%%%%%%%%%%%%
en = zero(itr,1);
W = zeros(M,itr);
%%%%%%%%%%%%迭代计算%%%%%%%%%%%%
for k = M:itr            %第k次迭代x = xn(k:-1:k-M-+1); %滤波器M个抽头的输入y = W(:,k-1).'* x;   %滤波器的输出en(k) = dn(k) - y;   %第k次迭代的误差%%%%%滤波器权值计算的迭代值%%%%%W(:,k) = W(:,k-1) + 2*mu*en(k)*x;
end
%%%%%%%%%%%%求最优时滤波器的输出序列%%%%%%%%%%
yn = inf *ones(size(xn));
for k = M:length(xn)x = xn(k:-1:k-M+1);yn(k) = W(:,end).'* x;
end

LMS算法测试

%%%%%%%%%%%%%验证所设计的自适应滤波器%%%%%%%%%%%
%%%%%function main()
close all;
%%%调用周期信号%%%%%
t = 0:99;
xs = 10*sin(0.5*t);
figure(1);
subplot(2,1,1);
plot(t,xs);grid;
ylabel('幅值');
title('it{输入周期性信号}');
%%%%%噪声产生%%%%%
randn('state',sum(100*clock));
xn = randn(1,100);
subplot(2,1,2);
plot(t,xn);grid;
ylabel('幅值');
xlabel('时间');
title('it{随机噪声信号}');
%%%%%信号滤波%%%%%
xn = xs + xn;
xn = xn.';       %输入信号序列
dn = xs.';
M = 20';
lambda_max = max(eig(xn*xn.'));
mu = rand()*(1/lambda_max);
[yn,W,en] = LMS(xn,dn,M,mu);
%%%%绘制滤波器输入信号%%%%%%
figure(2);
subplot(2,1,1);
plot(t,xn);grid;
ylabel('幅值');
xlabel('时间');
title('it{滤波器输入信号}');
%%%%绘制滤波器输出信号%%%%%%
subplot(2,1,2);
plot(t,yn);grid;
ylabel('幅值');
xlabel('时间');
title('it{自适应滤波器输出信号}');
%%%绘制自适应滤波器输出信号，预期输出信号和两者的误差%%%%%
figure(3);
plot(t,yn,'b',t,dn,'g',t,dn-yn,'r');grid;
legend('自适应滤波器输出','预期输出误差');
ylabel('幅值');
xlabel('时间');
title('it{自适应滤波器}');

RLS算法matlab实现：

[primary,Fs,Nbits] = wavread('RLSprimsp.wav');
x = primary';               %含噪声的信号
ref = wavread('RLSrefns.wav'); r = ref';                   % 输入噪声序列
M = 32 ;                    %设置序列长度
L = length (x);             %数据的长度
Lambda = 0.999;             %遗忘因子
Delta = 0.001 ;             %初始化因子
ab = rand(M,1);             % 产生均匀分布的随机数矩阵%RLS算法
y = filter (ab,1,x);
E = eye (M,M);              %生成单位矩阵
IR = Delta * E ;
wn = zeros (M,1) ;          %把权值初始化
for n = M : L    xn = x(n);                        %当前x(n)rn = r(n:-1:n-M+1) ;                %当前r(n)en = xn - rn*wn;g = rn*IR;IR = 1/Lambda*[IR - (IR*rn'*rn*IR)/(Lambda+rn*IR*rn')];wn = wn + g' .* en;                 v = rn * wn;e(n) = xn - v;
end
for n = M : L    xn = x(n);                      %当前x(n)rn = r(n:-1:n-M+1) ;            %当前r(n)v = rn * wn;e(n) = xn - v;
end
figure(1)           %建立图形窗口
subplot(2,1,1);
plot(x);        %建立线性图形
subplot(2,1,2);
plot(r);figure(2)
subplot(3,1,1);
plot(e);
subplot(3,1,2);
spy1=2*abs(fft(x))/L;
spy1=spy1(1:L/2);
df=Fs/L;
f=(0:L/2-1)*df;
plot (f,spy1);
subplot(3,1,3);
N = length(e);
spy2=2*abs(fft(e))/N;
spy2=spy2(1:N/2);
df=Fs/N;
f=(0:N/2-1)*df;
figure(3);
plot (f,spy2);
xlabel ('f(Hz)');
ylabel('原始语音的幅度谱');
audiowrite('eaudio.wav',e,Fs);

【算法】麦克风阵列的自适应降噪算法相关推荐

基于麦克风阵列的声源定位算法之GCC-PHAT
目前基于麦克风阵列的声源定位方法大致可以分为三类: 基于最大输出功率的可控波束形成技术基于高分辨率谱图估计技术和基于声音时间差(time-delay estimation,TDE). 基于TDE的算 ...
麦克风声源定位原理_基于麦克风阵列的声源定位算法之GCC-PHAT
目前基于麦克风阵列的声源定位方法大致可以分为三类:基于最大输出功率的可控波束形成技术.基于高分辨率谱图估计技术和基于声音时间差(time-delay estimation,TDE)的声源定位技术. 基 ...
android降噪算法,小米科普自研降噪算法：可精准分离人声与环境噪音并实现降噪...
集微网8月25日消息(文/数码控),小米手机官方微博在昨天发表长文科普小米10 至尊纪念版搭载的自研降噪算法,官方称它可精准分离人声与环境噪音并实现降噪,再喧闹的车水马龙,与噪音"正面交手& ...
OpenCV —— 阈值分割（直方图技术法，熵算法，Otsu，自适应阈值算法）
阈值分割 1. 全局阈值分割直方图技术法熵算法 Otsu算法 2. 局部阈值分割自适应阈值阈值的分割的核心就是如何选取阈值,选取正确的阈值时分割成功的关键.可以使用手动设置阈值,也可以采用直方 ...
阵列算法matlab,阵列信号处理的常见算法
阵列处理算法常用程序 AR argamse.m armaorder.m lsar.m lsarma.m myprogramme.m mywarma.m yulewalker.m ESPRIT TAM算 ...
麦克风阵列语音增强beamforming算法
delay and sum 关键步骤在于计算延时, 可以通过GCC-PHAT方法进行计算, 即广义互相关-相位变换方法. GCC-PHAT(广义互相关-相位变换) x(n) x(n) 和y(n) y( ...
麦克风阵列语音增强算法——固定波束形成算法
与单一麦克风不同,麦克风阵列除了能区分接收到的语音信号的时域和频域特性之外,还能区分空间特效,能在嘈杂的语音环境中在特定的方向上形成波束来获得特定声源发出的语音信号,并且能有效抑制噪声.波束形成方法可 ...
图像降噪算法——DnCNN / FFDNet / CBDNet / RIDNet / PMRID / SID
图像降噪算法--DnCNN / FFDNet / CBDNet / RIDNet / PMRID / SID 图像降噪算法--DnCNN / FFDNet / CBDNet / RIDNet / PM ...
图像降噪算法——维纳滤波
图像降噪算法--维纳滤波图像降噪算法--维纳滤波 1. 基本原理 2. C++代码实现 3. 结论图像降噪算法--维纳滤波维纳滤波是在频域中处理图像的一种算法,是一种非常经典的图像增强算法,不仅 ...
图像降噪算法——低秩聚类：WNNM算法
图像降噪算法--低秩聚类:WNNM算法图像降噪算法--低秩聚类:WNNM算法 1. 基本原理 2. matlab代码 3. 结论图像降噪算法--低秩聚类:WNNM算法同样是为了完善自己知识版图的 ...

【算法】麦克风阵列的自适应降噪算法

一、理论基础

1、基于LMS算法

1.1▽(k)\triangledown(k)▽(k)的近似计算

1.2μ\muμ的选择

【算法】麦克风阵列的自适应降噪算法相关推荐

最新文章

热门文章