DSPRelated.com
   Code Snippets

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

Random Number Generator

● February 27, 2011●2 comments ● Coded in C for the SHARC
// Random number generator that sounds pleasant in audio algorithms.
// The argument and return value is a 30-bit (positive nonzero) integer.
// This is Duttaro's 30-bit pseudorandom number generator, p.124 J.AES, vol 50 no 3 March 2002; 
// I found this 30-bit pseudorandom number generator sounds far superior to the 31-bit 
// pseudorandom number generator presented in the same article.
// To make result a signed fixed-point value, do not shift left into the sign bit;
// instead, subtract 0x20000000 then multiply by a scale factor.
#define NextRand(rr) (((((rr >> 16) ^ (rr >> 15) ^ (rr >> 1) ^ rr) & 1) << 29) | (rr >> 1))

TI F281x and F2833x ADC Configuration

Emmanuel ● February 17, 2011 ● Coded in C for the TI F28x
/**********************************************************************/
#include "DSP281x_Device.h"		//if using F2812
//#include "DSP2833x_Device.h" //if using F28335
/**********************************************************************
void InitAdc(void)
{
	
/*** ADC Reset*/

	AdcRegs.ADCTRL1.bit.RESET = 1;

	asm(" RPT #22 || NOP");		//Wait 2 clocks assuming ADCCLK=12.5 Mhz and SYSCLKOUT=150 Mhz
	
																			
/*** ADC Power-On	*/

	AdcRegs.ADCTRL3.all = 0x00ED;
			
// bit 15-9      0's:    reserved
// bit 8		 0:		 EXTREF
// bit 7-6       11:     ADCBGRFDN, Reference voltage, 00=off, 11=on
// bit 5         1:      ADCPWDN, Main Power Source, 0=off, 1=on
// bit 4-1       0110:   ADCCLKPS, clock prescaler, FCLK=HSPCLK/(2*ADCCLKPS)
// bit 0         1:      SMODE_SEL, 0=sequencial sampling, 1=simultaneous sampling

	DelayUs(5000);			//TI's asm function	
							//Wait 5 ms

/*** ADCMAXCONV	Register configuration*/

	AdcRegs.ADCMAXCONV.all = 0x0000;
	
// bit 15-7      0's:    reserved
// bit 6-4       000:    MAX_CONV2 
// bit 3-0       0000:   MAX_CONV1  (0 = 1 conversion)

/*** Select channels for sampling in SEQ1	*/

	AdcRegs.ADCCHSELSEQ1.bit.CONV00 = 0;	//Channels A0 and B0
		
/***  ADCTRL1 Register configuration	*/

	AdcRegs.ADCTRL1.all = 0x0700;
	
// bit 15        0:      reserved
// bit 14        0:      RESET
// bit 13-12     00:     SUSMOD, 00=ignore emulation suspend
// bit 11-8      0111:   ACQ_PS Acquisition window pre-scaler
// bit 7         0:      CPS Core clock
// bit 6         0:      CONT_RUN
// bit 5         0:      SEQ_OVRD
// bit 4         0:      SEQ_CASC
// bit 3-0       0000:   reserved

/*** ADCTRL2 Register configuration	*/

	AdcRegs.ADCTRL2.all = 0x0800;
	
// bit 15        0:      EVB_SOC_SEQ: SOC trigger by EVB
// bit 14        0:      RST_SEQ1: SEQ1 reset
// bit 13        0:      SOC_SEQ1: trigger signal for SEQ1
// bit 12        0:      reserved
// bit 11        1:      INT_ENA_SEQ1: SEQ1 interruption enable
// bit 10        0:      INT_MOD_SEQ1: interruption mode
// bit 9         0:      reserved
// bit 8         0:      EVA_SOC_SEQ1, 1=SEQ1 starts with EVA_SOC trigger
// bit 7         0:      EXT_SOC_SEQ1, 1=SEQ1 starts with ADCSOC input signal
// bit 6         0:      RST_SEQ2
// bit 5         0:      SOC_SEQ2
// bit 4         0:      reserved
// bit 3         0:      INT_ENA_SEQ2: SEQ1 interruption enabled
// bit 2         0:      INT_MOD_SEQ2: interruption mode
// bit 1         0:      reserved
// bit 0         0:      EVB_SOC_SEQ2

/*** ADC interruption enable	*/
PieCtrlRegs.PIEIER1.bit.INTx6= 1;	                 
IER |= 0x0001;              						

}

CRC calculation

Emmanuel ● February 15, 2011 ● Coded in C 
unsigned char crc(unsigned char data[],int L)
{
	
    unsigned char P=0x83;
    unsigned char *data_ptr;
	unsigned char crc_result;
	int bit;
	crc_result=0;
	for(data_ptr=data;data_ptr<data+L;data_ptr++)
	{
		for(bit=7;bit>=0;bit--)
		{
            crc_result=(crc_result<<1)+((*data_ptr&0x1<<bit)>>bit);
			if (crc_result>=128)
			{
				crc_result=crc_result^P;
			}
		}
	}
	//multiply by 2^n-k;
	for(bit=6;bit>=0;bit--)
		{
            crc_result=(crc_result<<1);
			if (crc_result>=128)
			{
				crc_result=crc_result^P;
			}
		}

	return crc_result;
		
}

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

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