DSPRelated.com
   Code Snippets

Decimation or Downsampling in scilab

Senthilkumar ● March 21, 2011●2 comments ● Coded in Scilab 
//Multirate Signal Processing in scilab
//Downsampling a sinusoidal signal by a factor of 2
clear;
clc;
n = 0:%pi/200:2*%pi;
x = sin(%pi*n);   //original signal
downsampling_x = x(1:2:length(x)); //downsampled by a factor of 2
subplot(2,1,1)
plot(1:length(x),x);
xtitle('original singal')
subplot(2,1,2)
plot(1:length(downsampling_x),downsampling_x);
xtitle('Downsampled Signal by a factor of 2');

Discrete Plot in scilab

Senthilkumar ● March 21, 2011●1 comment ● Coded in Scilab 
//This Program Illustrates the discrete plot in scilab
//using plot2d3 function
clear;
clc;
close;
a =1.5; 
n =1:10;  
x = (a)^n;
a=gca();
a.thickness = 2;
plot2d3('gnn',n,x) 
xtitle('Graphical Representation of Exponentially Increasing Signal','n','x[n]');

Continuous Time Fourier Transform of a Square Waveform

Senthilkumar ● March 21, 2011 ● Coded in Scilab 
// Continuous Time Fourier Transform and Frequency Response of a Square Waveform
// x(t)= A, from -T1 to T1
clear;
clc;
close;
// CTS Signal
A =1;    //Amplitude
Dt = 0.005;
T1 = 4;  //Time in seconds
t = -T1:Dt:T1;
for i = 1:length(t)
  xt(i) = A;
end
//
// Continuous-time Fourier Transform
Wmax = 2*%pi*1;        //Analog Frequency = 1Hz
K = 4;
k = 0:(K/1000):K;
W = k*Wmax/K;
xt = xt';
XW = xt* exp(-sqrt(-1)*t'*W) * Dt;
XW_Mag = real(XW);
W = [-mtlb_fliplr(W), W(2:1001)]; // Omega from -Wmax to Wmax
XW_Mag = [mtlb_fliplr(XW_Mag), XW_Mag(2:1001)];
//
subplot(2,1,1);
a = gca();
a.data_bounds=[-4,0;4,2];
a.y_location ="origin";
plot2d2(t,xt);
xlabel('t in msec.');
title('Continuous Time Signal x(t)')
subplot(2,1,2);
a = gca();
a.y_location ="origin";
a.x_location ="origin";
plot(W,XW_Mag);
xlabel('Frequency in Radians/Seconds');
title('Continuous-time Fourier Transform  X(jW)')

Continuous Time Fourier Transform of a Exponential Signal

Senthilkumar ● March 21, 2011 ● Coded in Scilab 
//Continuous Time Fourier Transform of a
//Continuous Time Exponential Signal x(t)= exp(-A*t)u(t), A>0
//And Plotting its Magnitude response and phase response
clear;
clc;
close;
// Analog Signal
A =1;    //Amplitude
Dt = 0.005;
t = 0:Dt:10;
xt = exp(-A*t);
//
// Continuous-time Fourier Transform
Wmax = 2*%pi*1;        //Analog Frequency = 1Hz
K = 4;
k = 0:(K/1000):K;
W = k*Wmax/K;
XW = xt* exp(-sqrt(-1)*t'*W) * Dt;
XW_Mag = abs(XW);
W = [-mtlb_fliplr(W), W(2:1001)]; // Omega from -Wmax to Wmax
XW_Mag = [mtlb_fliplr(XW_Mag), XW_Mag(2:1001)];
[XW_Phase,db] = phasemag(XW);
XW_Phase = [-mtlb_fliplr(XW_Phase),XW_Phase(2:1001)];
//Plotting Continuous Time Signal
figure
a = gca();
a.y_location = "origin";
plot(t,xt);
xlabel('t in sec.');
ylabel('x(t)')
title('Continuous Time Signal')
figure
//Plotting Magnitude Response of CTS
subplot(2,1,1);
a = gca();
a.y_location = "origin";
plot(W,XW_Mag);
xlabel('Frequency in Radians/Seconds---> W');
ylabel('abs(X(jW))')
title('Magnitude Response (CTFT)')
//Plotting Phase Reponse of CTS
subplot(2,1,2);
a = gca();
a.y_location = "origin";
a.x_location = "origin";
plot(W,XW_Phase*%pi/180);
xlabel('                         Frequency in Radians/Seconds---> W');
ylabel('                                                <X(jW)')
title('Phase Response(CTFT) in Radians')

CTFS coefficients of a periodic square waveform

Senthilkumar ● March 21, 2011 ● Coded in Scilab 
//CTFS coefficients of a periodic square waveform 
//x(t) = 1, |t|<T1, and 0, T1<|t|<T/2
clear;
close;
clc;
T =4;
T1 = T/4;
t = -T1:T1/100:T1;
Wo = 2*%pi/T;
xt =ones(1,length(t)); //Square Wave Generation
//Finding the Fourier Series Coefficients
for k =0:5
  C(k+1,:) = exp(-sqrt(-1)*Wo*t.*k);
  a(k+1) = xt*C(k+1,:)'/length(t);
  if(abs(a(k+1))<=0.1) 
    a(k+1)=0;
  end
end
a =a';
a_conj = real(a(:))-sqrt(-1)*imag(a(:));
ak = [a_conj($:-1:1)',a(2:$)];
k = 0:5;
k = [-k($:-1:1),k(2:$)];
Spectrum_ak = (1/2)*real(ak);
//Plotting the Continuous Time Signal
figure
a = gca();
a.y_location = "origin";
a.x_location = "origin";
a.data_bounds=[-2,0;2,2];
plot2d2(t,xt,5)
poly1 = a.children(1).children(1);
poly1.thickness = 3; 
title('x(t)')
xlabel('                                                       t')
//Plotting the Spectral coefficients of CT Signal
figure
a = gca();
a.y_location = "origin";
a.x_location = "origin";
plot2d3('gnn',k,Spectrum_ak,5)
poly1 = a.children(1).children(1);
poly1.thickness = 3; 
title('abs(ak)')
xlabel('                                                       k')
//Result
//Spectrum_ak  =
//0.0633127  0. -0.1055559 0. 0.3167197 0.5 0.3167197 0. -0.1055559    0.    0.0633127

Continuous Time Fourier Series of sine signal

Senthilkumar ● March 21, 2011 ● Coded in Scilab 
//Continuous Time Fourier Series Spectral Coefficients of 
//a periodic sine signal x(t) = sin(Wot)
clear;
close;
clc;
t = 0:0.01:1;
T = 1;
Wo = 2*%pi/T;
xt = sin(Wo*t);
for k =0:5
  C(k+1,:) = exp(-sqrt(-1)*Wo*t.*k);
  a(k+1) = xt*C(k+1,:)'/length(t); //fourier series is done
  if(abs(a(k+1))<=0.01) 
    a(k+1)=0;
  end
end
a =a';
ak = [-a($:-1:1),a(2:$)];
disp(ak,'Continuous Time Fourier Series Coefficients are:')
//Result
//Continuous Time Fourier Series Coefficients are:   
// column  1 to 11
// 0 0 0 0 -1.010D-18+0.4950495i 0 1.010D-18-0.4950495i 0 0 0 0

Continuous Time Fourier Series of Cosine signal

Senthilkumar ● March 21, 2011 ● Coded in Scilab 
//Continuous Time Fourier Series Spectral Coefficients of 
//a periodic Cosine signal x(t) = cos(Wot)
clear;
close;
clc;
t = 0:0.01:1;
T = 1;
Wo = 2*%pi/T;
xt = cos(Wo*t);
for k =0:5
  C(k+1,:) = exp(-sqrt(-1)*Wo*t.*k);
  a(k+1) = xt*C(k+1,:)'/length(t); //fourier series is done
  if(abs(a(k+1))<=0.01) 
    a(k+1)=0;
  end
end
a =a';
ak = [a($:-1:1),a(2:$)];
disp(ak,'Continuous Time Fourier Series Coefficients are:')
//Result
//Continuous Time Fourier Series Coefficients are:   
// column  1 to 11
// 0 0 0 0 0.5049505+1.010D-18i  0  0.5049505-1.010D-18i  0  0 0 0

2D Linear Convolution for 2D Image Processing

Senthilkumar ● February 27, 2011●2 comments ● Coded in Scilab 
function [y,X,H] = conv2d2(x,h)

[x1,x2] = size(x);
[h1,h2] = size(h);
X = zeros(x1+h1-1,x2+h2-1); 
H = zeros(x1+h1-1,x2+h2-1); 
Y = zeros(x1+h1-1,x2+h2-1); 
for i = 1:x1
    for j = 1:x2
        X(i,j)=x(i,j);
    end
end
for i =1:h1
    for j = 1:h2
        H(i,j)=h(i,j);
    end
end
disp(X,'x=') 
disp(H,'h =')
Y = fft2d(X).*fft2d(H);
disp(Y)
y = ifft2d(Y);
endfunction
/////////////////////////////////////////////
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
/////////////////////////////////////////////////
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

Circular convolution

Senthilkumar ● February 8, 2011●1 comment ● Coded in Scilab 
//Performing Circular COnvolution
//Using DFT
clear all;
clc;
close;
L = 4; //Length of the Sequence
N = 4;  // N -point DFT
x1 = [2,1,2,1];
x2 = [1,2,3,4];
//Computing DFT 
X1 = dft(x1,-1)
X2 = dft(x2,-1)
//Multiplication of 2 DFTs
X3 = X1.*X2
//Circular Convolution Result
x3 =abs(dft(X3,1))
disp(x3,'x3=')

linear convolution in scilab

Senthilkumar ● February 8, 2011●2 comments ● Coded in Scilab 
//Program for Linear Convolution
clc;
clear all;
close ;
x = input('enter x seq');
h = input('enter h seq');
m = length(x);
n = length(h);
//Method : Using Direct Convolution Sum Formula
for i = 1:n+m-1
    conv_sum = 0;
    for j = 1:i
        if (((i-j+1) <= n)&(j <= m))
            conv_sum = conv_sum + x(j)*h(i-j+1);
        end;
        y(i) = conv_sum;
    end;    
end;
disp(y,'y=')