function [sig_out]=adaptive_filter(main_sig,ref_sig,slector,M)
% 自适应滤波程序
% main_sig 主通道信号
% ref_sig 参考通道信号
% slector 滤波算法选择，可选LMS、NLMS、RLS
% M 滤波器阶数
% sig_out 滤波后的输出信号
% by milttuso 2011.10.17
% 当信号的功率发生改变时，此时各算法的参数可能不是最适合的了，需要试验一下
u=ref_sig.'; %滤波器输入信号
d=main_sig.'; %滤波器期望信号
% slector='LMS';
switch slector
%%%%%%%%%%% ======== LMS算法 =======%%%%%%%%%%%%
%%%% 存在误差放大问题。。。。关键是要取一个合适的步长
case 'LMS'
mu=0.005e-4; %步长因子
w_lms=zeros(M,1); %权向量初始化
N=length(u);
e_lms=zeros(1,N);
sig_out=zeros(1,N); %输出序列初始化
for n=1:N
if n < M
ui = [ref_sig(n:-1:1), zeros(1, M - n)].';
else
ui = ref_sig(n:-1:n - M + 1).';
end
sig_out(n)=d(n)-w_lms'*ui;
w_lms=w_lms+mu*conj(sig_out(n))*ui;
end
w_opt=w_lms;
e_nlms=sig_out;
% figure
% plot(abs(e_nlms))
% xlabel('迭代次数')
% ylabel('误差曲线')
for n=1:N
if n < M
ui = [ref_sig(n:-1:1), zeros(1, M - n)].';
else
ui = ref_sig(n:-1:n - M + 1).';
end
sig_out(n)=d(n)-w_opt'*ui;
end
%%%%%%%% ========NLMS算法=========%%%%%%%%%%%%%%%
case 'NLMS'
mu=0.16; %步长因子
omega=0.00004;
w_nlms=zeros(M,1); %权向量初始化
N=length(u);
e_nlms=zeros(1,N);
sig_out=zeros(1,N); %输出序列初始化
for n=1:N
if n < M
ui = [ref_sig(n:-1:1), zeros(1, M - n)].';
else
ui = ref_sig(n:-1:n - M + 1).';
end
sig_out(n)=d(n)-w_nlms'*ui;
w_nlms=w_nlms+mu*conj(sig_out(n))*ui/(omega+ui'*ui);
end
w_opt=w_nlms;
e_nlms=sig_out;
% figure
% plot(abs(e_nlms))
% xlabel('迭代次数')
% ylabel('误差曲线')
for n=1:N
if n < M
ui = [ref_sig(n:-1:1), zeros(1, M - n)].';
else
ui = ref_sig(n:-1:n - M + 1).';
end
sig_out(n)=d(n)-w_opt'*ui;
end
%%%%%%%%%%%%==========RLS算法===========%%%%%%%%%%%%%%%%%
case 'RLS'
delta=0.004;
lamda=1;
w_rls=zeros(M,1); %权向量初始化
N=length(u);
sig_out=zeros(1,N); %输出序列初始化
P=eye(M)/delta ; %初始化P
e_rls=zeros(N,1);
for n=1:N
if n < M
ui = [ref_sig(n:-1:1), zeros(1, M - n)].';
else
ui = ref_sig(n:-1:n - M + 1).';
end
kn=1/lamda*P*ui/(1+1/lamda*ui'*P*ui); %k(n)
e_rls(n,:)=d(n)-w_rls'*ui;
w_rls=w_rls+kn*e_rls(n,:)';
P=P/lamda-kn*ui'*P/lamda;
end
w_opt=w_rls;
% figure
% plot(abs(e_rls))
% xlabel('迭代次数')
% ylabel('误差曲线')
for n=1:N
if n < M
ui = [ref_sig(n:-1:1), zeros(1, M - n)].';
else
ui = ref_sig(n:-1:n - M + 1).';
end
sig_out(n)=d(n)-w_opt'*ui;
end
end
