DSPRelated.com
   Code Snippets

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

Template to develop audio DSP using the (OSS) API + GNU/Linux

Gabriel Rivas ● February 6, 2011 ● Coded in C 
#include <stdio.h>
#include <fcntl.h>
#include <sys/soundcard.h>

#define BUF_SIZE 2

/*buffers for sound device*/
unsigned char devbuf[BUF_SIZE];

/*File descriptor for device*/
int audio_fd;

/*The sound card is known as the DSP*/
char DEVICE_NAME[]="/dev/dsp";

/*Device Settings*/
int format;
int channels;
int fs;
int len;
int frag;
int devcaps;

int main()
{
    unsigned int temp;
    unsigned int yout;

    /*Samples format is 16 bit unsigned*/
    format = AFMT_U16_LE;
    /*1 Channel, MONO*/
    channels = 1;
    /*Sampling Rate is 16 KHz*/
    fs = 16000;
    /*This is a parameter used to set the DSP for low latencies*/
    frag = 0x00020004;

	/******************************************************
          Set parameters in sound card device
    *****************************************************/

	/*Open sound card device*/
	if ((audio_fd = open(DEVICE_NAME,O_RDWR, 0)) == -1) {
		/* Open of device failed */
		perror(DEVICE_NAME);
		exit(1);
	}	

	if (ioctl (audio_fd, SNDCTL_DSP_SETFRAGMENT, &frag) == -1)
	{
	    perror ("SNDCTL_DSP_SETFRAGMENT");
	    exit (-1);
	}

	/*Set audio format*/
	if (ioctl(audio_fd, SNDCTL_DSP_SETFMT, &format) == -1) {
		/* fatal error */
		perror("SNDCTL_DSP_SETFMT");
		exit(1);
	}
	/*Set number of channels*/
	if (ioctl(audio_fd, SNDCTL_DSP_CHANNELS, &channels) == -1) {
		/* Fatal error */
		perror("SNDCTL_DSP_CHANNELS");
		exit(1);
	}
	
	/*Set sampling rate*/
	if (ioctl(audio_fd, SNDCTL_DSP_SPEED, &fs)==-1) {
		/* Fatal error */
		perror("SNDCTL_DSP_SPEED");
		exit(1);
	}

	if (ioctl(audio_fd, SNDCTL_DSP_SYNC) == -1) {
		/* fatal error */
		perror("SNDCTL_DSP_SYNC");
		exit(1);
	}	

	/******************************************************
          This is the infinite loop to:
          1. Read a sample
          2. Process the sample
          3. Write the sample to soundcard output
    *****************************************************/
    while(1) {
        /* This is a blocking read function, the application will wait
           until the sound card captures a sample
        */
        if ((len = read(audio_fd,devbuf, sizeof(devbuf))) == -1) {
	        perror("audio read");
	        exit(1);
	    }

        /* We can pass this variable to a temporary value, 
           in this case as unsigned 16 value, and then use this value 
           for processing
        */  
	    temp = (devbuf[0])|(devbuf[1]<<8);

        /* In this case no processing is done so the value is passed to the output. You can also use this sample to build an audio file, or send it trought a communications interface, etc.
        */
	    yout = temp;

	    devbuf[0] = yout&0xFF;
	    devbuf[1] = (yout>>8)&0xFF;
		
        /* This is the way the output samples are sent to the soundcard            
        */
        if ((len = write(audio_fd,devbuf,sizeof(devbuf)) == -1)) {
	        perror("audio write");
	        exit(1);
        }
    }	
    return 0;
}

TI F281x SCI Configuration

Emmanuel ● February 4, 2011 ● Coded in C for the TI F28x
/**********************************************************************/
#include "DSP281x_Device.h"
/**********************************************************************
* Function: InitSci()
*
* Description:  This function initializes F281x SCI. In this case the desired rate
*				is 9600 bps, 8 databits, 1 stop bit and no parity.
**********************************************************************/
void InitSci(void)
{

/*** SCI reset	*/

	ScibRegs.SCICTL1.bit.SWRESET= 0;
	
/*** Setting the baud rate to 9600 bps assuming LSPCLK=150 Mhz */

	ScibRegs.SCIHBAUD = 0x07;		
	ScibRegs.SCILBAUD = 0xA0;
	

/*** Configuration for no stop bits, no parity and 8 data bits		*/

	ScibRegs.SCICCR.all = 0x07;
	
// bit 7 		0:		STOP BITS
// bit 6 		0:		PARITY
// bit 5 		0:		PARITY ENABLE
// bit 4 		0:		LOOP BACK ENA
// bit 3 		0:		ADDR/IDLE MODE
// bit 2-0 		111:	SCI CHAR

/*** Transmitter and receiver enabled	*/

	ScibRegs.SCICTL1.all = 0x03;
	
// bit 7 		0:		Reserved
// bit 6 		0:		RX ERR INT
// bit 5 		0:		SW RESET
// bit 4 		0		Reserved
// bit 3 		0:		TXWAKE
// bit 2 		0:		SLEEP
// bit 1 		1:		TXENA, Habilitación del transmisor
// bit 0 		1:		RXENA, Habilitación del receptor

/*** Tx and Rx interruption enabled 	*/

	ScibRegs.SCICTL2.all = 0x0003;
	
// bit 7 		0:		TXRDY
// bit 6 		0:		TX EMPTY
// bit 5-2 		0's:	Reserved
// bit 1 		1:		RX/BK INT
// bit 0 		1:		TX INT
	
/***Emulation mode configuration	*/

	ScibRegs.SCIPRI.all = 0x0010;

// bit 7-5 		0's:	Reserved
// bit 4-3 		10:		Emulation Mode
// bit 2-0 		0's:	Reserved

/*** SCI Reset		*/

	ScibRegs.SCICTL1.bit.SWRESET= 1;
	
/*** SCI interruptions enabled in PIE */

	PieCtrlRegs.PIEIER9.bit.INTx4= 1;	  
	PieCtrlRegs.PIEIER9.bit.INTx3= 1;	                
    IER |= 0x0100;              						

}

/**End of file*/

Frequency Estimate using FFT

● February 2, 2011●1 comment ● Coded in Matlab 
%For a real-valued, single tone sine or cosine wave, estimate the frequency 
%using an fft. Resolution will be limited by the sampling rate of the tone
%and the number of points in the fft.

%Procedure:
%1. Generate the test tone with a given Fs and N
%2. Take the fft and keep the first half
%3. Detect the maximum peak and its index
%4. Equate this index to a frequency

%Error Estimate:
% (Fs/N)= Hz/bin

%Things to keep in mind:
%Adding noise to the input will skew the results
%windowing will help with the fft

%% Clear and close everything
clc;
close all;
clear all; %remove this line if you are trying to use breakpoints

%% Just copy into m file, save and run
Fs = 1e6; %1MHz
fi = 1.333e3; %1kHz
t = 0:1/Fs:.5;
y = sin(2*pi*fi*t);

%Plot the time and frequency domain. Be sure to zoom in to see the waveform
%and spectrum
figure;
subplot(2,1,1);
temp = (2/fi)*Fs;
plot(y);
xlabel('time (s)');
subplot(2,1,2);
sX=fft(y);
N=length(sX);
sXM = abs(sX(1:floor(end/2))).^2; %take the magnitude and only keep 0:Fs/2
plot(0:Fs/N:(Fs/2)-1,sXM);
xlabel('Frequency')
ylabel('Magnitude')

[vv, ii]=max(sXM); %find the index of the largest value in the spectrum

freqEst = (ii-1)*Fs/N; %units are Hz
resMin = Fs/N; %units are Hz

display(freqEst);%display to the command window the results
display(resMin);

TMX Transmultiplexer (MATLAB Version)

David Valencia ● February 2, 2011●6 comments ● Coded in Matlab 
% UPIITA IPN 2010
% Procesamiento Digital de Señales
% Grupo 6TM1  Equipo #8
% Team members
%   Rogelio Baeza Nieves
%   Ponce Mosso Catarina
%   Valencia Pesqueira José David

% TMX Transmultiplexer -> Chain processing

close all;clear all;clc;

n_stages = 2;
n_branches = 2^n_stages;
dec_factor = 2^n_stages;
niter = 500;

N = 16;
input = zeros(2^n_stages,N);
for(k = 0:2^n_stages-1)
    input(k+1,:) = (1 + k * N ):((k+1)*N);
end
input;
in_cnt = 1;

typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db2';

if(typeofbasis == 'b')
    [Rf,Df] = biorwavf(typbior);
    [h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
    [h0,h1,g0,g1] = wfilters(typor);
end;

L = length(h0);

try
    sinc_factor = getsinc(n_stages,L) - 1;
catch
    disp('Error');
    return
end

for i = 1:1
hx = fliplr(formafiltrosdwpt(n_stages,2^n_stages,h0,h1));
Lhx = length(hx);
H = zeros(2^n_stages, Lhx);
for i=0:(2^n_stages)-1
    H(i+1,:)=fliplr(formafiltrosdwpt(n_stages,2^n_stages-i,h0,h1));
end
La = length(H(1,:));

G = zeros(2^n_stages, Lhx);
for i=0:(2^n_stages)-1
    G(i+1,:)=fliplr(formafiltrosdwpt(n_stages,2^n_stages-i,g0,g1));
end
Ls = length(G(1,:));

double2cfloat(H,'h', G,'g',sinc_factor);

analysisbuffer = zeros(1,La);

y = zeros(2^n_stages,1);

sbuffer = zeros(2^n_stages, Ls);
xnbuff = zeros(2^n_stages,1);
chanbuffer = [0];
outputbuffer = zeros(2^n_stages,1);
end

for i=1:niter %General iteration counter
    
    %% Synthesis stage
    i
    % Vector shifting (CLK/1)
    for j = 1:(2^n_stages)
        sbuffer(j,:) = recorreder(sbuffer(j,:),1);
    end
    
% Input interpolation
    if (mod(i,2^n_stages) == 1)&&(in_cnt<N-1)
        sbuffer(:,1) = input(:,in_cnt);
        in_cnt = in_cnt + 1;
    else
        sbuffer(:,1) = zeros(2^n_stages,1);
    end
    
    disp('inputs: ');
    disp(sbuffer)
    
    xnbuff = zeros(2^n_stages,1);
    for j=0:(2^n_stages)-1
        for k = 1 : Ls %convolution
            xnbuff(j+1,1) = xnbuff(j+1,1) + sbuffer(j+1,k)*G(j+1,k);
        end
    end
    
    chanbuffer = sum(xnbuff)   
    
    %% Analysis stage
    analysisbuffer = recorreder(analysisbuffer,1);    
    analysisbuffer(1)=chanbuffer;
    analysisbuffer;
   
    if (i>sinc_factor && mod(i-sinc_factor,2^n_stages) == 1)
    % If the counter is a multiply of the decimating factor, then compute
        y = zeros(2^n_stages,1);
        for j = 0:(2^n_stages)-1
            for k=1:La %convolución
                y(j+1,1) = y(j+1,1) + analysisbuffer(k)*H(j+1,k);
            end;
        end
        outputbuffer = y;
    end
    input;
    i
    outputbuffer
    pause;
end
  	return;

TMX formafiltrosdwpt.m - Integrated version

David Valencia ● February 2, 2011 ● Coded in Matlab 
function [hx] = formafiltrosdwpt(n_stages,branch,h0,h1)
p = branch;

% Integrated version for DWPT/TMX filter generation

hx = 0;
hx(1) = 1;

switch n_stages
    case 1
        if mod(branch,2) ~= 0
            hx = h0;
        else
            hx = h1;
        end
    case 2
        switch branch
            case 1
                hx = conv(h0,upsample(h0,2));
            case 2
                hx = conv(h0,upsample(h1,2));
            case 3
                hx = conv(h1,upsample(h0,2));
            case 4
                hx = conv(h1,upsample(h1,2));
            otherwise
                beep;
                fprintf('\nERROR');
        end
        
    otherwise
        
        for i=0:n_stages-2
            q = floor(p /(2^(n_stages-1-i)));
            if (q == 1)
                hx = conv(hx,upsample(h1,2^i));
            else
                hx = conv(hx,upsample(h0,2^i));
            end
            p = mod(p,2^(n_stages-1-i)); %For DWPT and TMX
        end
        
        t = mod(branch,2);
        if(t == 1)
            hx = conv(hx,upsample(h0,2^(n_stages-1)));
        else
            hx = conv(hx,upsample(h1,2^(n_stages-1)));
        end
        
        
end

db2_2stages.h - TMX Filter data in C

David Valencia ● February 2, 2011 ● Coded in C 
// Beginning of file

#define sinc 3 // sincronization factor
#define N 11 // Filter length
#define M 4 // Number of filters(branches)

float h[M][N]={
{
0.0000e+000,
1.6747e-002,
2.9006e-002,
-7.9247e-002,
1.1274e-001,
-2.9575e-001,
-7.9247e-002,
7.6226e-001,
-2.9575e-001,
-4.0401e-001,
2.3325e-001,
},
{
0.0000e+000,
-6.2500e-002,
-1.0825e-001,
2.9575e-001,
-4.2075e-001,
6.7075e-001,
-4.5425e-001,
2.0425e-001,
-7.9247e-002,
-1.0825e-001,
6.2500e-002,
},
{
0.0000e+000,
-6.2500e-002,
-1.0825e-001,
-1.3726e-001,
-1.7075e-001,
3.5377e-001,
7.2877e-001,
-4.5753e-002,
-5.1226e-001,
-1.0825e-001,
6.2500e-002,
},
{
0.0000e+000,
2.3325e-001,
4.0401e-001,
5.1226e-001,
6.3726e-001,
2.9575e-001,
7.9247e-002,
-1.2260e-002,
-1.3726e-001,
-2.9006e-002,
1.6747e-002,
},
};

float g[M][N]={
{
0.0000e+000,
2.3325e-001,
-4.0401e-001,
-2.9575e-001,
7.6226e-001,
-7.9247e-002,
-2.9575e-001,
1.1274e-001,
-7.9247e-002,
2.9006e-002,
1.6747e-002,
},
{
0.0000e+000,
6.2500e-002,
-1.0825e-001,
-7.9247e-002,
2.0425e-001,
-4.5425e-001,
6.7075e-001,
-4.2075e-001,
2.9575e-001,
-1.0825e-001,
-6.2500e-002,
},
{
0.0000e+000,
6.2500e-002,
-1.0825e-001,
-5.1226e-001,
-4.5753e-002,
7.2877e-001,
3.5377e-001,
-1.7075e-001,
-1.3726e-001,
-1.0825e-001,
-6.2500e-002,
},
{
0.0000e+000,
1.6747e-002,
-2.9006e-002,
-1.3726e-001,
-1.2260e-002,
7.9247e-002,
2.9575e-001,
6.3726e-001,
5.1226e-001,
4.0401e-001,
2.3325e-001,
},
};

// EOF

getsinc.m - TMX Transmultiplexer (MATLAB version)

David Valencia ● February 2, 2011 ● Coded in Matlab 
function [sf] = getsinc(n_etapas, L)

% Experimentally obtained synchronization factors

if(n_etapas == 1)
    sf = 2;

elseif (n_etapas == 2)
    sf = L;

elseif L == 4    
    if(n_etapas == 3) % 8 branches
        sf = 22;
    elseif(n_etapas == 4) % 16 branches
        sf = 36;
    elseif(n_etapas == 5) % 32 branches
        sf = 82;
    end
    
elseif L == 6
    if n_etapas == 3 % 8 branches
        sf = 20;
    elseif n_etapas == 4 % 16 branches
        sf = 30;
    elseif n_etapas == 5 % 32 branches
        sf = 80;
    end
    
elseif L == 8
    if(n_etapas == 3) % 8 branches
        sf = 26;
    elseif(n_etapas == 4) % 16 branches
        sf = 32;        
    end
    
elseif L == 10
    if(n_etapas == 1)
        sf = 2;
    elseif(n_etapas == 2)
        sf = 10;
    end
end

double2cfloat.m - TMX Filter formatting for C compilers

David Valencia ● February 2, 2011 ● Coded in Matlab 
% UPIITA IPN 2010
% José David Valencia Pesqueira

% This program is useful to convert numbers in MATLAB's Double
% format to a text file in float format for C compilers

% You should copy the resulting text to a file and save it with .h
% extension

% THIS PROGRAM RECEIVES AN UNIDIMENTIONAL VECTOR OF DOUBLE NUMBERS

% close all; clear all; clc;

function double2cfloat(x , namex, y, namey, sinc_factor)
[m,n] = size(x);
fprintf('// Begin of file\n\n');
fprintf('#define sinc %d // Sincronization factor\n',sinc_factor);
fprintf('#define N %d // Filter length\n',n);
if m~=0
    fprintf('#define M %d // Number of filters(branches)\n\n',m);
    fprintf('float %s[M][N]=',namex);    
else
    fprintf('float %s[N]=',namex);
end
disp('');
disp('{');
for i=0:m-1 %Rows
    if m~=1
        disp('{');
    end
    for j=0:n-1 %columns
        fprintf('%1.4e,\n',x(i+1,j+1));
    end
    if m~=1 
        fprintf('},\n');
    end
end
disp('};');
fprintf('\n\n');

if m~=0
    fprintf('float %s[M][N]=',namey);    
else
    fprintf('float %s[N]=',namey);
end
disp('');
disp('{');
for i=0:m-1 %Rows
    if m~=1
        disp('{');
    end
    for j=0:n-1 %columns
        fprintf('%1.4e,\n',y(i+1,j+1));
    end
    if m~=1 
        fprintf('},\n');
    end
end
disp('};');

fprintf('\n// EOF\n\n');

Interrupt.c - DSK6713 External Interruptions via GPIO

David Valencia ● February 2, 2011●1 comment ● Coded in C for the TI C67x
/* Interruption.c 
JOSE DAVID VALENCIA PESQUEIRA
UPIITA-IPN 
THIS PROGRAM DETECTS AN EXTERNAL INTERRUPTION ON A GPIO PIN*/ 

#include <stdio.h> 
#include <string.h> 
#include <stdlib.h> 
#include <time.h> 
#include <csl_gpio.h> 
#include <csl_gpiohal.h> 
#include <csl_irq.h> 

/* GLOBAL VARIABLES
The flags are marked as volatile, to tell the compiler
that these can be modified by an external event, avoiding
to put them on registers and rather putting them on RAM
so they're ready for immediate use
*/

volatile int startflag = 0; 

GPIO_Handle gpio_handle;  /* Handle for the GPIO */ 

// GPIO Registers configuration

GPIO_Config gpio_config = {          
	0x00000000,  
	// gpgc = Interruption passtrhrough mode and direct GPIO control mode
	0x0000FFFF, // gpen = All GPIO 0-15 pins enabled
	0x00000000, // gdir = All GPIO pins as inputs
	0x00000000, // gpval = Stores logical level of pins
	0x00000010, // IRQ enable for pin 4
	0x00000010, // Enable Event for pin 4
	0x00000000  // gppol -- default state  
}; 

irq_ext_enable(){ 
	/* First, globally disable interruptions
	in the CSR (Control Status Register).
	Bit 0 is the GIE (Global Interrupt Enable)*/
	CSR = (CSR)&(0xFFFFFFFE); 
	// Enable NMIE bit in the IER (bit 1)
	IER |= 0x00000002; 
	// Enable INT4 bit in the IER (4th bit)
	IER |= 0x00000010; 
	// Finally, globally enable interruptions in the CSR
	CSR |= 0x00000001; 
} 

main() 
{      
	// To configure the GPIO
	gpio_handle = GPIO_open( GPIO_DEV0, GPIO_OPEN_RESET ); 
	GPIO_config(gpio_handle,&gpio_config); 
	irq_ext_enable(); 
	comm_intr();      //init DSK 
	DSK6713_LED_init(); 
	while(startflag == 0);  
	printf(“La interrupción externa encendió el led 3 del dsk6713\n”); 
	while(1) // Infinite While
}    

interrupt void c_int04()           //ISR 
{ 
	startflag = 1; 
	iter = 0; 
	DSK6713_LED_on(0);   //To turn on the led
}