DSPRelated.com
   Code Snippets

auto correlation in scilab

Senthilkumar ● 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

Senthilkumar ● 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

Senthilkumar ● 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

Senthilkumar ● 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

Senthilkumar ● 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

Senthilkumar ● 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, 2011●14 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);