clear all;clc;
N=2000; %2000个样本点
M=150; %滤波器阶数
n=1:N;
s=5*sin(0.05*n)'; %期望信号
v=normrnd(0,5^0.5,1,N)'; %高斯噪声
x=s+v; %观测信号
y=zeros(1,N)'; %用于存储输出结果
e=y; %用于存储误差
Ee2=y; %用于存储均方误差
rxx=xcorr(x,'unbiased');
Rxx=toeplitz(rxx(N:N+N-1));
W=zeros(M,N); %加权矢量初始值，并按列存储
U=2/sum(diag(Rxx,0));
u=U/2; %收敛因子应该满足条件u<2/tr(Rxx)
for k=M:N
xk=x(k:-1:k-M+1); %第k次迭代时，滤波器的M个输入数据
y(k)=W(:,k)'*xk; %滤波器输出
e(k)=s(k)-y(k); %误差
W(:,k+1)=W(:,k)+u*e(k)*xk; %更新权值矢量,注意W按列存储
Ee2(k)=e(k)*e(k);
end
yn=zeros(1,N)';
for j=M:N
xj=x(j:-1:j-M+1);
yn(j)=W(:,N)'*xj; %最佳权矢量存储在W(:,N)中
end
figure(1)
subplot(221)
plot(s);
title('期望信号');
subplot(222)
plot(v)
title('干扰信号');
subplot(223)
plot(x)
title('观测信号');
%因为M之前的观测信号没有进行滤波处理,所以图像只画了前M之后的数据
subplot(224)
plot(yn);
title('滤波后信号');
%画加权矢量随递推次数变化的三维图像
figure(2)
p=M:N;
q=1:M;
[p,q]=meshgrid(p,q);
mesh(W(:,M:N));
x1=xlabel('递推次数'); %x轴标题
x2=ylabel('加权矢量下标'); %y轴标题
x3=zlabel('加权矢量W'); %z轴标题
set(x1,'rotation',15); %x轴标题旋转
set(x2,'rotation',-25); %y轴标题旋转
title('加权矢量变化图');
figure(3)
plot(Ee2);
xlabel('递推次数');
ylabel('均方误差');
title('均方误差');
评论0