NOTICE: The Processors Wiki will End-of-Life on January 15, 2021. It is recommended to download any files or other content you may need that are hosted on processors.wiki.ti.com. The site is now set to read only.
MexExample
This example shows how to take DSP C code that uses TI C66x intrinsic FFT functions and create a mex file for use in Matlab. This allows you to compare simulation results in Matlab from your DSP code and any straight matlab functions/model you have.
Contents
Library Requirements[edit]
This example makes use of intrinsic functions that are part of TI's DSPLib. Please find the DSPLib as part of the Multicore Software Development Kit (MCSDK). This can be downloaded at:
http://software-dl.ti.com/sdoemb/sdoemb_public_sw/bios_mcsdk/latest/index_FDS.html
or standalone DSPLib at:
http://software-dl.ti.com/sdoemb/sdoemb_public_sw/dsplib/latest/index_FDS.html
After installing, the intrinsic source can be found at:
C:\ti\dsplib_c66x_3_1_x_x\packages\ti\dsplib\src
The functions used for this example are:
DSP_fft16x16
DSP_fft16x16r
DSP_fft16x32
DSP_fft32x32
DSPF_dp_fftDPxDP
DSPF_sp_fftSPxSP
The intrinsic functions from DSPLib version 3.1.1.1 (and later) will compile with host intrinsic libraries. If using earlier versions of the DSPLib, edits to the instrinsic files will be needed (See http://processors.wiki.ti.com/index.php/Mex_old_DSPLIB_Instr)
Tool Requirements[edit]
This flow requires the C6x intrinsic simulator, available at:
http://processors.wiki.ti.com/index.php/Run_Intrinsics_Code_Anywhere.
An environmental variable C6XSIM_INSTALL_DIR indicating its installation directory is needed.
Executing the Example[edit]
This example was tested with Matlab 7.1.2.0 (R2011a) and VC++ 2008 Express Edition
Please note: The code here does not compile with matlab in-built compiler lcc, use mex -setup in matlab to choose a VC++ compiler.
run makeFftDspLib with proper argument flavor to create mex version of the corresponding DSPLIB FFT
runExample.m compares the mexed DSPLib FFT output with the output from the MatLab in-built FFTs.
makeFftDspLib.m[edit]
function []=makeFftDspLib(flavor)
% function to mex the desired FFT
% flavor - the type of FFT to be used
% currently supported flavors
% 'DSP_fft16x16'
% 'DSP_fft16x16r'
% 'DSP_fft16x32'
% 'DSP_fft32x32'
% 'DSP_fft32x32s'
% 'DSP_sp_fftSPxSP'
% 'DSP_dp_fftDPxDP'
%
% get directories
% C6x simulator installation directory
c6xSimDir = getenv('C6XSIM_INSTALL_DIR');
% DSPLIB directiry
dspLibDir = './';
mexDir = './';
defines = '-DTMS320C66X -D_TMS320C6600 -DLITTLE_ENDIAN_HOST -D_LITTLE_ENDIAN';
srcFiles = [c6xSimDir '/C6xSimulator.c ' ...
c6xSimDir '/C66_ag_intrins.c ' ...
c6xSimDir '/C66_data_sim.c ' ...
mexDir '/mexFftDspLib.c '];
includeDir = [' -I' c6xSimDir ' -I' dspLibDir 'packages'];
switch flavor
case 'DSP_fft16x16'
srcDir = [dspLibDir '/packages/ti/dsplib/src/DSP_fft16x16/c66/'];
srcFiles = [srcFiles ...
srcDir 'gen_twiddle_fft16x16.c '...
srcDir 'DSP_fft16x16.c'];
defines = [defines ' -DDSP_FFT16X16'];
outFile = ['DSP_fft16x16'];
case 'DSP_fft16x16r'
srcDir = [dspLibDir '/packages/ti/dsplib/src/DSP_fft16x16r/c66/'];
srcFiles = [srcFiles ...
srcDir 'gen_twiddle_fft16x16r.c '...
srcDir 'DSP_fft16x16r.c'];
defines = [defines ' -DDSP_FFT16X16R'];
outFile = ['DSP_fft16x16r'];
case 'DSP_fft16x32'
srcDir = [dspLibDir '/packages/ti/dsplib/src/DSP_fft16x32/c66/'];
srcFiles = [srcFiles ...
srcDir 'gen_twiddle_fft16x32.c '...
srcDir 'DSP_fft16x32.c'];
defines = [defines ' -DDSP_FFT16X32'];
outFile = ['DSP_fft16x32'];
case 'DSP_fft32x32'
srcDir = [dspLibDir '/packages/ti/dsplib/src/DSP_fft32x32/c66/'];
srcFiles = [srcFiles ...
srcDir 'gen_twiddle_fft32x32.c '...
srcDir 'DSP_fft32x32.c'];
defines = [defines ' -DDSP_FFT32X32'];
outFile = ['DSP_fft32x32'];
case 'DSPF_sp_fftSPxSP'
srcDir = [dspLibDir '/packages/ti/dsplib/src/DSPF_sp_fftSPxSP/c66/'];
srcFiles = [srcFiles ...
srcDir 'DSPF_sp_fftSPxSP_opt.c'];
defines = [defines ' -DDSPF_SP_FFTSPXSP'];
outFile = ['DSPF_sp_fftSPxSP'];
case 'DSPF_dp_fftDPxDP'
srcDir = [dspLibDir '/packages/ti/dsplib/src/DSPF_dp_fftDPxDP/c66/'];
srcFiles = [srcFiles ...
srcDir 'DSPF_dp_fftDPxDP.c'];
defines = [defines ' -DDSPF_DP_FFTDPXDP'];
outFile = ['DSPF_dp_fftDPxDP'];
otherwise
error(['~makeFft.m: flavor not defined for mexing']);
end;
% create the mex command
mexCommand = ['mex -output ' mexDir outFile ' ' srcFiles ...
includeDir ' ' defines];
% disp(mexCommand);
% evlaute mex command
eval(mexCommand);
return
mexFftDspLib.c[edit]
/**
* @file mexFftDspLib.c
* @brief The file contains mex implementation of DSPLib fft funcstions
* only input is input vector (complex or real)
* only output is fft of input data (complex)
* (C) Copyright 2012, Texas Instruments, Inc.
*/
#include "mex.h"
#ifdef DSP_FFT16X16
#include "ti/dsplib/src/DSP_fft16x16/DSP_fft16x16.h"
#include "ti/dsplib/src/DSP_fft16x16/c66/gen_twiddle_fft16x16.h"
#define IOTYPE short
#define TWIDDLETYPE short
#define MAXFFTSIZE 65536
#define MINFFTSIZE 16
#endif
#ifdef DSP_FFT16X16R
#include "ti/dsplib/src/DSP_fft16x16r/DSP_fft16x16r.h"
#include "ti/dsplib/src/DSP_fft16x16r/c66/gen_twiddle_fft16x16r.h"
#define IOTYPE short
#define TWIDDLETYPE short
#define MAXFFTSIZE 65536
#define MINFFTSIZE 16
#endif
#ifdef DSP_FFT16X32
#include "ti/dsplib/src/DSP_fft16x32/DSP_fft16x32.h"
#include "ti/dsplib/src/DSP_fft16x32/c66/gen_twiddle_fft16x32.h"
#define IOTYPE int
#define TWIDDLETYPE short
#define MAXFFTSIZE 65536
#define MINFFTSIZE 16
#endif
#ifdef DSP_FFT32X32
#include "ti/dsplib/src/DSP_fft32x32/DSP_fft32x32.h"
#include "ti/dsplib/src/DSP_fft32x32/c66/gen_twiddle_fft32x32.h"
#define IOTYPE int
#define TWIDDLETYPE int
#define MAXFFTSIZE 65536
#define MINFFTSIZE 16
#endif
#ifdef DSPF_SP_FFTSPXSP
#include "ti/dsplib/src/DSPF_sp_fftSPxSP/c66/DSPF_sp_fftSPxSP_opt.h"
#define IOTYPE float
#define TWIDDLETYPE float
#define MAXFFTSIZE 131072
#define MINFFTSIZE 8
void gen_twiddle_fftspxsp(float *w, int n);
unsigned char brev[64] = {
0x0, 0x20, 0x10, 0x30, 0x8, 0x28, 0x18, 0x38,
0x4, 0x24, 0x14, 0x34, 0xc, 0x2c, 0x1c, 0x3c,
0x2, 0x22, 0x12, 0x32, 0xa, 0x2a, 0x1a, 0x3a,
0x6, 0x26, 0x16, 0x36, 0xe, 0x2e, 0x1e, 0x3e,
0x1, 0x21, 0x11, 0x31, 0x9, 0x29, 0x19, 0x39,
0x5, 0x25, 0x15, 0x35, 0xd, 0x2d, 0x1d, 0x3d,
0x3, 0x23, 0x13, 0x33, 0xb, 0x2b, 0x1b, 0x3b,
0x7, 0x27, 0x17, 0x37, 0xf, 0x2f, 0x1f, 0x3f
};
#endif
#ifdef DSPF_DP_FFTDPXDP
#include "ti/dsplib/src/DSPF_dp_fftDPxDP/c66/DSPF_dp_fftDPxDP.h"
#define IOTYPE double
#define TWIDDLETYPE double
#define MAXFFTSIZE 131072
#define MINFFTSIZE 8
void gen_twiddle_fftdpxdp(double *w, int n);
#endif
/* The gateway routine. */
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[] )
{
double *pInputR, *pInputI, * pOutputR, *pOutputI;
IOTYPE *fftInput, *fftOutput;
TWIDDLETYPE *twiddle;
int fftSize, index, size;
int m, n;
/* get input vector sizes */
m = mxGetM(prhs[0]);
n = mxGetN(prhs[0]);
/* input error checking */
if(m==1)
{
fftSize = n;
}
else if(n==1)
{
fftSize = m;
}
else
{
mexErrMsgTxt("Input must be a vector");
}
/* check for fftSize to be power of 2 */
if((fftSize & (fftSize-1)) != 0)
{
mexErrMsgTxt("size of input vector must be a power of 2");
}
if(fftSize >= MAXFFTSIZE)
{
mexErrMsgTxt("size of input vector exceeds maximum allowed value");
}
if(fftSize <= MINFFTSIZE)
{
mexErrMsgTxt("size of input vector is less than minimum allowed value");
}
/* allocate memory (twice FFT size to accomodate real amd imaginary) */
fftInput = (IOTYPE *)mxCalloc(sizeof(IOTYPE), 2*fftSize);
fftOutput = (IOTYPE *)mxCalloc(sizeof(IOTYPE), 2*fftSize);
twiddle = (TWIDDLETYPE *)mxCalloc(sizeof(TWIDDLETYPE), 2*fftSize);
if(mxIsComplex(prhs[0])) /* input data is complex */
{
pInputR = mxGetPr(prhs[0]);
pInputI = mxGetPi(prhs[0]);
for (index = 0; index < fftSize; index++)
{
/* fill in real data */
fftInput[2*index] = (IOTYPE)(pInputR[index]);
/* fill in imag data */
fftInput[2*index+1] = (IOTYPE)(pInputI[index]);
}
}
else if(mxIsDouble(prhs[0])) /* input data is real */
{
pInputR = mxGetPr(prhs[0]);
for (index = 0; index < fftSize; index++)
{
/* fill in real data */
fftInput[2*index] = (IOTYPE)(pInputR[index]);
/* imaginary part filled with zero */
}
}
else
{
mexErrMsgTxt("Input must be a complex or real");
}
#ifdef DSP_FFT16X16
/* generate twiddle factor */
size = gen_twiddle_fft16x16(twiddle, fftSize);
/* call FFT function */
DSP_fft16x16(twiddle, fftSize, fftInput, fftOutput);
#endif
#ifdef DSP_FFT16X16R
{ int i, j, rad;
j = 0;
for (i = 0; i <= 31; i++)
if ((fftSize & (1 << i)) == 0)
j++;
else
break;
if (j % 2 == 0)
rad = 4;
else
rad = 2;
/* generate twiddle factor */
size = gen_twiddle_fft16x16r(twiddle, fftSize);
/* call FFT function */
DSP_fft16x16r(fftSize, fftInput, twiddle, fftOutput, rad, 0, fftSize);
}
#endif
#ifdef DSP_FFT16X32
/* generate twiddle factor */
size = gen_twiddle_fft16x32(twiddle, fftSize);
/* call FFT function */
DSP_fft16x32(twiddle, fftSize, fftInput, fftOutput);
#endif
#ifdef DSP_FFT32X32
/* generate twiddle factor */
size = gen_twiddle_fft32x32(twiddle, fftSize, pow(2.0, 31.0) - 0.5);
/* call FFT function */
DSP_fft32x32(twiddle, fftSize, fftInput, fftOutput);
#endif
#ifdef DSPF_SP_FFTSPXSP
{ int i, j, rad;
j = 0;
for (i = 0; i <= 31; i++)
if ((fftSize & (1 << i)) == 0)
j++;
else
break;
if (j % 2 == 0)
rad = 4;
else
rad = 2;
/* generate twiddle factor */
gen_twiddle_fftspxsp(twiddle, fftSize);
/* call FFT function */
DSPF_sp_fftSPxSP_opt(fftSize, fftInput, twiddle, fftOutput, brev, rad, 0, fftSize);
}
#endif
#ifdef DSPF_DP_FFTDPXDP
{ int i, j, rad;
j = 0;
for (i = 0; i <= 31; i++)
if ((fftSize & (1 << i)) == 0)
j++;
else
break;
if (j % 2 == 0)
rad = 4;
else
rad = 2;
/* generate twiddle factor */
gen_twiddle_fftdpxdp(twiddle, fftSize);
/* call FFT function */
DSPF_dp_fftDPxDP(fftSize, fftInput, twiddle, fftOutput, rad, 0, fftSize);
}
#endif
/* put data into matlab output */
/* create memory */
plhs[0] = mxCreateDoubleMatrix(m, n, mxCOMPLEX);
/* get memory pointers */
pOutputR = mxGetPr(plhs[0]);
pOutputI = mxGetPi(plhs[0]);
/* data into output memory */
for (index = 0; index < fftSize; index++)
{
pOutputR[index] = (double) fftOutput[2*index];
pOutputI[index] = (double) fftOutput[2*index+1];
}
/* free all allocated memory */
mxFree(fftInput);
mxFree(fftOutput);
mxFree(twiddle);
}
#ifdef DSPF_SP_FFTSPXSP
#include math.h
/* Function for generating Specialized sequence of twiddle factors */
void gen_twiddle_fftspxsp(float *w, int n)
{
int i, j, k;
const double PI = 3.141592654;
for (j = 1, k = 0; j <= n >> 2; j = j << 2)
{
for (i = 0; i < n >> 2; i += j)
{
#ifdef _LITTLE_ENDIAN
w[k] = (float) sin (2 * PI * i / n);
w[k + 1] = (float) cos (2 * PI * i / n);
w[k + 2] = (float) sin (4 * PI * i / n);
w[k + 3] = (float) cos (4 * PI * i / n);
w[k + 4] = (float) sin (6 * PI * i / n);
w[k + 5] = (float) cos (6 * PI * i / n);
#else
w[k] = (float) cos (2 * PI * i / n);
w[k + 1] = (float) -sin (2 * PI * i / n);
w[k + 2] = (float) cos (4 * PI * i / n);
w[k + 3] = (float) -sin (4 * PI * i / n);
w[k + 4] = (float) cos (6 * PI * i / n);
w[k + 5] = (float) -sin (6 * PI * i / n);
#endif
k += 6;
}
}
}
#endif
#ifdef DSPF_DP_FFTDPXDP
/* Function for generating Specialized sequence of twiddle factors */
void gen_twiddle_fftdpxdp(double *w, int n)
{
int i, j, k;
const double PI = 3.141592654;
for (j = 1, k = 0; j <= n >> 2; j = j << 2)
{
for (i = 0; i < n >> 2; i += j)
{
w[k] = cos (2 * PI * i / n);
w[k + 1] = sin (2 * PI * i / n);
w[k + 2] = cos (4 * PI * i / n);
w[k + 3] = sin (4 * PI * i / n);
w[k + 4] = cos (6 * PI * i / n);
w[k + 5] = sin (6 * PI * i / n);
k += 6;
}
}
}
#endif
runExample.m[edit]
% first mex all FFT versions
disp('mexing DSP_fft16x16');
makeFftDspLib('DSP_fft16x16');
disp('mexing DSP_fft16x16r');
makeFftDspLib('DSP_fft16x16r');
disp('mexing DSP_fft16x32');
makeFftDspLib('DSP_fft16x32');
disp('mexing DSP_fft32x32');
makeFftDspLib('DSP_fft32x32');
disp('mexing DSPF_sp_fftSPxSP');
makeFftDspLib('DSPF_sp_fftSPxSP');
disp('mexing DSPF_dp_fftDPxDP');
makeFftDspLib('DSPF_dp_fftDPxDP');
listN=[128 256 512 1024];
for cnt=1:length(listN),
N=listN(cnt);
disp('running test For N =');
disp(N);
% scaled input (to ensure no saturation for 16x16 implementations)
x=fix(2^(15-log(N)/log(4)-1)*rand(N,1));
% Matlab FFT
xfRef=fft(x);
xfOut=DSP_fft16x16(x);
scale=2^(ceil(log(N)/log(4))-1);
disp('DSP_fft16x16 SNR:');
disp(20*log10(norm(scale*xfOut-xfRef)/norm(xfRef)));
xfOut=DSP_fft16x16r(x);
scale=2^(ceil(log(N)/log(4))-1);
disp('DSP_fft16x16r SNR:');
disp(20*log10(norm(scale*xfOut-xfRef)/norm(xfRef)));
xfOut=DSP_fft16x32(x);
scale=1;
disp('DSP_fft16x32 SNR:');
disp(20*log10(norm(scale*xfOut-xfRef)/norm(xfRef)));
xfOut=DSP_fft32x32(x);
scale=1;
disp('DSP_fft32x32 SNR:');
disp(20*log10(norm(scale*xfOut-xfRef)/norm(xfRef)));
xfOut=DSPF_sp_fftSPxSP(x);
scale=1;
disp('DSPF_sp_fftSPxSP SNR:');
disp(20*log10(norm(scale*xfOut-xfRef)/norm(xfRef)));
xfOut=DSPF_dp_fftDPxDP(x);
scale=1;
disp('DSPF_dp_fftDPxDP SNR:');
disp(20*log10(norm(scale*xfOut-xfRef)/norm(xfRef)));
end;
