## auto correlation in scilab

February 8, 2011 Coded in Scilab
//Autocorrelation of a given Input Sequence
//Finding out the period of the signal using autocorrelation technique
clear;
clc;
close;
x = input('Enter the given discrete time sequence');
L = length(x);
h = zeros(1,L);
for i = 1:L
h(L-i+1) = x(i);
end
N = 2*L-1;
Rxx = zeros(1,N);
for i = L+1:N
h(i) = 0;
end
for i = L+1:N
x(i) = 0;
end
for n = 1:N
for k = 1:N
if(n >= k)
Rxx(n) = Rxx(n)+x(n-k+1)*h(k);
end
end
end
disp('Auto Correlation Result is')
Rxx
disp('Center Value is the Maximum of autocorrelation result')
[m,n] = max(Rxx)
disp('Period of the given signal using Auto Correlation Sequence')
n

## sinc function in scilab

February 8, 2011 Coded in Scilab
function [y]=sinc_new(x)
i=find(x==0);
x(i)= 1;      // don't need this is /0 warning is off
y = sin(%pi*x)./(%pi*x);
y(i) = 1;
endfunction

## Discrete Cosine Tranform Matrix generation

February 6, 2011 Coded in Scilab
function[DCT] = dct_mtx(n)
[cc,rr] = meshgrid(0:n-1);
//disp(cc)
//disp(rr)
DCT = sqrt(2 / n) * cos(%pi * (2*cc + 1) .* rr / (2 * n));
DCT(1,:) = DCT(1,:) / sqrt(2);
endfunction

## Discrete Fourier Tranform Matrix generation

February 6, 2011 Coded in Scilab
function [D] = dft_mtx(n)
f = 2*%pi/n;                 // Angular increment.
w = (0:f:2*%pi-f/2).' *%i;   //Column.
//disp(w)
x = 0:n-1;                  // Row.
D = exp(-w*x);              // Exponent of outer product.
for i = 1:n
for j = 1:n
if((abs(real(D(i,j)))<0.0001)&(abs(imag(D(i,j)))<0.0001))
D(i,j)=0;
elseif(abs(real(D(i,j)))<0.0001)
D(i,j)= 0+%i*imag(D(i,j));
elseif(abs(imag(D(i,j)))<0.0001)
D(i,j)= real(D(i,j))+0;
end
end
end
endfunction

## 2D IFFT

February 2, 2011 Coded in Scilab
function [a] =ifft2d(a2)
//a2 = 2D-DFT of any real or complex 2D matrix
//a = 2D-IDFT of a2
m=size(a2,1)
n=size(a2,2)
//Inverse Fourier transform along the rows
for i=1:n
a1(:,i)=exp(2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a2(:,i)
end
//Inverse fourier transform along the columns
for j=1:m
atemp=exp(2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).'
a(j,:)=atemp.'
end
a = a/(m*n)
a = real(a)
endfunction

## 2D FFT

February 2, 2011 Coded in Scilab
function [a2] = fft2d(a)
//a = any real or complex 2D matrix
//a2 = 2D-DFT of 2D matrix  'a'
m=size(a,1)
n=size(a,2)
// fourier transform along the rows
for i=1:n
a1(:,i)=exp(-2*%i*%pi*(0:m-1)'.*.(0:m-1)/m)*a(:,i)
end
// fourier transform along the columns
for j=1:m
a2temp=exp(-2*%i*%pi*(0:n-1)'.*.(0:n-1)/n)*(a1(j,:)).'
a2(j,:)=a2temp.'
end
for i = 1:m
for j = 1:n
if((abs(real(a2(i,j)))<0.0001)&(abs(imag(a2(i,j)))<0.0001))
a2(i,j)=0;
elseif(abs(real(a2(i,j)))<0.0001)
a2(i,j)= 0+%i*imag(a2(i,j));
elseif(abs(imag(a2(i,j)))<0.0001)
a2(i,j)= real(a2(i,j))+0;
end
end
end

## Basic Kalman filter

January 20, 201114 comments Coded in Scilab
// Clear initialisation
clf()

// Time vector
t =[0:0.5:100];
len_t = length(t);
dt = 100 / len_t;

// Real values x: position v: velocity a: acceleration
x_vrai = 0 * t;
v_vrai = 0 * t;
a_vrai = 0 * t;

// Meseared values.
x_mes = 0 * t;
v_mes = 0 * t;
a_mes = 0 * t;

// Initialisation
for i = 0.02 * len_t: 0.3 * len_t,
a_vrai(i) = 1;
end

for i = 0.5 * len_t: 0.7 * len_t,
a_vrai(i) = 2;
end

// State of kalman filter
etat_vrai = [ x_vrai(1);  v_vrai(1); a_vrai(1) ];

A = [ 1 dt dt*dt/2; 0 1 dt; 0 0 1];//transition matrice
for i = 2:length(t),
etat_vrai = [ x_vrai(i-1);  v_vrai(i-1); a_vrai(i) ];
etat_vrai = A * etat_vrai;
x_vrai(i) = etat_vrai(1);
v_vrai(i) = etat_vrai(2);
end

max_brt_x = 200;
max_brt_v = 10;
max_brt_a = 0.3;

// Random noise
x_brt = grand(1, len_t, 'unf', -max_brt_x, max_brt_x);
v_brt = grand(1, len_t, 'unf', -max_brt_v, max_brt_v);
a_brt = grand(1, len_t, 'unf', -max_brt_a, max_brt_a);

// Compute meas
x_mes = x_vrai + x_brt;
v_mes = v_vrai + v_brt;
a_mes = a_vrai + a_brt;

// START
x_est = 0 * t;
v_est = 0 * t;
a_est = 0 * t;

etat_mes = [ 0; 0; 0 ];
etat_est = [ 0; 0; 0 ];
inn = [ 0; 0; 0 ];
// Matrix of covariance of the bruit.
// independant noise -> null except for diagonal,
// and on the diagonal == sigma^2.
// VAR(aX+b) = a^2 * VAR(X);
// VAR(X) = E[X*X] - E[X]*E[X]
R = [ max_brt_x*max_brt_x 0 0; 0 max_brt_v*max_brt_v 0; 0 0 max_brt_a*max_brt_a ];

Q = [1 0 0; 0 1 0; 0 0 1];
P = Q;
for i = 2:length(t),
// Init : measured state
etat_mes = [ x_mes(i);  v_mes(i); a_mes(i) ];

// Innovation: measured state -  estimated state
inn = etat_mes - etat_est;

// Covariance of the innovation
S = P + R;

// Gain
K = A * inv(S);

// Estimated state
etat_est = A * etat_est + K * inn;

// Covariance of prediction error
// C = identity
P = A * P * A' + Q - A * P * inv(S) * P * A';

// End: estimated state
x_est(i) = etat_est(1);
v_est(i) = etat_est(2);
a_est(i) = etat_est(3);

end

// Blue : real
// Gree : noised
// Red: estimated
subplot(311)
plot(t, a_vrai, t, a_mes, t, a_est);
subplot(312)
plot(t, v_vrai, t, v_mes, t, v_est);
subplot(313)
plot(t, x_vrai, t, x_mes, t, x_est);