MATLAB自适应滤波去噪.docx
- 文档编号:6337405
- 上传时间:2023-01-05
- 格式:DOCX
- 页数:19
- 大小:113.95KB
MATLAB自适应滤波去噪.docx
《MATLAB自适应滤波去噪.docx》由会员分享,可在线阅读,更多相关《MATLAB自适应滤波去噪.docx(19页珍藏版)》请在冰豆网上搜索。
MATLAB自适应滤波去噪
《MATLAB自适应滤波去噪》
课程设计报告
1.课程设计目的
此次课程设计目的是为了让我们学会使用MATLAB进行计算机仿真,使用自适应滤波法设计一个语音去噪声电路。
培养我们的电路设计思路及其算法,明白理论与实践相结合的重要性,培养了我们的实际操作能力以及锻炼我们对实际问题的分析与解决的能力。
2.课程设计内容
2.1LMS自适应算法原理
自适应过程一般采用典型LMS自适应算法,但当滤波器的输入信号为有色随机过程时,特别是当输入信号为高度相关时,这种算法收敛速度要下降许多,这主要是因为输入信号的自相关矩阵特征值的分散程度加剧将导致算法收敛性能的恶化和稳态误差的增大。
此时若采用变换域算法可以增加算法收敛速度。
变换域算法的基本思想是:
先对输入信号进行一次正交变换以去除或衰减其相关性,然后将变换后的信号加到自适应滤波器以实现滤波处理,从而改善相关矩阵的条件数。
因为离散傅立叶变换DFT本身具有近似正交性,加之有FFT快速算法,故频域分块LMSFBLMS算法被广泛应用。
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]=refere[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算法源程序
clearall
closeall
%channelsystemorder
sysorder=5;
%Numberofsystempoints
N=2000;
inp=randn(N,1);
n=randn(N,1);
[b,a]=butter(2,0.25);
Gz=tf(b,a,-1);
%ThisfunctionissubmittedtomakeinverseZ-transform(Matlabcentralfileexchange)
%Thefirstsysorderweightvalue
%h=ldiv(b,a,sysorder)';
%ifyouuseldivthiswillgiveh:
filterweightstobe
h=[0.0976;
0.2873;
0.3360;
0.2210;
0.0964;];
y=lsim(Gz,inp);
%addsomenoise
n=n*std(y)/(10*std(n));
d=y+n;
totallength=size(d,1);
%Take60pointsfortraining
N=60;
%beginofalgorithm
w=zeros(sysorder,1);
forn=sysorder:
N
u=inp(n:
-1:
n-sysorder+1);
y(n)=w'*u;
e(n)=d(n)-y(n);
%Startwithbigmuforspeedingtheconvergencethenslowdowntoreachthecorrectweights
ifn<20
mu=0.32;
else
mu=0.15;
end
w=w+mu*u*e(n);
end
%checkofresults
forn=N+1:
totallength
u=inp(n:
-1:
n-sysorder+1);
y(n)=w'*u;
e(n)=d(n)-y(n);
end
holdon
plot(d)
plot(y,'r');
title('Systemoutput');
xlabel('Samples')
ylabel('Trueandestimatedoutput')
figure
semilogy((abs(e)));
title('Errorcurve');
xlabel('Samples')
ylabel('Errorvalue')
figure
plot(h,'k+')
holdon
plot(w,'r*')
legend('Actualweights','Estimatedweights')
title('Comparisonoftheactualweightsandtheestimatedweights');
axis([060.050.35])
%RLS算法
randn('seed',0);
rand('seed',0);
NoOfData=8000;%Setnoofdatapointsusedfortraining
Order=32;%Settheadaptivefilterorder
Lambda=0.98;%Settheforgettingfactor
Delta=0.001;%RinitializedtoDelta*I
x=randn(NoOfData,1);%Inputassumedtobewhite
h=rand(Order,1);%Systempickedrandomly
d=filter(h,1,x);%Generateoutput(desiredsignal)
%InitializeRLS
P=Delta*eye(Order,Order);
w=zeros(Order,1);
%RLSAdaptation
forn=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;
%Plotresults
figure;
plot(20*log10(abs(e)));
title('LearningCurve');
xlabel('IterationNumber');
ylabel('OutputEstimationErrorindB');
figure;
semilogy(w_err);
title('WeightEstimationError');
xlabel('IterationNumber');
ylabel('WeightErrorindB');
2.3去噪声前后信号波形
Systemoutput
Errorcurve
Comparison
LearningCurve
WeightEstimationError
附录资料:
MATLAB的30个方法
1内部常数
pi
圆周率
exp
(1)
自然对数的底数e
i或j
虚数单位
Inf或inf
无穷大
2数学运算符
a+b
加法
a-b
减法
a*b
矩阵乘法
a.*b
数组乘法
a/b
矩阵右除
a\b
矩阵左除
a./b
数组右除
a.\b
数组左除
a^b
矩阵乘方
a.^b
数组乘方
-a
负号
'
共轭转置
.'
一般转置
3关系运算符
==
等于
<
小于
>
大于
<=
小于或等于
>=
大于或等于
~=
不等于
4常用内部数学函数
指数函数
exp(x)
以e为底数
对数函数
log(x)
自然对数,即以e为底数的对数
log10(x)
常用对数,即以10为底数的对数
log2(x)
以2为底数的x的对数
开方函数
sqrt(x)
表示x的算术平方根
绝对值函数
abs(x)
表示实数的绝对值以及复数的模
三角函数
(自变量的单位为弧度)
sin(x)
正弦函数
cos(x)
余弦函数
tan(x)
正切函数
cot(x)
余切函数
sec(x)
正割函数
csc(x)
余割函数
反三角函数
asin(x)
反正弦函数
acos(x)
反余弦函数
atan(x)
反正切函数
acot(x)
反余切函数
asec(x)
反正割函数
acsc(x)
反余割函数
双曲函数
sinh(x)
双曲正弦函数
cosh(x)
双曲余弦函数
tanh(x)
双曲正切函数
coth(x)
双曲余切函数
sech(x)
双曲正割函数
csch(x)
双曲余割函数
反双曲函数
asinh(x)
反双曲正弦函数
acosh(x)
反双曲余弦函数
atanh(x)
反双曲正切函数
acoth(x)
反双曲余切函数
asech(x)
反双曲正割函数
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因式分解
syms表达式中包含的变量
factor(表达式)
8代数式展开
syms表达式中包含的变量
expand(表达式)
9合并同类项
syms表达式中包含的变量
collect(表达式,指定的变量)
10进行数学式化简
syms表达式中包含的变量
simplify(表达式)
11进行变量替换
syms表达式和代换式中包含的所有变量
subs(表达式,要替换的变量或式子,代换式)
12进行数学式的转换
调用Maple中数学式的转换命令,调用格式如下:
maple(‘Maple的数学式转换命令’) 即:
maple(‘convert(表达式,form)’’)将表达式转换成form的表示方式
maple(‘convert(表达式,form,x)’)指定变量为x,将依赖于变量
x的函数转换成form的表示方式(此指令仅对form为exp与sincos的转换式有用)
13解方程
solve(’方程’,’变元’)
注:
方程的等号用普通的等号:
=
14解不等式
调用maple中解不等式的命令即可,调用形式如下:
maple('maple中解不等式的命令')
具体说,包括以下五种:
maple('solve(不等式)')
maple('solve(不等式,变元)')
maple('solve({不等式},变元)')
maple('solve(不等式,{变元})')
maple('solve({不等式},{变元})')
15解不等式组
调用maple中解不等式组的命令即可,调用形式如下:
maple('maple中解不等式组的命令')
即:
maple('solve({不等式组},{变元组})')
16 画图
方法1:
先产生横坐标x的取值和相应的纵坐标y的取值,然后执行命令:
plot(x,y)
方法2:
fplot('f(x)',[xmin,xmax])
fplot('f(x)',[xmin,xmax,ymin,ymax])
方法3:
ezplot('f(x)')
ezplot('f(x)',[xmin,xmax])
ezplot('f(x)',[xmin,xmax,ymin,ymax])
17求极限
(1) 极限:
symsx
limit(f(x),x,a)
(2)单侧极限:
左极限:
symsx
limit(f(x),x,a,’left’)
右极限:
symsx
limit(f(x),x,a,’right’)
18求导数
diff('f(x)')
diff('f(x)','x')
或者:
symsx
diff(f(x))
symsx
diff(f(x),x)
19求高阶导数
diff('f(x)',n)
diff('f(x)','x',n)
或者:
symsx
diff(f(x),n)
symsx
diff(f(x),x,n)
20在MATLAB中没有直接求隐函数导数的命令,但是我们可以根据数学中求隐函数导数的方法,在中一步一步地进行推导;也可以自己编一个求隐函数导数的小程序;不过,最简便的方法是调用Maple中求隐函数导数的命令,调用格式如下:
maple('implicitdiff(f(x,y)=0,y,x)')
在MATLAB中,没有直接求参数方程确定的函数的导数的命令,只能根据参数方程确定的函数的求导公式
一步一步地进行推导;或者,干脆自己编一个小程序,应用起来会更加方便。
21求不定积分
int('f(x)')
int('f(x)','x')
或者:
symsx
int(f(x))
symsx
int(f(x),x)
22求定积分、广义积分
int('f(x)',a,b)
int('f(x)','x',a,b)
或者:
symsx
int(f(x),a,b)
symsx
int(f(x),x,a,b)
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对数列和级数进行求和
symsn
symsum(f(n),na,b)
26进行连乘
maple('product(f(n),n=a..b)')
27展开级数
symsx
taylor(f(x),x,n,a)
28进行积分变换
symsst
laplace(f(t),t,s)拉普拉斯变换
ilaplace(F(s),s,t)拉普拉斯变换的逆变换
symstω
fourier(f(t),t,ω)傅立叶变换
ifourier(F(ω),ω,t)傅立叶变换的逆变换
symsnz
ztrans(f(n),n,z)Z变换
iztrans(F(z),z,n)Z变换的逆变换
在matlab中,矩形法、梯形法和辛普森法求近似积分
可以用自身的命令,也可调用Maple的相应命令。
调用方法如下:
maple('with(student)')
maple('Maple中求定积分近似值的命令')
29解微分方程
dsolve('微分方程','自变量')
dsolve('微分方程','初始条件或边界条件','自变量')
30解微分方程组
dsolve('微分方程组','自变量')
dsolve('微分方程组','初始条件或边界条件','自变量')
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 自适应 滤波