%Generation of Time-varying channel co-efficient (real)
mu=1;
a=0.2;
sigma2=1;
sigmaw2=0.1;
sigmav2=1;
h=sqrt(sigmaw2+(aˆ2)*sigma2)*randn(1)+a*mu;
LEN=100;
for i=2:1:LEN
h(i)=a*h(i-1)+randn(1)*sqrt(sigmaw2);
end
%a,sigmaw2 and sigmav2 are known
%We would like to track the time varying channel co-efficients based on the
%output corresponding to the pilot input signal and the corresponding output.
x=(round(rand(1,LEN))*2-1)*2;
y=[];
for i=1:1:LEN
y=[y h(i)*x(i)+sqrt(sigmav2)*randn(1)];
end
%Initialize
p_n(1)=1; hcapn(1)=1.5;
for k=2:1:LEN
p_n_n_1(k)=(a^2)*p_n(k-1)+sigmaw2;
kalman(k)=x(k)*p_n_n_1(k)/((x(k)^2)*p_n_n_1(k)+sigmav2);
hcapn_n_1(k)=a*hcapn(k-1);
residual(k)=y(k)-x(k)*hcapn_n_1(k);
hcapn(k)=kalman(k)*residual(k)+hcapn_n_1(k);
p_n(k)=(1-x(k)*kalman(k))*p_n_n_1(k);
end
figure
subplot(3,1,1)
plot(x(2:1:end))
title('Input sequence')
subplot(3,1,2)
plot(y(2:1:end))
title('Corresponding output sequence')
subplot(3,1,3)
plot(2:1:LEN,h(2:1:end)/sqrt(sum(h(2:1:end).^2)))
hold on
plot(2:1:LEN,hcapn(2:1:end)/sqrt(sum(hcapn(2:1:end).^2)),'r')
title('Actual and the estimated time varying filter co-efficent')
figure
subplot(5,1,1)
plot(2:1:LEN,p_n_n_1(2:1:end))
title('Variance of the estimation of h(n) given y(1) y(1)...y(n-1)');
subplot(5,1,2)
plot(2:1:LEN,kalman(2:1:end))
title('Kalman gain used to update the estimation of h(n) given y(1)...y(n)')
subplot(5,1,3)
plot(2:1:LEN,hcapn_n_1(2:1:end))
title('Estimation of h(n) given y(1)...y(n-1)');
subplot(5,1,4)
plot(2:1:LEN,residual(2:1:end))
title('Residual to update the estimation of h(n) given y(1)...y(n)')
subplot(5,1,4)
plot(2:1:LEN,hcapn(2:1:end))
title('Estimation of h(n) given y(1)...y(n-1) y(n)');
subplot(5,1,5)
plot(2:1:LEN,p_n(2:1:end))
title('Variance of the estimation of h(n) given y(1)...y(n)');
