[SOLVED] Kalman filter for channel estimation

[Lekshmi B S]

Can anyone please provide me with the matlab code for channel estimation using a kalman filter (in presence of noise)...
It is urgent...
Thank u...

 

[Lekshmi B S]

i have formulated the code for the above..but i would like to add gaussian noise to the system.,..can anyone explain how a statistical gaussian noise is generated in matlab.....i think the in-built functions such as randn, wgn, awgn etc are used for generating random noise...am not much clear about this..please help me to solve this issue....

Thank you...

 

[phongphanp]

I think the type of test data set is important , also.

 

[Lekshmi B S]

Isn't dis code enough?
this is the matlab code i formulated using the equations. the noise added is gaussian. but the problem is that i am getting a constant error instead of error getting reduced with the number of iterations.
can u find the mistakes in this.. Do help...
Thank you..

N=20; % number of iterations
F=0.4; % state transition model
H=0.07; % observation model
Q=0.4; % covariance of process noise
R=0.8; % covariance of observation noise

x=zeros(1,N); % state estimate vector
z=zeros(1,N); % output vector
residual=zeros(1,N); % measurement residual
K=zeros(1,N); % Kalman gain
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;

w(1)=(sqrt(1/(((2*pi)^N)*det(Q))))*(exp((-x_0'*inv(Q)*x_0)/2)); % process noise
v(1)=(sqrt(1/(((2*pi)^N)*det(R))))*(exp((-x_0'*inv(R)*x_0)/2)); % observation noise

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)*residua l(1);
aposteriori_error(1)=apriori_error(1)*(1-(H*K(1)));

for i=2:N
w(i)=(sqrt(1/(((2*pi)^N)*det(Q))))*(exp((-x(i-1)'*inv(Q)*x(i-1))/2));
v(i)=(sqrt(1/(((2*pi)^N)*det(R))))*(exp((-x(i-1)'*inv(R)*x(i-1))/2));
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)*residua l(i);
aposteriori_error(i)=apriori_error(i)*(1-(H*K(i)));
end

error1=x-apriori_state;
error2=x-aposteriori_state;

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,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(3);
subplot(2,1,1);
plot(error1,K);
xlabel('Error1');
ylabel('Kalman gain');
title('Error1 Vs Kalman gain');
subplot(2,1,2);
plot(error2,K);
xlabel('Error2');
ylabel('Kalman gain');
title('Error2 Vs Kalman gain');

- - - Updated - - -

i used the equation for gaussian distribution to generate the noise in this case..

 

[Lekshmi B S]

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..

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');