2007-07-30 02:16:01 

關于自適應濾波的，LMS和RLS。
 
clear;
close all;
N=1500;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        以下利用同余法生成白噪的過程              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=6;
x_0=1;
M=255;
f=2;
for i=1:N
 x_A=A*x_0;
 xi=mod(x_A,M);
 vi=xi/256;
 v(i)=(vi-0.5)*f;
 x_0=xi;
end                                          %噪聲v

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        以下利用System Identification Toolbox     %
%       工具箱提供系統辨識的輸入信號函數idinput    %
%                  生成M序列的過程                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
input=idinput([15 1 100],'prbs',[0 1],[-1 1]);        %u采用15拍M序列
input=input';                                         %輸入input

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      以下利用Model Predictive Control Toolbox    %
%         工具箱建立數學模型并模擬輸入的過程       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d=zeros(1,N);                    
A=[1];
B=[1 0.4 1.5 2.2];
C=[1];
m=idpoly(A,B,C);                         %建立數學模型
d=sim(m,[input'v']);                     %得到輸出值
d=d';
m
order=4;                     %系統階數為4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        LMS         %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w=zeros(1,order);                %初始化w
mu=0.003043;                     %設置最優的mu，由mu.txt運行得到

for i=order:N                    %主要迭代過程，解釋見原理部分
  u=input(i:-1:i-order+1);
  y(i)=w*u';
  e(i)=d(i)-y(i);
  w=w+2*mu*e(i)*u;
  w_delta(i)=norm(w-B);         %w_delta用來記錄每次迭代w接近B的程度
end

figure
subplot(211)
plot(d,'b')
hold on
plot(y,'r')
legend('期望輸出','估計輸出值y')
title('運行時LMS算法輸出值擬和')
xlabel('時間采樣')
ylabel('輸出值')

subplot(212)
semilogy(abs(e));
title('LMS算法輸出誤差值e')
xlabel('時間采樣')
ylabel('誤差值e')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%       RLS       %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w_1=zeros(order,1);

   %初始化w
Lambda=0.99999;       %設置Lambda
Delta=0.00000001;        %設置Delta

c=(1/Delta)*eye(order,order);   %設置c(0)

for i=order:N             %主要迭代過程，解釋見原理部分
u=input(i:-1:i-order+1);

g=c*u'/(Lambda+u*c*u');

c=(c-g*u*c)/Lambda;

y_1(i)=w_1'*u';

e_1(i)=d(i)-y_1(i);

w_1=w_1+g*e_1(i);

w_delta_1(i)=norm(w_1-B');   %w_delta_1用來記錄每次迭代w接近B的程度
end

 

figure;
subplot(211);
plot(d,'b')

hold on;

plot(y_1,'r');

legend('期望輸出','估計輸出值y')
title('運行時RLS算法輸出值擬和') ;

xlabel('時間采樣')

ylabel('輸出值')

 

subplot(212);
semilogy((abs(e_1))) ;

title('RLS算法輸出擬和誤差值e') ;

xlabel('采樣間隔') ;

ylabel('誤差值e') ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        LMS_Advance      %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w_2=zeros(1,order);        %初始化w

for i=order:N              %主要迭代過程，解釋見原理部分
  u=input(i:-1:i-order+1);
  y_2(i)=w_2*u';
  e_2(i)=d(i)-y_2(i);
  if i<33                  %設置最佳拐點，由拐點.txt運行得到
  mu=0.05;                 %設置最佳拐點值，由拐點值.txt運行得到
  else
  mu=0.003043;             %設置最優的mu，由mu.txt運行得到
  end
  w_2=w_2+2*mu*e_2(i)*u;
  w_delta_2(i)=norm(w_2-B);    %w_delta_2用來記錄每次迭代w接近B的程度
end

figure
subplot(211)
plot(d,'b')
hold on
plot(y_2,'r')
legend('期望輸出','估計輸出值y')
title('運行時LMS_Advance算法輸出值擬和')
xlabel('時間采樣')
ylabel('輸出值')

subplot(212)
semilogy(abs(e_2));
title('LMS_Advance算法輸出誤差值e')
xlabel('時間采樣')
ylabel('誤差值e')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%  三種算法比較     %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m_1=idpoly(1,w,1);                         %建立數學模型
d_1=sim(m_1,[input'v']);                   %得到輸出值
figure
subplot(311)
plot(d);                                   %繪圖：LMS算法擬和輸出值
hold on;
plot(d_1,'r');
title('LMS算法擬和輸出值')
legend('期望輸出','擬和輸出值')
xlabel('時間采樣')
ylabel('輸出值')

m_2=idpoly(1,w_1',1);                      %建立數學模型
d_2=sim(m_2,[input'v']);                   %得到輸出值
subplot(312)
plot(d);                                     %繪圖：RLS算法擬和輸出值
hold on;
plot(d_2,'r');
title('RLS算法擬和輸出值')
legend('期望輸出','擬和輸出值')
xlabel('時間采樣')
ylabel('輸出值')

m_3=idpoly(1,w_2,1);                         %建立數學模型
d_3=sim(m_3,[input'v']);                    %得到輸出值
subplot(313)
plot(d);                                     %繪圖：LMS_Advance算法擬和輸出值
hold on;
plot(d_3,'r');
title('LMS_Advance算法擬和輸出值')
legend('期望輸出','擬和輸出值')
xlabel('時間采樣')
ylabel('輸出值')

figure                              %繪圖：擬和輸出誤差比較，并計算各個算方法誤差矩陣的范數
plot(abs(d-d_1'));
norm_error_output_LMS=norm(abs(d-d_1'))
hold on;
plot(abs(d-d_2'),'r');
error_output_RLS=norm(abs(d-d_2'))
hold on;
plot(abs(d-d_3'),'g');
norm_error_output_LMS_Advance=norm(abs(d-d_3'))
hold on;
title('擬和輸出誤差比較')
legend('LMS估計的擬和輸出誤差','RLS估計的擬和輸出誤差','LMS_Advance估計的擬和輸出誤差')
xlabel('時間采樣')
ylabel('輸出誤差')

figure                           %繪圖：擬和輸出誤差比較，并顯示各個算法的權值估計結果
plot(B,'bo')
Actual_weights=B
hold on
plot(w,'r*')
w_LMS=w
hold on
plot(w_1,'g*')
w_RLS=w_1'
hold on
plot(w_2,'y*')
w_LMS_Advance=w_2
legend('實際權值','LMS估計的權值','RLS估計的權值','LMS_Advance估計的權值')
title('三種算法估計權值的比較')

figure                           %繪圖：三種算法權值收斂性的比較，并計算穩定后各個算法權值估                                        計誤差矩陣的范數
semilogy(abs(w_delta));

norm_error_w_LMS=norm(w_delta(800:900))
hold on;
semilogy(abs(w_delta_1),'r');

norm_error_w_RLS=norm(w_delta_1(800:900))
hold on;
semilogy(abs(w_delta_2),'g');

norm_error_w_LMS_Advance=norm(w_delta_2(800:900))
legend('LMS算法估計與真實權值的誤差變化','RLS算法估計與真實權值的誤差變化','LMS_Advance算法估計與真實權值的誤差變化')
title('三種算法權值收斂性的比較') ;