rls.m

振动抑制的matlab程序源码 RLS算法 matlab程序

function RLS

% IIR_rls : IIR RLS algorithm demo
clear all
close all
clc

inporder=2;		%X的阶数
outorder=2;		%Y的阶数
sysorder = inporder + outorder ;	%总阶数，w用				

N = 8000;
T = 1:1:N;
f = 5;                              %信号频率
fc = 1000;                          %采样频率，代表每秒采样多少个点 

noise1 = 0.3*randn(1,N);            %噪声信号
% noise1=sin(2*pi*f/fc*T);          %正弦信号
s=sin(2*pi*f/fc*T);                 %输入的有效信号(期望的信号)

inp= noise1 + s;					%实际输入信号
d0 = inp';							%转置
noise = noise1';
d = d0;

%主通道传递函数，可选，不选时可认为是1
% b = [0.05 -0.4];
% a = [1 -1.131 0.25];
% Gz = filt(b,a,- 1); 
% d = lsim(Gz,inp); 		%输出

lamda = 0.999 ;

%initial P matrix
delta = 1e2 ;
P = delta * eye (sysorder ) ;
w = zeros ( sysorder , 1 ) ;		%权衡系数
y = zeros(N,1);						

for n = max(inporder,outorder) + 1 : N
%  	u=d0(n:- 1:n- inporder+1) ;         %整个混合信号作为输入
   u=0.8*noise(n : -1 : n- inporder+1) ;     %仅将噪声作为输入
	outp= y(n-1 : - 1 : n-outorder) ;
    
	u=[u ; outp];   %x(n),x(n-1), y(n-1),y(n-2)
	y(n)=w' * u;
	e(n) = d(n) - y(n) ;
	
	phi = u' * P ;  % row * cof, the diag of P is equal to a1 a2, b0, b1, so phi is [b0*x(n) b1*x(n-1) a1*y(n-1) a2*y(n-2)]
	k = phi' /(lamda + phi * u );   %phi * u = [b0*x(n)^2 b1*x(n-1)^2 a1*y(n-1)^2 a2*y(n-2)^2]
	
	%更新w和p
	w = w + k * e(n) ;
	P = ( P - k * phi ) / lamda ;
 	Recordedw(1:sysorder,n)=w;
end

%b = [0.05 -0.4];		a0,a1
%a = [1 -1.131 0.25];	1, b1,b2
w 	%显示w，可以看到辨识出来的数值跟上述a，b非常接近

hold on 
plot(d0)             %代表如果不采用自适应滤波时的振动情况
plot(0.5*e,'r');    %代表采用了自适应滤波后的振动情况
title('System output') ; 
xlabel('Samples') 
ylabel('True and estimated output') 

figure
plot(Recordedw(1:sysorder,sysorder:N)' );
title(' 估计权值收敛' ) ;
xlabel(' 抽样' );
ylabel(' 权值' );
axis([1 N- sysorder min(min(Recordedw(1:sysorder,sysorder:N)' )) max(max(Recordedw(1:sysorder,sysorder:N)' )) ]);
hold off