自適應(yīng)濾波的MATLAB實現(xiàn)

1、自適應(yīng)濾波的MATLAB實現(xiàn)2204090 發(fā)表于 : 2006-12-06 08:56來源 : 中國振動聯(lián)盟1.LMS 算法的仿真程序：%lms 算法clear allclose allhold off% 系統(tǒng)信道權(quán)數(shù)sysorder = 5 ;% 抽頭數(shù)N=1000;% 總采樣次數(shù)inp = randn(N,1);% 產(chǎn)生高斯隨機系列n = randn(N,1);b,a = butter(2,0.25);Gz = tf(b,a,-1);% 逆變換函數(shù)h= 0.0976;0.2873;0.3360;0.2210;0.0964;% 信道特性向量y = lsim(Gz,inp);% 加入噪聲n

2、= n * std(y)/(10*std(n);%噪聲信號d = y + n;% 期望輸出信號totallength=size(d,1);% 步長N=60 ; %60 節(jié)點作為訓(xùn)練序列%算法的開始w = zeros ( sysorder , 1 ) ;% 初始化for n = sysorder : Nu = inp(n:-1:n-sysorder+1) ;% u的矩陣y(n)= w' * u;% 系統(tǒng)輸出e(n) = d(n) - y(n) ;% 誤差if n < 20mu=0.32;elsemu=0.15;endw = w + mu * u * e(n) ;%迭代方程end%檢

3、驗結(jié)果for n = N+1 : totallengthu = inp(n:-1:n-sysorder+1) ;y(n) = w' * u ;e(n) = d(n) - y(n) ;% 誤差endhold onplot(d)plot(y,'r');title(' 系統(tǒng)輸出 ') ;xlabel(' 樣本 ')ylabel(' 實際輸出 ')figuresemilogy(abs(e) ;% e 的絕對值坐標(biāo)title(' 誤差曲線 ') ;xlabel(' 樣本 ')ylabel('

4、誤差矢量 ')figure% 作圖plot(h, 'k+')hold onplot(w, 'r*')legend('實際權(quán)矢量 ','估計權(quán)矢量 ')title(' 比較實際和估計權(quán)矢量') ;axis(0 6 0.05 0.35)2.NLMS 算法的仿真程序：%lms 算法clear allclose allhold off% 系統(tǒng)信道權(quán)數(shù)sysorder = 5 ;% 抽頭數(shù)N=1000;% 總采樣次數(shù)inp = randn(N,1);% 產(chǎn)生高斯隨機系列n = randn(N,1);b,a = but

5、ter(2,0.25);Gz = tf(b,a,-1);% 逆變換函數(shù)h= 0.0976;0.2873;0.3360;0.2210;0.0964;% 信道特性向量y = lsim(Gz,inp);% 加入噪聲n = n * std(y)/(10*std(n);%噪聲信號d = y + n;% 期望輸出信號totallength=size(d,1);% 步長N=60 ; %60 節(jié)點作為訓(xùn)練序列%算法的開始w = zeros ( sysorder , 1 ) ;% 初始化for n = sysorder : Nu = inp(n:-1:n-sysorder+1) ;% u的矩陣y(n)= w&#

6、39; * u;% 系統(tǒng)輸出r(n)=u'*u;% 自相關(guān)矩陣e(n) = d(n) - y(n) ;% 誤差fai=.0001;% 修正參數(shù)，為防止u'*u 過小導(dǎo)致步長值太大而設(shè)置的if n < 20mu=0.32;elsemu=0.15;endw = w + mu * u * e(n)/(r(n)+fai) ;%迭代方程end%檢驗結(jié)果for n = N+1 : totallengthu = inp(n:-1:n-sysorder+1) ;y(n) = w' * u ;e(n) = d(n) - y(n) ;% 誤差endhold onplot(d)plot

7、(y,'r');title(' 系統(tǒng)輸出 ') ;xlabel(' 樣本 ')ylabel(' 實際輸出 ')figuresemilogy(abs(e) ;% e 的絕對值坐標(biāo)title(' 誤差曲線 ') ;xlabel(' 樣本 ')ylabel(' 誤差矢量 ')figure% 作圖plot(h, 'k+')hold onplot(w, 'r*')legend('實際權(quán)矢量 ','估計權(quán)矢量 ')title(

8、9; 比較實際和估計權(quán)矢量') ;axis(0 6 0.05 0.35)3.RLS 算法的仿真程序：% RLS 算法randn('seed', 0) ;rand('seed', 0) ;NoOfData = 8000 ; % Set no of data points used for trainingOrder = 32 ; %自適應(yīng)濾波權(quán)數(shù)Lambda = 0.98 ; %遺忘因子Delta = 0.001 ; %相關(guān)矩陣R 的初始化x = randn(NoOfData, 1) ;% 高斯隨機系列h = rand(Order, 1) ; %系統(tǒng)隨機抽

9、樣d = filter(h, 1, x) ; %期望輸出% RLS 算法的初始化P = Delta * eye ( Order, Order ) ;% 相關(guān)矩陣w = zeros ( Order, 1 ) ;% 濾波系數(shù)矢量的初始化% RLS Adaptationfor n = Order : NoOfData ;u = x(n:-1:n-Order+1) ;% 延時函數(shù)pi_ = u' * P ;% 互相關(guān)函數(shù)k = Lambda + pi_ * u ;K = pi_'/k;% 增益矢量e(n) = d(n) - w' * u ;% 誤差函數(shù)w = w + K * e

10、(n) ;% 遞歸公式PPrime = K * pi_ ;P = ( P - PPrime ) / Lambda ;% 誤差相關(guān)矩陣w_err(n) = norm(h - w) ;% 真實估計誤差end ;% 作圖表示結(jié)果figure ;plot(20*log10(abs(e) ;%| e | 的誤差曲線title(' 學(xué)習(xí)曲線 ') ;xlabel(' 迭代次數(shù) ') ;ylabel(' 輸出誤差估計') ;figure ;semilogy(w_err) ;% 作實際估計誤差圖title(' 矢量估計誤差') ;xlabel(&

11、#39; 迭代次數(shù) ') ;ylabel(' 誤差權(quán)矢量 ') ;4.自適應(yīng)均衡器的仿真程序：% Illustration of the conventional RLS algorithm close all;W=2.9;Nexp=10;N=2000;Nmc=1; % Number of ensemble realizations M=11;% 抽頭系數(shù)lambda=0.99;% 遺忘因子varv=0.001;% 噪聲方差h=zeros(3,1);%h 的初始化er=zeros(N,Nmc);%er 的初始化h(1)=0.5*(1+cos(2*pi*(1-2)/W);

12、h(2)=0.5*(1+cos(2*pi*(2-2)/W);h(3)=0.5*(1+cos(2*pi*(3-2)/W);% 學(xué)習(xí)曲線hc=0 h(1) h(2) h(3)'n0=7;t=(1:N)'for i=1:Nmcy=sign(rand(N,1)-0.5);% 輸入信號v=sqrt(varv)*randn(N,1);% 噪聲信號x=filter(hc,1,y)+v;%信號混合x=zeros(M-1,1);x;%x矩陣yd=zeros(N+M-1,1); % 延遲信號初始化e=yd;yd(n0+M-1:N+M-1)=y(1:N-n0+1);% CRLS 算法% Initia

13、lization lambda=0.98; P=(10-3)*eye(M,M); c=zeros(M,1);g=c;glambda=g;% 迭代范圍for n=M:M+N-1xn=flipud(x(n-M+1:n);glambda=P*xn;alphal=lambda+conj(glambda')*xn;g=glambda/lambda;a(n)=1-conj(g')*xn;P=(P-g*conj(glambda')/lambda;P=(P+P')/2;e(n)=yd(n)-conj(c')*xn;c=c+g*conj(e(n);endendeplot=

14、e(M:M+N-1).2;subplot(2,1,1), plot(t,abs(eplot)ylabel('|e(n)|2');xlabel('n');subplot(2,1,2), plot(t,a(M:M+N-1);ylabel('alpha(n)');xlabel('n');5.自適應(yīng)陷波器的仿真程序：N=400; % 總采樣長度t=0:N-1; % 時間的變化范圍s=sin(2*pi*t/20); % 輸入信號A=0.5; % 干擾信號的幅值fai=pi/3;% 干擾信號的相移n=A*cos(2*pi*t/10+fai);%干擾信號x=s+n;% 信號混合subplot(2,2,1);% 作第一子圖plot(t,s);subplot(2,2,2); % 作第二子圖plot(t,x);x1=cos(2*pi*t/10);x2=sin