基于LMS算法的多麦克风降噪

武汉理工大学《信息处理课群综合训练》课程设计说明书绪论本次课程设计是基于LMS算法的多麦克风降噪，多麦克风降噪是降噪技术的一个重要应用。我国的降噪技术采用的手段主要有三种，其中的动态降噪技术(DNR)又可以分为模拟动态降噪技术和数字动态降噪技术。目前国内外解决噪声问题最普遍的方法是采用模拟动态降噪技术，数字降噪技术的研究尚处于初期阶段。数字降噪技术比模拟降噪技术具有更大的优点。模拟降噪技术全采用硬件实施，修改和调试十分困难，对元器件参数的变化也很敏感，技术指标受元器件的误差影响较大，降噪效果不稳定，不利于产品的批量生产。而数字降噪技术由于采用计算机技术实现自适应滤波，通过修改软件算法就可以达到不同的降噪效果，不用更改硬件结构，调试和维修都非常方便；数字降噪技术采用自适应滤波技术，可以实时跟踪噪声的变化进一步进行处理，因此降噪效果较好。另外，数字降噪技术抗干扰能力强，本身具有自恢复能力，并且在整个音频带内降噪比较均衡，而模拟降噪技术偏重于低频段，高频段效果较差。因此降噪技术未来的发展方向是数字降噪技术，以数字信号处理(DSP)及其相关算法为技术支撑的数字降噪技术代表着当今降噪技术的发展。目前市场上的麦克风降噪产品主要是模拟降噪，因此数字降噪的设计在国内属于领先技术。多麦克风数字降噪的系统原理是通过麦克风装置直接检测出噪声信号和音频信号的混合信号，然后将混合信号通过DSP数字降噪模块进行噪声分离并产生降噪信号来抵消噪声，因此人耳就可以只听到较纯净的音频信号而不受环境噪声的干扰。本文采用最小均方误差(LMS)算法，实现了数字降噪DSP中消除噪声的模块自适应滤波器的设计，介绍了其在MATLAB中编程及仿真输出，并通过程序实现了设计。2语音信号的输入2.1语音信号的采样相关原理2.1.1采样频率和采样位数采样频率是指计算机每秒钟采集多少个声音样本，是描述声音文件的音质、音调，衡量声卡、声音文件的质量标准。采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。采样频率与声音频率之间有一定的关系，根据奈奎斯特理论，只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。这就是说采样频率是衡量声卡采集、记录和还原声音文件的质量标准。采样位数即采样值或取样值，用来衡量声音波动变化的参数，是指声卡在采集和播放声音文件时所使用数字声音信号的二进制位数。采样位数和采样频率对于音频接口来说是最为重要的两个指标，也是选择音频接口的两个重要标准。无论采样频率如何，理论上来说采样的位数决定了音频数据最大的力度范围。每增加一个采样位数相当于力度范围增加了6dB。采样位数越多则捕捉到的信号越精确。对于采样率来说你可以想象它类似于一个照相机，44.1kHz意味着音频流进入计算机时计算机每秒会对其拍照达441000次。显然采样频率越高，计算机摄取的图片越多，对于原始音频的还原也越加精确。2.1.2采样定理在进行模拟/数字信号的转换过程中，当采样频率fsmax大于等于信号最高频率fmax的2倍时，即：fsmax>=2fmax,则采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的5～10倍，采样定理又称奈奎斯特定理。1924年奈奎斯特(Nyquist)就推导出：在理想低通信道的最大码元传输速率=2W*log2N(其中W是理想低通信道的带宽,N是电平强度)。2.2语音信号的采集语音信号的采集可以利用PC机上的声卡，WINDOWS操作系统可以进行数字信号的采集。本次设计下载了一段.wav格式的录音当做语音输入信号。并命名为”OriSound”2.3语音信号的时频分析在MATLAB软件平台下，利用wavread函数对语音信号进行采样，记住采样频率和采样点数，Wavread函数调用格式:y=wavread(file)%读取file所规定的wav文件，返回采样值放在向量y中[y,fs,bits]=wavread(file)%采样值放在向量y中，fs表示采样频率（hz），bits表示采样位数y=wavread（file，N）%读取前N点的采样值放在向量y中y=wavread（file，[N1,N2]）%读取从N1到N2点的采样值放在向量y中对语音信号OriSound.wav进行采样其程序如下：[y,fs,nbits]=wavread(‘OriSound’);%语音信号加载入MATLAB仿真软件平台中画出语音信号的时域波形，再对语音信号进行频谱分析。MATLAB提供了快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：Xk=fft(xn,N)（式2-3-1）参数xn为被变换的时域序列向量，N是DFT变换区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。，当N小于xn的长度时，fft函数计算xn的前N个元素，忽略其后面的元素。在本次设计中，我们利用fft对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。程序如下：[y,fs,nbits]=wavread(‘OriSound’);%语音信号的采集sound(y,fs,nbits);%语音信号的播放n=length(y);Y=fft(y,n);%快速傅里叶变换figure;subplot(2,1,1);plot(y);title(‘原始信号波形’,’fontweight’,’bold’);axis([7800080000-11]);grid;subplot(2,1,2);plot(abs(Y));(‘原始信号频谱’,’fontweight’,’bold’);axis([015000004000]);grid;程序结果如下图：图2-1原始信号采集波形图2.4语音信号加噪与频谱分析在MATLAB中产生高斯白噪声非常方便，我们可以直接应用两个函数：一个是WGN，另一个是AWGN。WGN用于产生高斯白噪声，AWGN则用于在某一信号中加入高斯白噪声。也可直接用randn函数产生高斯分布序列。在本次设计中，我们是利用MATLAB中的随机函数(rand或randn)产生噪声加入到语音信号中，模仿语音信号被污染，并对其频谱分析。Randn函数有两种基本调用格式：Randn(n)和Randn（m,n）,前者产生n×n服从标准高斯分布的随机数矩阵，后者产生m×n的随机数矩阵。在这里，我们选用Randn（m,n）函数。语音信号添加噪声及其频谱分析的主要程序如下：[y,fs,nbits]=wavread('OriSound');sound(y,fs,nbits);n=length(y);Noise=0.2*randn(n,2);s=y+Noise;sound(s);figure;subplot(2,1,1);plot(s);title('加噪语音信号的时域波形','fontweight','bold');axis([7800080000-11]);grid;S=fft(s);subplot(2,1,2);plot(abs(S));title('加噪语音信号的频域波形','fontweight','bold');axis([015000004000]);grid;程序结果如下图：图图2-2信号加噪时域波形图与频谱图3基于LMS自适应滤波器的设计在实际应用中，常常无法得到信号和噪声统计特性的先验知识。在这种情况下，自适应滤波技术能够获得极佳的滤波性能，因而具有很好的应用价值。常用的自适应滤波技术有：最小均方（LMS）自适应滤波器、递推最小二乘（RLS）滤波器、格型滤波器和无限冲激响应（IIR）滤波器等。这些自适应滤波技术的应用又包括：自适应噪声抵消、自适应谱线增强和陷波等。LMS自适应滤波器是使滤波器的输出信号与期望响应之间的误差的均方值为最小，因此称为最小均方误差（LMS）自适应滤波器。3.1LMS算法构成自适应数字滤波器的基本部件是自适应线性组合器，如图3-1的所示。设线性组合器的M个输入为x(k-1)…,x(k-M)，其输出y(k)是这些输入加权后的线性组合，即y(k)=Wx(k-i)（式3-1-1）图3-1自适应线性组合器定义权向量，W=[W1,W2,W3,…,Wm]，且X(k)=[X((k-1)T,…,X((k-M)T)]（式3-1-2）在图3-1中，令d(k)代表“所期望的响应”，并定义误差信号ε(k)=d(k)-y(k)=d(k)-WX(k-i)（式3-1-3）式（3-1-3）写成向量形式ε(k)=d(k)-WX(k)=d(k)-X(k)W（式3-1-4）误差平方为ε(k)=d(k)-2d(k)X(k)W+WX(k)X(k)W上式两边取数学期望后，得均方误差E{ε(k)}=E{d(k)}-2E{d(k)X(k)}W+WE{X(k)X(k)}W（式3-1-5）定义互相关函数行向量R：R=E{d(k)X(k)}（式3-1-6）和自相关函数矩阵R=E{X(k)X(k)}（式3-1-7）则均方误差（3-1-5）式可表述为E{ε(k)}=E{d(k)}-2RW+WRW（式3-1-8）这表明，均方误差是权系数向量W的二次函数，它是一个中间向上凹的抛物形曲面，是具有唯一最小值的函数。调节权系数使均方误差为最小，相当于沿抛物形曲面下降找最小值。可以用梯度来求该最小值。将式（3-1-8）对权系数W求导数，得到均方误差函数的梯度(k)=-2R+2RW(式3-1-9)令(k)=0，即可求出最佳权系数向量W=RR（式3-1-10）它恰好是研究Wiener滤波器遇到过的Wiener-Hopf方程。因此，最佳权系数向量通常也叫作Wiener权系数向量。将W代入式（3-1-8）得最小均方误差E{ε(k)}=E{d(k)}-RW（式3-1-11）利用式（3-1-10）求最佳权系数向量的精确解需要知道R和R的先验统计知识，而且还需要进行矩阵求逆等运算。WidrowandHoff(1960)提出了一种在这些先验统计知识未知时求W的近似值的方法，习惯上称为WidrowandHoffLMS算法。这种算法的根据是最优化方法中的最速下降法。根据最速下降法，“下一时刻”权系数向量W(k+1)应该等于“现时刻”权系数向量W(k)加上一个负均方误差梯度−(k)的比例项，即W(k+1)=W(k)-μ(k)（式3-1-12）式中，μ是一个控制收敛速度与稳定性的常数，称之为收敛因子。不难看出，LMS算法有两个关键：梯度(k)的计算以及收敛因子μ的选择。（一）(k)的近似计算精确计算梯度(k)是十分困难的，一种粗略的但是却十分有效的计算(k)的近似方法是：直接取ε(k)作为均方误差E{ε(k)}的估计值，即(k)=[ε(k)]=2ε(k)[ε(k)]（式3-1-13）得到梯度估值(k)=-2ε(k)X(k)于是，Widrow–HoffLMS算法最终为W(k+1)=W(k)+2με(k)X(k)(式3-1-14）式（3-1-14）的实现方框图如图3-2所示图3-2LMS算法的实现方框图下面分析梯度估值(k)的无偏性。(k)的数学期望为（式3-1-15）在上面的推导过程中，利用了d(k)和ε(k)二者皆为标量的事实。在得到最后的结果时，利用了式（3-1-9）。式（3-1-15）表明，梯度估值(k)是无偏估（二）μ的选择对权系数向量更新公式（3-1-14）两边取数学期望，得E{W(k+1)}=E{W(k)}+2μE{ε(k)X(k)}=(I-2μR)E{W(k)}+2μR（式3-1-16）式中，I为单位矩阵，R=E{d(k)X(k)}和R=E{X(k)X(k)}。当时，k=0时，E{W(1)}=(I-2μR)E{W(0)}+2μR对于k=1，利用上式结果，则有E{W(2)}=(I-2μR)E{W(1)}+2μR(I-2μR)E{W(0)}+2μ(I-2μR)R起始时，E{W(0)}=W(0)故重复以上迭代至k+1，则有E{W(k+1)}=(I-2μR)W(0)+2μ(I-2μR)R（式3-1-17)由于R是实值的对称阵，我们可以写出其特征值分解式R=QΣQ=QΣQ（式3-1-18）这里，我们利用了正定阵Q的性质Q=Q，且Σ=diag(λ,…λ)是对角阵，其对角元素λi是R的特征值。将式（3-1-18）代入式（3-1-19）后得E{W(k+1)}=(I-2μQΣQ)W+2μ(I-2μQΣQ)R（式3-1-19）注意到以下恒等式及关系式：（1）(I-2μQΣQ)=Q(I-2μΣ)Q（2）(I-2μQΣQ)=Q[(2μΣ)]Q（3）假定所有的对角元素的值均小于1（这可以通过适当选择μ实现），则(I-2μΣ)=0（4）R=QΣQ将上式代入式（8-1-19），结果有E{W(k+1)}=QΣQR=RR=W（式3-1-20）由此可见，当迭代次数无限增加时，权系数向量的数学期望值可收敛至Wiener解，其条件是对角阵(I-2μΣ)的所有对角元素均小于1，即0<μ<(式3-1-21)其中λmax是R的最大特征值。μ称为收敛因子，它决定达到式（3-1-20）的速率。事实上，W(k)收敛于W由比值d=λmax/λmin决定，该比值叫做谱动态范围。大的d值喻示要花费很长的时间才会收敛到最佳权值。克服这一困难的方法之一是产生正交数据。基本LMS自适应算法如下：初始化：W(0)=0;R(0)=I;选择μ：0<μ<Fork=1tonfinaldo:W(k)=W(k-1)+2μ[x(k)-W(k-1)X(k)]X(k)LMS自适应滤波器如图3-3所示:图3-3LMS自适应滤波器3.2自适应噪声抵消原理自适应噪声抵消的目的是要去除主信号中的背景噪声。主信号由有用信号和背景噪声组成，而背景噪声与参考信号中的噪声相关。因此，自适应噪声抵消技术主要依赖于从主信号和噪声中获取参考信号。Widrow和Hoff发展了最小均方误差（LMS）自适应算法和称为自适应线性阈值逻辑单元（ADALINE）的模式识别方法。1965年，基于最小均方误差准则(LMS)的自适应噪声抵消首次得以实现，随后，自适应噪声抵消在信号处理、地震和生物医学领域均获得成功应用。基于维纳理论的自适应噪声抵消需要无限加权滤波器，以极小化输出误差。为了实现维纳滤波方案，必须使用有限加权滤波器。换句话说，自适应滤波器必须假定维纳滤波器是一个有限冲激响应（FIR）滤波器。图3-4自适应噪声抵消原理方框图如图3-4（a）所示是基于维纳滤波器的自适应噪声抵消原理方框图。主信号由有用信号x(n)和背景噪声v(n)构成，其中s(n)和v(n)不相关。参考信号r(n)可与s(n)或v(n)相关。(n)是背景噪声的最佳估计。(n)可以通过选择最佳FIR维纳滤波器的最佳加权(n)计算得出，即(n)=(n)r(n-m)0≤m≤M(式3-2-1)其中，M表示滤波器的阶；r(n-m)由延时获得。具有M个权重滤波器的估计误差e(n)由下式定义：e(n)=x(n)-(n)=x(n)-(n)(n)（式3-2-2）由正交原理有，e(n)和r(n)正交。对式（3–2-2）两边取平方和数学期望，可得E[e(n)]=E[x(n)]-2P+R（式3-2-3）其中，输入信号s(n)和参考矢量(n)之间的互相关用表示，即=E[x(n)(n)]（式8-2-4）R表示输入自相关矩阵，即R=[(n)(n)]（式3-2-5）令均方估计误差函数的梯度等于0，可得最佳FIR滤波器（维纳滤波器）权重如下，=R（式3-2-6）实际上，通常和R的统计量是未知的。然而，用Widrow和Hoff提出的方法迭代求解式（3-2-6）能够克服这一限制。如果参考信号和主信号中的噪声相关，则自适应滤波器将在输出端去除其相关性，具体方法是：从参考信道的噪声中产生一个主信道中背景噪声的估计值(n)，然后从主信道中减去这个估计噪声(n)，那么自适应滤波器的输出就是有用信号的估计(n)。用最速下降法（或梯度下降法）可得到式（3-2-6）的解。自适应滤波器的加权值被更新的第(n+1)步迭代式为(n+1)=(n)+Δ(n)=(n)-（式3-2-7）3.3基于最小均方误差准则(LMS)的自适应噪声抵消如果没有关于参考信号向量P和输入自相关矩阵R的先验信息，要实现最优滤波器加权是不可能的。因此，Widrow和Hoff提出了另一种可迭代的维纳FIR滤波实现方法。在这种方法中，滤波器的权重可被更新为：(n+1)=(n)+2μe(n)(n)（式3-3-1）综上所述，基于最小均方误差准则(LMS)的自适应噪声抵消算法可按以下步骤实现：第一步：设一个初值；(0)；第二步：计算自适应FIR滤波器的输出(n)，(n)=(n)r(n-m)（式3-3-2）其中，M表示滤波器的阶。第三步：估计当前时刻n的误差e(n),e(n)=x(n)-(n)≈(n)（式3-3-3）第四步：用最速下降LMS算法更新滤波器权重：(n):(n+1)=(n)+2μe(n)r(n-m)0（式3-3-4）第五步：校验误差是否满足标准。若满足，则停止迭代，否则进行下一步；第六步：nn+1，到下一个时刻，重复以上步骤，直至满足要求为止。收敛参数μ必须是正数，并且满足：0<μ<（式3-3-5）0<μ<（式3-3-6）其中，λmax表示自相关矩阵R的最大特征值。然而，在实际应用中，R的具体值是不知道的，参数μ的值也需要试探性地选择。若μ取值小，能保证收敛，但需要注意的是，如果取得过小，收敛速度将非常慢；相反，若μ取值大，可以提高收敛速度，却是以噪声收敛为代价的。如果参考输入信号r(n)是频率为ω的正弦信号，自适应滤波器将从主信号中滤除所有的频率为ω的正弦成分。基于最小均方误差准则(LMS)的自适应噪声抵消算法的程序见下节。3.4LMS算法程序function[yn,W,en]=LMS(xn,dn,M,mu,itr)LMS(LeastMeanSqure)算法%输入参数:%xn输入的信号序列(列向量)%dn所期望的响应序列(列向量)%M滤波器的阶数(标量)%mu收敛因子(步长)(标量)要求大于0,小于xn的相关矩阵最大特征值的倒数%itr迭代次数(标量)默认为xn的长度,M<itr<length(xn)%输出参数:%W滤波器的权值矩阵(矩阵)大小为Mxitr,%en误差序列(itrx1)(列向量)%yn实际输出序列(列向量)%参数个数必须为4个或5个ifnargin==4%4个时递归迭代的次数为xn的长度itr=length(xn);elseifnargin==5%5个时满足M<itr<length(xn)ifitr>length(xn)|itr<Merror('迭代次数过大或过小!');endelseerror('请检查输入参数的个数!');end%初始化参数en=zeros(itr,1);%误差序列,en(k)表示第k次迭代时预期输出与实际输入的误差W=zeros(M,itr);%每一行代表一个加权参量,每一列代表-次迭代,初始为0%迭代计算fork=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));fork=M:length(xn)x=xn(k:-1:k-M+1);yn(k)=W(:,end).'*x;end4滤波前后信号的波形及频谱4.1验证所设计的自适应滤波器为了验证滤波器的可使用性，我们用常用的sin函数来进行验证。其具体程序及运行结果如下：调用LMS算法：%functionmain()closeall%周期信号的产生t=0:99;xs=10*sin(0.5*t);figure;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;%滤波器的阶数rho_max=max(eig(xn*xn.'));%输入信号相关矩阵的最大特征值mu=rand()*(1/rho_max);%收敛因子0<mu<1/rho[yn,W,en]=LMS(xn,dn,M,mu);%绘制滤波器输入信号figure;subplot(2,1,1);plot(t,xn);grid;ylabel('幅值');xlabel('时间');title('it{滤波器输入信号}');%绘制自适应滤波器输出信号subplot(2,1,2);plot(t,yn);grid;ylabel('幅值');xlabel('时间');title('it{自适应滤波器输出信号}');%绘制自适应滤波器输出信号,预期输出信号和两者的误差figureplot(t,yn,'b',t,dn,'g',t,dn-yn,'r');grid;legend('自适应滤波器输出','预期输出','误差');ylabel('幅值');xlabel('时间');title('it{自适应滤波器}');运行结果：图图4-1验证自适应滤波器可行性1图4-2验证自适应滤波器可行性2图4-3验证自适应滤波器可行性3由所得结果可知，所设计的滤波器符合要求。4.2对主麦克风音频信号滤波4.2.1程序流程图图4-4程序流程图4.2.2LMS自适应滤波用自己设计的滤波器对加噪的主麦克风语音信号进行滤波，主麦克风录制的带噪声的语音信号为LMSprimsp.wav，参考麦克风录制的噪声语音信号为LMSrefns.wav。1）主麦克风录制的语音信号的读取[primary,fs,nbits]=wavread('LMSprimsp.wav');sound(primary,fs,nbits);%语音信号的播放2）绘制主麦克风录制的语音信号波形及频谱n=length(primary);P=fft(primary,n);%快速傅里叶变换figure;subplot(2,1,1);plot(primary);ylabel('幅值');xlabel('时间');title('主麦克风录制的受噪声污染的语音信号','fontweight','bold');grid;subplot(2,1,2);plot(abs(P));title('主麦克风语音信号频谱','fontweight','bold');grid;运行结果：图图4-5主麦克风录制的语音信号3）参考麦克风录制的噪声信号的读取，并绘制其波形及频谱fref=wavread('LMSrefns.wav');n=length(fref);F=fft(fref,n);figure;subplot(2,1,1);plot(fref);grid;ylabel('幅值');xlabel('时间');title('参考麦克风录制的噪声信号','fontweight','bold');subplot(2,1,2);plot(abs(F));title('噪声信号频谱','fontweight','bold');grid;运行结果：图4-6参考麦克风录制的噪声信号4）调用LMS自适应滤波函数对主麦克风录制的语音信号进行滤波，得到增强后的语音信号。xs=primary-fref;xn=primary.';%输入信号序列dn=xs.';%预期结果序列M=20;%滤波器阶数rho_max=max(eig(xn*xn.'));%输入信号相关矩阵的最大特征值mu=rand()*(1/rho_max);%收敛因子0<mu<1/rho[yn,W,en]=LMS(xn,dn,M,mu);5）算法仿真收敛后，得到增强的语音信号，绘制滤波增强后的语音信号，并分别对增强前后的语音信号作频谱分析。figure;subplot(2,1,1);plot(primary);axis([1length(primary)min(primary)max(primary)]);title('输入主麦克风信号');subplot(2,1,2);plot(t,yn);grid;ylabel('幅值');xlabel('时间');title('自适应滤波器输出信号');%绘制滤波增强后信号频谱n=length(yn);Y=fft(yn,n);figuresubplot(2,1,2);plot(abs(Y));grid;title('自适应滤波器输出信号频谱','fontweight','bold');运行结果：图4-6滤波增强前后信号频谱7）回放增强后的语音信号wavplay(yn)

；

%语音输出降噪后的语音wavwrite(yn,’output’)；

%生成语音文件output.wav文件4.3调试分析通过对滤波前后主麦克风语音信号波形的分析可知，滤波前信号波形受噪声干扰严重，有用信号被噪声覆盖，几乎都是噪声；滤波后滤除了信号中的噪声分量，有用信号被提取出来，起到了信号增强的目的。通过对滤波前后主麦克风语音信号频谱的观察分析可知，滤波前后频谱分量明显减弱，滤波后频谱分量较集中，有明显改善。经过以上的降噪处理后，可在Matlab中用函数sound对声音进行回放。其调用格式：sound(y,Fs)，sound(y)和sound(y,Fs,bits)。可以察觉滤波前后的声音有明显的变化。对主麦克风输入的受噪声干扰的语音信号，进行基于LMS自适应滤波的源程序代码见附录。心得体会参考文献[1]刘泉，江雪梅.《信号与系统》.高等教育出版社出版社，2005[2]刘泉，阙大顺，郭志强.《数字信号处理原理与实现（第2版）》.电子工业出版社，2009[3]樊昌信，曹丽娜.《通信原理（第6版）》.国防工业出版社，2006[4]李正周.《MATLAB数字信号处理与应用》.清华大学出版社，2008[5]徐明远，刘增力.《MATLAB仿真在信号处理中的应用》.电子科技大学出版社，2007附录程序代码：source=wavread('input.wav');%读入wav格式的音频文件noise=0.19*randn(length(source),1);source1=source+noise;wa