【討論】自適應濾波的MATLAB實現(xiàn)
1.LMS算法的仿真程序：
%lms 算法
clear all
close all
hold off%系統(tǒng)信道權數(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 = n * std(y)/(10*std(n));%噪聲信號
d = y + n;%期望輸出信號
totallength=size(d,1);%步長
N=60 ; %60節(jié)點作為訓練序列
%算法的開始
w = zeros ( sysorder , 1 ) ;%初始化
for n = sysorder : N 
u = inp(n:-1:n-sysorder+1) ;% u的矩陣
y(n)= w' * u;%系統(tǒng)輸出
e(n) = d(n) - y(n) ;%誤差
if n < 20
mu=0.32;
else
mu=0.15;
end
w = w + mu * u * e(n) ;%迭代方程
end 
%檢驗結果
for n = N+1 : totallength
u = inp(n:-1:n-sysorder+1) ;
y(n) = w' * u ;
e(n) = d(n) - y(n) ;%誤差
end 
hold on
plot(d)
plot(y,'r');
title('系統(tǒng)輸出') ;
xlabel('樣本')
ylabel('實際輸出')
figure
semilogy((abs(e))) ;% e的絕對值坐標
title('誤差曲線') ;
xlabel('樣本')
ylabel('誤差矢量')
figure%作圖
plot(h, 'k+')
hold on
plot(w, 'r*')
legend('實際權矢量','估計權矢量')
title('比較實際和估計權矢量') ;
axis([0 6 0.05 0.35])
2.NLMS算法的仿真程序：
%lms 算法
clear all
close all
hold off%系統(tǒng)信道權數(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 = n * std(y)/(10*std(n));%噪聲信號
d = y + n;%期望輸出信號
totallength=size(d,1);%步長
N=60 ; %60節(jié)點作為訓練序列
%算法的開始
w = zeros ( sysorder , 1 ) ;%初始化
for n = sysorder : N 
u = inp(n:-1:n-sysorder+1) ;% u的矩陣
y(n)= w' * u;%系統(tǒng)輸出
r(n)=u'*u;%自相關矩陣
e(n) = d(n) - y(n) ;%誤差
fai=.0001;%修正參數(shù)，為防止u'*u過小導致步長值太大而設置的
if n < 20
mu=0.32;
else
mu=0.15;
end
w = w + mu * u * e(n)/(r(n)+fai) ;;%迭代方程
end 
%檢驗結果
for n = N+1 : totallength
u = inp(n:-1:n-sysorder+1) ;
y(n) = w' * u ;
e(n) = d(n) - y(n) ;%誤差
end 
hold on
plot(d)
plot(y,'r');
title('系統(tǒng)輸出') ;
xlabel('樣本')
ylabel('實際輸出')
figure
semilogy((abs(e))) ;% e的絕對值坐標
title('誤差曲線') ;
xlabel('樣本')
ylabel('誤差矢量')
figure%作圖
plot(h, 'k+')
hold on
plot(w, 'r*')
legend('實際權矢量','估計權矢量')
title('比較實際和估計權矢量') ;
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 training
Order = 32 ; % 自適應濾波權數(shù)
Lambda = 0.98 ; % 遺忘因子
Delta = 0.001 ; % 相關矩陣R的初始化
x = randn(NoOfData, 1) ;%高斯隨機系列
h = rand(Order, 1) ; % 系統(tǒng)隨機抽樣
d = filter(h, 1, x) ; % 期望輸出
% RLS算法的初始化
P = Delta * eye ( Order, Order ) ;%相關矩陣
w = zeros ( Order, 1 ) ;%濾波系數(shù)矢量的初始化
% RLS Adaptation
for n = Order : NoOfData ; 
u = x(n:-1:n-Order+1) ;%延時函數(shù)
pi_ = u' * P ;%互相關函數(shù)
k = Lambda + pi_ * u ;
K = pi_'/k;%增益矢量
e(n) = d(n) - w' * u ;%誤差函數(shù)
w = w + K * e(n) ;%遞歸公式
PPrime = K * pi_ ;
P = ( P - PPrime ) / Lambda ;%誤差相關矩陣
w_err(n) = norm(h - w) ;%真實估計誤差
end ;
% 作圖表示結果
figure ;
plot(20*log10(abs(e))) ;%| e |的誤差曲線
title('學習曲線') ;
xlabel('迭代次數(shù)') ;
ylabel('輸出誤差估計') ;
figure ;
semilogy(w_err) ;%作實際估計誤差圖
title('矢量估計誤差') ;
xlabel('迭代次數(shù)') ;
ylabel('誤差權矢量') ;
4.自適應均衡器的仿真程序：
% 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));
h(2)=0.5*(1+cos(2*pi*(2-2)/W));
h(3)=0.5*(1+cos(2*pi*(3-2)/W));
% 學習曲線
hc=[0 h(1) h(2) h(3)]';
n0=7;
t=(1:N)';
for i=1:Nmc
y=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 算法
% Initialization
lambda=0.98;
P=(10^-3)*eye(M,M);
c=zeros(M,1);
g=c;
glambda=g;
% 迭代范圍
for n=M:M+N-1
xn=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));
end
end
eplot=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.自適應陷波器的仿真程序：
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(2*pi*t/10);
%初始化
w1=0.1;
w2=0.1;
e=zeros(1,N);
y=0;
u=0.05;%迭代步長
for i=1:N
y=w1*x1(i)+w2*x2(i);
e(i)=x(i)-y;%誤差信號
w1=w1+u*e(i)*x1(i);%迭代方程
w2=w2+u*e(i)*x2(i);%迭代方程
end
subplot(2,2,3); %作第三子圖
plot(t,e);
subplot(2,2,4); %作第四子圖
plot(t,s-e);