MATLAB自适应滤波去噪

《MATLAB自适应滤波去噪》
课程设计报告
1.课程设计目的
此次课程设计目的是为了让我们学会使用MATLAB进行计算机仿真，使用自适应滤波法设计一个语音去噪声电路。

培养我们的电路设计思路及其算法，明白理论与实践相结合的重要性，培养了我们的实际操作能力以及锻炼我们对实际问题的分析与解决的能力。

2.课程设计内容
2.1 LMS自适应算法原理
自适应过程一般采用典型LMS自适应算法，但当滤波器的输入信号为有色随机过程时，特别是当输入信号为高度相关时，这种算法收敛速度要下降许多，这主要是因为输入信号的自相关矩阵特征值的分散程度加剧将导致算法收敛性能的恶化和稳态误差的增大。

此时若采用变换域算法可以增加算法收敛速度。

变换域算法的基本思想是：先对输入信号进行一次正交变换以去除或衰减其相关性，然后将变换后的信号加到自适应滤波器以实现滤波处理，从而改善相关矩阵的条件数。

因为离散傅立叶变换DFT本身具有近似正交性，加之有FFT快速算法，故频域分块LMS FBLMS算法被广泛应用。

FBLMS算法本质上是以频域来实现时域分块LMS算法的，即将时域数据分组构成N个点的数据块，且在每块上滤波权系数保持不变。

其原理框图如图2所示。

FBLMS 算法在频域内可以用数字信号处理中的重叠保留法来实现，其计算量比时域法大为减少，也可以用重叠相加法来计算，但这种算法比重叠保留法需要较大的计算量。

块数据的任何重叠比例都是可行的，但以50％的重叠计算效率为最高。

对FBLMS算法和典
型LMS算法的运算量做了比较，并从理论上讨论了两个算法中乘法部分的运算量。

本文从实际工程出发，详细分析了两个算法中乘法和加法的总运算量，其结果为：
复杂度之比＝FBLMS实数乘加次数／LMS实数乘加次数＝（25Nlog2N＋2N－4）／[2N(2N－1)]
采用ADSP的C语言来实现FBLMS算法的程序如下：
for(i＝0;i＜＝30;i＋＋)
{for(j＝0;j＜＝n－1;j＋＋)
{in[j]＝input[i×N＋j;]
rfft(in,tin,nf,wfft,wst,n);
rfft(w,tw,wf,wfft,wst,n);
cvecvmlt(inf,wf,inw,n);
ifft(inw,t,O,wfft,wst,n);
for(j＝0,j＜＝N－1;j＋＋)
{y[i×N＋j]＝O[N＋j]．re;
e[i×N＋j]＝re fere[i×N＋j]－y[i×N＋j];
temp[N＋j]＝e[i×N＋j;}
rfft(temp,t,E,wfft,wst,n);
for(j＝0;j＜＝n－1;j＋＋)
{inf_conj[j]＝conjf(inf[j]);}
cvecvmlt(E,inf_conj,Ein,n);
ifft(Ein,t,Ein,wfft,wst,n);
for(j＝0;j＜＝N－1;j＋＋)
{OO[j]＝Ein[j]．re;
w[j]＝w[j]＋2＊u＊OO[j];}
}
在EZ－KIT测试板中，笔者用汇编语言和C语言程序分别测试了典型LMS算法的运行速度，并与FBLMS算法的C语言运行速度进行了比较，表2所列是其比较结果，从表2可以看出滤波器阶数为64时，即使是用C语言编写的FBLMS算法也比用汇编编写的LMS算法速度快20％以上，如果滤波器的阶数更大，则速度会提高更多。

2.2 语音信号去噪声源程序
%lms算法源程序
clear all
close all
%channel system order
sysorder = 5 ;
% Number of system points
N=2000;
inp = randn(N,1);
n = randn(N,1);
[b,a] = butter(2,0.25);
Gz = tf(b,a,-1);
%This function is submitted to make inverse Z-transform (Matlab central file exchange)
%The first sysorder weight value
%h=ldiv(b,a,sysorder)';
% if you use ldiv this will give h :filter weights to be
h= [0.0976;
0.2873;
0.3360;
0.2210;
0.0964;];
y = lsim(Gz,inp);
%add some noise
n = n * std(y)/(10*std(n));
d = y + n;
totallength=size(d,1);
%Take 60 points for training
N=60 ;
%begin of algorithm
w = zeros ( sysorder , 1 ) ;
for n = sysorder : N
u = inp(n:-1:n-sysorder+1) ;
y(n)= w' * u;
e(n) = d(n) - y(n) ;
% Start with big mu for speeding the convergence then slow down to reach the correct weights
if n < 20
mu=0.32;
else
mu=0.15;
end
w = w + mu * u * e(n) ;
end
%check of results
for n = N+1 : totallength
u = inp(n:-1:n-sysorder+1) ;
y(n) = w' * u ;
e(n) = d(n) - y(n) ;
end
hold on
plot(d)
plot(y,'r');
title('System output') ;
xlabel('Samples')
ylabel('True and estimated output')
figure
semilogy((abs(e))) ;
title('Error curve') ;
xlabel('Samples')
ylabel('Error value')
figure
plot(h, 'k+')
hold on
plot(w, 'r*')
legend('Actual weights','Estimated weights')
title('Comparison of the actual weights and the estimated weights') ; axis([0 6 0.05 0.35])
% RLS 算法
randn('seed', 0) ;
rand('seed', 0) ;
NoOfData = 8000 ; % Set no of data points used for training Order = 32 ; % Set the adaptive filter order
Lambda = 0.98 ; % Set the forgetting factor
Delta = 0.001 ; % R initialized to Delta*I
x = randn(NoOfData, 1) ;% Input assumed to be white
h = rand(Order, 1) ; % System picked randomly
d = filter(h, 1, x) ; % Generat
e output (desired signal)
% Initialize RLS
P = Delta * eye ( Order, Order ) ;
w = zeros ( Order, 1 ) ;
% RLS Adaptation
for n = Order : NoOfData ;
u = x(n:-1:n-Order+1) ;
pi_ = u' * P ;
k = Lambda + pi_ * u ;
K = pi_'/k;
e(n) = d(n) - w' * u ;
w = w + K * e(n) ;
PPrime = K * pi_ ;
P = ( P - PPrime ) / Lambda ;
w_err(n) = norm(h - w) ;
end ;
% Plot results
figure ;
plot(20*log10(abs(e))) ;
title('Learning Curve') ;
xlabel('Iteration Number') ;
ylabel('Output Estimation Error in dB') ;
figure ;
semilogy(w_err) ;
title('Weight Estimation Error') ;
xlabel('Iteration Number') ;
ylabel('Weight Error in dB') ;
2.3 去噪声前后信号波形
System output
Error curve
Comparison Learning Curve
Weight Estimation Error
附录资料：MATLAB的30个方法1 内部常数
2 数学运算符
3 关系运算符
4 常用内部数学函数
acsch(x) 反双曲余割函数
求角度函数atan2(y,x) 以坐标原点为顶点，x轴正半轴为始边，从原点到点（x，y）的射线为终边的角，其单位为弧度，范围为（，]
数论函数gcd(a,b) 两个整数的最大公约数lcm(a,b) 两个整数的最小公倍数
排列组合函数factorial(n)
阶乘函数，表示n的阶乘
复数函数real(z) 实部函数
imag(z) 虚部函数
abs(z) 求复数z的模
angle(z)
求复数z的辐角，其范围是
（，]
conj(z) 求复数z的共轭复数
求整函数与截尾函数ceil(x)
表示大于或等于实数x的最小整
数
floor(x)
表示小于或等于实数x的最大整
数
round(x) 最接近x的整数
最大、最小max([a，b，求最大数
函数c，．．．])
min([a，b，
求最小数
c，．．])
符号函数
sign(x)
5 自定义函数-调用时：“[返回值列]=M文件名（参数列）”
function 返回变量=函数名（输入变量）
注释说明语句段（此部分可有可无）
函数体语句
6．进行函数的复合运算
compose(f,g) 返回值为f(g(y))
compose(f,g,z) 返回值为f(g(z))
compose(f,g,x,.z) 返回值为f(g(z))
compose(f,g,x,y,z) 返回值为f(g(z))
7 因式分解
8 代数式展开
9 合并同类项
10 进行数学式化简
11 进行变量替换
12 进行数学式的转换
调用Maple中数学式的转换命令，调用格式如下：maple(‘Maple 的数学式转换命令’) 即：
13 解方程
solve(’方程’，’变元’)
注：方程的等号用普通的等号： =
14 解不等式
调用maple中解不等式的命令即可，调用形式如下：
具体说，包括以下五种：
15 解不等式组
调用maple中解不等式组的命令即可，调用形式如下：
16 画图
17 求极限
（1）极限：
（2）单侧极限：
左极限：
右极限：
18 求导数
或者：
19 求高阶导数
或者：
20
在
MA TLAB中没有直接求隐函数导数的命令，但是我们可以根据数学
中求隐函数导数的方法，在中一步一步地进行推导；也可以自己编一个求隐函数导数的小程序；不过，最简便的方法是调用Maple 中求隐函数导数的命令，调用格式如下：
maple('implicitdiff(f(x,y)=0,y,x)')
在MATLAB中，没有直接求参数方程确定的函数的导数的命令，只能根据参数方程确定的函数的求导公式
一步一步地进行推导；或者，干脆自己编一个小程序，应用起来会更加方便。

21 求不定积分
int('f(x)')
int ('f(x)','x')
或者：
syms x
22 求定积分、广义积分
或者：
23 进行换元积分的计算
自身没有提供这一功能，但是可以调用Maple函数库中的
changevar命令，调用方法如下：
maple(' with(student)' ) 加载student函数库后，才能使用changevar命令
maple(' changevar( m(x)=p(u), Int(f(x),x) ) ' ) 把积分
表达式中的m(x)代换成p(u)
24 进行分部积分的计算
自身没有提供这一功能，但是可以调用Maple函数库中的
intparts命令，调用方法如下：
maple(' with(student)' ) 加载student函数库后，才能使用intparts命令
maple('intparts(Int(f(x),x),u)' ) 指定u，用分部积分公式进行计算
25 对数列和级数进行求和
syms n
symsum(f(n), n a ,b )
26 进行连乘
27 展开级数
28 进行积分变换
在matlab中，矩形法、梯形法和辛普森法求近似积分
可以用自身的命令，也可调用Maple的相应命令。

调用方法如下：
29 解微分方程
30 解微分方程组。