MATLAB自適應(yīng)濾波去噪

1、MATLAB自適應(yīng)濾波去噪課程設(shè)計(jì)報(bào)告1. 課程設(shè)計(jì)目的此次課程設(shè)計(jì)目的是為了讓我們學(xué)會(huì)使用 MATLA進(jìn)行計(jì)算機(jī)仿真，使用自適應(yīng) 濾波法設(shè)計(jì)一個(gè)語(yǔ)音去噪聲電路。培養(yǎng)我們的電路設(shè)計(jì)思路及其算法，明白理論和實(shí) 踐相結(jié)合的重要性，培養(yǎng)了我們的實(shí)際操作能力以及鍛煉我們對(duì)實(shí)際問(wèn)題的分析和解 決的能力。2. 課程設(shè)計(jì)內(nèi)容2.1 LMS自適應(yīng)算法原理自適應(yīng)過(guò)程一般采用典型LMS自適應(yīng)算法，但當(dāng)濾波器的輸入信號(hào)為有色隨機(jī)過(guò) 程時(shí)，特別是當(dāng)輸入信號(hào)為高度相關(guān)時(shí)，這種算法收斂速度要下降許多，這主要是因 為輸入信號(hào)的自相關(guān)矩陣特征值的分散程度加劇將導(dǎo)致算法收斂性能的惡化和穩(wěn)態(tài) 誤差的增大。此時(shí)若采用變換域算法可以

2、增加算法收斂速度。變換域算法的基本思想 是：先對(duì)輸入信號(hào)進(jìn)行一次正交變換以去除或衰減其相關(guān)性，然后將變換后的信號(hào)加 到自適應(yīng)濾波器以實(shí)現(xiàn)濾波處理，從而改善相關(guān)矩陣的條件數(shù)。因?yàn)殡x散傅立葉變換 DFT本身具有近似正交性，加之有 FFT快速算法，故頻域分塊LMS FBLMS算法 被廣泛使用。FBLMSS法本質(zhì)上是以頻域來(lái)實(shí)現(xiàn)時(shí)域分塊 LMS算法的，即將時(shí)域數(shù)據(jù)分組構(gòu)成 N個(gè)點(diǎn)的數(shù)據(jù)塊，且在每塊上濾波權(quán)系數(shù)保持不變。其原理框圖如圖2所示。FBLMS算法在頻域內(nèi)可以用數(shù)字信號(hào)處理中的重疊保留法來(lái)實(shí)現(xiàn)，其計(jì)算量比時(shí)域法大為減少，也可以用重疊相加法來(lái)計(jì)算，但這種算法比重疊保留法需要較大的計(jì)算量。塊數(shù) 據(jù)的

3、任何重疊比例都是可行的，但以50%的重疊計(jì)算效率為最高。對(duì)FBLM算法和典 型LMS算法的運(yùn)算量做了比較，并從理論上討論了兩個(gè)算法中乘法部分的運(yùn)算量。本 文從實(shí)際工程出發(fā)，詳細(xì)分析了兩個(gè)算法中乘法和加法的總運(yùn)算量，其結(jié)果為：復(fù)雜度之比二FBLM麒數(shù)乘加次數(shù)/ LMS實(shí)數(shù)乘加次數(shù)=(25Nlog2N+ 2N 4)/ 2N(2N 1)采用ADSP勺C語(yǔ)言來(lái)實(shí)現(xiàn)FBLMSJ法的程序如下:for(i = 0;i v = 30;i + +)for(j 二0;j v= n 1;j + +)inj = inputi x N+ j;rfft(i n,ti n,n f,wfft,wst, n);rfft(w,t

4、w,wf,wfft,wst, n);cvecvmlt(i nf,wf,i nw, n);ifft(i nw,t,O,wfft,wst, n);for(j = 0,j v = N- 1;j + + )yi x N+ j = ON+ j . re;ei x N+ j = refereix N+ j yi x N+ j;tempN + j = ei x N+ j;rfft(temp,t,E,wfft,wst, n);for(j = 0;j v = n - 1;j + + )in f_conjj= conjf(in fj);cvecvmlt(E,i nf_conj,Ei n,n);ifft(Ei n,

5、t,Ei n,wfft,wst, n);for(j = 0;j v = N- 1;j + + )OOj二 Einj . re;wj二 wj + 2* u* OOj;在EZ KIT測(cè)試板中，筆者用匯編語(yǔ)言和C語(yǔ)言程序分別測(cè)試了典型LMS算法的 運(yùn)行速度，并和FBLMSf法的C語(yǔ)言運(yùn)行速度進(jìn)行了比較，表 2所列是其比較結(jié)果， 從表2可以看出濾波器階數(shù)為64時(shí)，即使是用C語(yǔ)言編寫(xiě)的FBLMSJ法也比用匯編 編寫(xiě)的LMS算法速度快20%以上，如果濾波器的階數(shù)更大，則速度會(huì)提高更多。2.2語(yǔ)音信號(hào)去噪聲源程序%lms算法源程序clear allclose all%cha nnel system orde

6、rsysorder = 5 ;% Number of system pointsN=2000;inp = randn( N,1);n = randn( N,1);b,a = butter(2,0.25);Gz = tf(b,a,-1);%This fun ctio n is submitted to make inv erse Z-tra nsform (Matlab cen tral file excha nge)%The first sysorder weight value%h=ldiv(b,a,sysorder)'% if you use ldiv this will give

7、 h :filter weights to beh= 0.0976;0.2873;0.3360;0.2210;0.0964;y = lsim(G z,in p);%add some no isen = n * std(y)/(10*std(n);d = y + n;totalle ngth=size(d,1);%Take 60 poi nts for trai ningN=60 ;%beg in of algorithmw = zeros ( sysorder , 1 );for n = sysorder : Nu = inp(n :-1: n-sysorder+1);y(n )= w'

8、; * u;e(n) = d(n) - y(n);% Start with big mu for speeding the convergenee then slow down to reach the correct weightsif n < 20mu=0.32;elsemu=0.15;endw = w + mu * u * e(n);end%check of resultsfor n = N+1 : totalle ngthu = inp(n :-1: n-sysorder+1);y( n) = w' * u ;e( n) = d(n) - y( n);endhold on

9、plot(d)plot(y,'r');title('System output');xlabel('Samples')ylabel('True and estimated output') figuresemilogy(abs(e);title('Error curve'); xlabel('Samples')ylabel('Error value')figureplot(h, 'k+')hold onplot(w, 'r*')lege nd('

10、;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);ran d('seed', 0);NoOfData = 8000 ;% Set no of data poi nts used for trainingOrder = 32 ;% Set the adaptive filter

11、orderLambda = 0.98 ;% Set the forgetti ng factorDelta = 0.001 ;% R in itialized to Delta*Ix = randn(No OfData, 1) ;% In put assumed to be whiteh = ran d(Order, 1) ;% System picked ran domlyd = filter(h, 1, x) ;% Gen erate output (desired sig nal)% In itialize RLSP = Delta * eye ( Order, Order );w =

12、zeros ( Order, 1 );% RLS Adaptati onfor 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) = no rm(h - w);end ;% Plot resultsfigure ;plot(20*log10(abs(e)

13、;title('Lear ning Curve');xlabel('lteration 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去噪聲前后信號(hào)波形System outputSystem output2004006008001000Sampl!e312001400160018002G00o5"5嘗 p£bui-¥ puE«un 匚Error curveComparis onComparison of the actual weights grid the astimated weights1II+1+ Actual weights 十 Estimated eights0.3+*-0254-+02-0 15*0 10.05+II+ i123465Lear ning CurveLearning! Curve10002