马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
clear
clc
M=3;L=12;yz=0.01;yu=0.90;dB=25;
H0=[0.62,-0.25,0.64,0.82,0,0,0,0,-0.92,0,-0.92,0]';
% syms k p i;
power1=0.5; N=5000;
y1=zeros(1,N-99);
u=rand(1,N);
a1=sqrt(12*power1);
u1=a1*u;
SNR=dB;
power2=power1/(10.^(SNR/20));
noise1=sqrt(power2)*randn(1,N);
u=u1+noise1;
for k=100:N
y(k)=[0.62 -0.25 0.64 0.82 1.84]*[u(k) u(k-1) u(k-2) u(k).^2 u(k-1).*u(k-2)]';
X(:,k)=[u(k),u(k-1),u(k-2),u(k).^2,u(k).*u(k-1),u(k).*u(k-2),u(k-1).*u(k),u(k-1).^2,u(k-1).*u(k-2),u(k-2).*u(k),u(k-2).*u(k-1),u(k-2).^2]';
end
power_y=var(y);
power3=power_y/(10.^(SNR/20));
noise2=sqrt(power3).*randn(1,N);
y=y+noise2;
P0=eye(L,L);
P100=[X(:,100) P0(:,2:L)]';R100=inv(P100'*P100);Y=[y(100) zeros(1,L-1)];
for k=100
Pk=P100;Yk=Y;
end
for k=100:N-1
count=0;
for j=1:L
if abs(X(j,k+1)-X(j,k))<yz
count=count+1;j=j+1;
else j=j+1;
end
end
if count/L<yu
Xp=X(:,k);Xc=Pk(1:L-1,:);Xl=Pk(L,:); Yc=Yk(:,1:L-1);
Rk=inv(Pk'*Pk);
Sk=Rk-Rk*Xp*Xp'*Rk/(1+Xp'*Rk*Xp);
k=k+1;
Xp=X(:,k);Yp=y(k);Pk=[Xp Xc']';Yk=[Yp Yc];
Rk=Sk+Sk*Xl'*Xl*Sk/(1-Xl*Sk*Xl');
Hk=[Sk+Sk*Xl'*Xl*Sk/(1-Xl*Sk*Xl')]*Pk'*Yk';
else
k=k+1;
end
H(:,k)=Hk;
e(k)=20*log10((norm(H0-H(:,k)))^2);
end
plot(e(100:600),'r');grid on
plot(H(1,100:600),'c');grid on
plot(H(2,100:600));grid on
plot(H(3,100:600));grid on |