Lekshmi B S
Member level 1
- Joined
- Apr 22, 2013
- Messages
- 34
- Helped
- 2
- Reputation
- 4
- Reaction score
- 2
- Trophy points
- 1,288
- Location
- India
- Activity points
- 1,514
this is the code i formulated for a kalman filter channel estimation. i would like to add a sine input to the same. can anyone help me through this please.
thank you..
thank you..
Code:
clc;
clear all;
close all;
N=20;
F=2.8;
H=0.07;
Q=7;
R=4;
x=zeros(1,N);
z=zeros(1,N);
residual=zeros(1,N);
K=zeros(1,N);
apriori_state=zeros(1,N);
apriori_error=ones(1,N);
aposteriori_state=zeros(1,N);
aposteriori_error=ones(1,N);
x_0=1;
aposteriori_state_0=1.5;
aposteriori_error_0=1;
randn('state',0);
w = wgn(1,N,5)*sqrt(Q);
v = wgn(1,N,5)*sqrt(R);
x(1)=F*x_0+w(1);
z(1)=H*x(1)+v(1);
apriori_state(1)=F*aposteriori_state_0;
apriori_error(1)=F*F'*aposteriori_error_0+Q;
residual(1)=z(1)-H*apriori_state(1);
K(1)=H*apriori_error(1)/(H*H'*apriori_error(1)+R);
aposteriori_state(1)=apriori_state(1)+K(1)*residual(1);
aposteriori_error(1)=apriori_error(1)*(1-(H*K(1)));
error1(1)=apriori_state(1)-x(1);
error2(1)=aposteriori_state_0-x_0;
c=1;
for i=2:N
x(i)=F*x(i-1)+w(i);
z(i)=H*x(i)+v(i);
apriori_state(i)=F*aposteriori_state(i-1);
apriori_error(i)=F*F'*aposteriori_error(i-1)+Q;
residual(i)=z(i)-H*apriori_state(i);
K(i)=H*apriori_error(i)/(H*H'*apriori_error(i)+R);
aposteriori_state(i)=apriori_state(i)+K(i)*residual(i);
aposteriori_error(i)=apriori_error(i)*(1-(H*K(i)));
error1(i)=apriori_state(i)-x(i)
error2(i)=aposteriori_state(i-1)-x(i-1)
c=c+1;
if (abs(error1(i))<=1)
break
end
end
abs(error1)
figure(1);
t=1:N;
plot(t,K);
xlabel('Observation number');
ylabel('Kalman gain');
title('Kalman gain');
figure(2);
t=1:N;
subplot(2,1,1);
plot(t,w);
xlabel('Observation number');
ylabel('w');
title('Process noise');
subplot(2,1,2);
plot(t,v);
xlabel('Observation number');
ylabel('v');
title('Observation noise');
figure(3);
t=1:c;
subplot(2,1,1);
plot(t,error1);
xlabel('Observation number');
ylabel('Error1');
title('Error between actual value and predicted value');
subplot(2,1,2);
plot(t,error2);
xlabel('Observation number');
ylabel('Error2');
title('Error between actual value and updated value');
figure(4);
m=1:c;
subplot(2,1,1);
plot(error1,K(m));
xlabel('Error1');
ylabel('Kalman gain');
title('Error1 Vs Kalman gain');
subplot(2,1,2);
plot(error2,K(m));
xlabel('Error2');
ylabel('Kalman gain');
title('Error2 Vs Kalman gain');