|
Hello Martin I was hoping I could find or remember the authors name, but I cant. I also saw the post for the other filter generation tool and suggest looking into it (as I will). In the meantime, here is something new that I have prepared for the new DSK tools that may be of interest or help. Basically it configures and runs an SDFT filter collecting its impulse response for N samples. The idea(s) being that - Using the SDFT to generate FIR coefficients is a new method - This code can be embedded into the DSP If anyone has seen this code before, this version has a new feature to it. By fully convolving the window coefficients in the F domain, the generated filters are no longer restricted to being brick wall filters. That is, complex shapes are now allowed. Best regards, Keith Larson /*======================================================================= FIR_C3K.C Keith Larson TMS320 DSP Applications (C) Copyright 1995-2003 Texas Instruments Incorporated This is unsupported code with no implied warranties or liabilities. See the C3x DSK disclaimer document for details This application shows how to interface the PCM3003 on the VC33 DSK to the TMS320VC33. L/R Channels are supported ========================================================================*/ #include <math.h> #include "C3MMR.H" #define pi 3.14159265 #define TIM0_prd 3 // AIC reference clock is TIM0 #define S0xctrl 0x00000111 // #define S0rctrl 0x00000111 // #define bigval 000010000h // Used in overflow mode saturation #define XCLKSRCE ( 0 << 6 ) // 0=external 1= internal #define RCLKSRCE ( 0 << 7 ) // #define XVAREN ( 0 << 8 ) // VAREN 0 1 #define RVAREN ( 0 << 9 ) // FS ____---------_____ #define XFSM ( 0 <<10 ) // 0=burst 1=continuous #define RFSM ( 0 <<11 ) // #define CLKXP ( 1 <<12 ) // #define CLKRP ( 1 <<13 ) // #define DXP ( 0 <<14 ) // #define DRP ( 0 <<15 ) // #define FSXP ( 1 <<16 ) // #define FSRP ( 1 <<17 ) // #define XLEN ( 11b<<18 ) // 00=8 01 10$ 112 #define RLEN ( 11b<<20 ) // #define EXTINT ( 0 <<22 ) // #define EXINT ( 1 <<23 ) // #define ERTINT ( 0 <<24 ) // #define ERINT ( 1 <<25 ) // #define SPBITS 32 #define BITSEL ((SPBITS/8)-1) // 0=8, 1, 2$ and 32 bits #define LP_STATUS (volatile long *) 0x80A000 // pointer CPLD status bits unsigned long S0gctrl_rst = XCLKSRCE|RCLKSRCE |XVAREN|RVAREN |XFSM|RFSM |CLKXP|CLKRP |DXP|DRP |FSXP|FSRP |EXTINT|EXINT|ERTINT|ERINT |(BITSEL<<18)|(BITSEL<<20); unsigned long S0gctrl_run = 0x0C000000 |XCLKSRCE|RCLKSRCE |XVAREN|RVAREN |XFSM|RFSM |CLKXP|CLKRP |DXP|DRP |FSXP|FSRP |EXTINT|EXINT|ERTINT|ERINT |(BITSEL<<18)|(BITSEL<<20); /*=================================================================== The ADC data is read and buffered here ====================================================================*/ #define Max_FIR_size 256 #define FIR_size 256 volatile float LEFT_OUT; volatile float LEFT_IN; volatile float RIGHT_OUT; volatile float RIGHT_IN; volatile long ch; //================================================== void InitPeripherals(void) { long T; *T0_ctrl = 0; // Halt TIM0 *T0_count = 0; // Set counts to 0 *T0_prd = TIM0_prd; // Set period *T0_ctrl = 0x2C1; // Restart both timers // - - - - - - - - - - - - - - - - - - - - - *S0_gctrl = S0gctrl_rst; // global control *S0_xctrl = S0xctrl; // transmit control *S0_rctrl = S0rctrl; // receive control *S0_xdata = 0; // DXR data value *S0_gctrl = S0gctrl_run; // global control T = *S0_rdata; *S0_xdata = T; // DXR data value asm(" or 080h,ST"); // Set overflow mode for fast saturate } //================================================== float FastRand(void) { static unsigned long A=0; float f; asm(" push ST "); asm(" andn 080h,ST "); // clear OVM bit A *= 0x107465L; A += 0x234567L; f = A & 0xFFFFF; f -= 0x80000; asm(" pop ST "); return f; } volatile float CTRLA_VAL[]={0.5, 0.000}; // nom value volatile float CTRLA_MAX[]={2.0, 0.050}; // max gain volatile float CTRLA_MIN[]={0.0, 0.000}; // min gain //================================================================ // Inline assembler versions of the FIR filter //================================================================ float FIR_R(float); // be sure to give a prototypes! float FIR_L(float); // int FIR_SIZE = FIR_size ; int FIR_SIZE2 = FIR_size-2; #if 1 asm("COEFADDR .word _COEF ; "); asm("RDATAADDR .word _D_RIGHT ; "); asm("_FIR_R: ldi @RDATAADDR,AR1 ; "); asm(" ldi @_FIR_SIZE,BK ; "); asm(" stf R0,*AR1++(1)% ; save input data "); asm(" ldi @COEFADDR,AR0 ; "); asm(" ldf 0,R2 ; initialize sum=0 "); asm(" mpyf3 *AR0++,*AR1++(1)%,R0 ; FIR MAC operation "); asm(" rpts @_FIR_SIZE2 ; "); asm(" mpyf3 *AR0++,*AR1++(1)%,R0 ; FIR MAC operation "); asm(" || addf3 R0,R2,R2 ; "); asm(" addf r2,r0 ; "); asm(" sti AR1,@RDATAADDR ; "); asm(" rets ; "); asm("LDATAADDR .word _D_LEFT "); asm("_FIR_L: ldi @LDATAADDR,AR1 ; "); asm(" ldi @_FIR_SIZE,BK ; "); asm(" stf R0,*AR1++(1)% ; save input data "); asm(" ldi @COEFADDR,AR0 ; "); asm(" ldf 0,R2 ; initialize sum=0 "); asm(" mpyf3 *AR0++,*AR1++(1)%,R0 ; FIR MAC operation "); asm(" rpts @_FIR_SIZE2 ; "); asm(" mpyf3 *AR0++,*AR1++(1)%,R0 ; FIR MAC operation "); asm(" || addf3 R0,R2,R2 ; "); asm(" addf r2,r0 ; "); asm(" sti AR1,@LDATAADDR ; "); asm(" rets ; "); #else //====================================================================== // An FIR filter written in C cannot fully use the C3x // circular addressing hardware and is therefor less efficient //====================================================================== float FIR_L(float f) { int n; static int i = 0; float *data=D_LEFT; float *coef=COEF; data[i]=f; f=0; for(n=0;n<FIR_size;n++) { f+=(*coef++ * data[i]); i++; if(i>=FIR_size)i=0; } i++; if(i>=FIR_size)i=0; return f; } float FIR_R(float f) { int n; static int i = 0; float *data=D_RIGHT; float *coef=COEF; data[i]=f; f=0; for(n=0;n<FIR_size;n++) { f+=(*coef++ * data[i]); i++; if(i>=FIR_size)i=0; } i++; if(i>=FIR_size)i=0; return f; } #endif //============================================================ // void sdft_gen_coef(void); // // Frequency and Time Domain Filter Coefficients // --------- // An SDFT filter is used to invert the frequency domain filter coefficients // in FREQ_RESP[] to a time domain set of FIR coefficients in COEF[]. This // is accomplished by initializing the bin data to that of an impulse and // then allowing the SDFT to run for N samples. The SDFT filter output // is therefor the FIR filter or impulse response. You can define the FIR // filter response in the tables below where each value is the response at // a particular frequency // // NOTES // ----- // Control bank CTRLA_VAL[0..1] is used to set the gain on the ADC analog // input jack/microphone [0], and random white noise[1] (internal function). // // The SDFT past history buffer is not required since the only input is // a single impulse at T=0 // // Time domain windowing can also be accomplished in the frequency domain // using convolution. Interestingly in a bandpass this becomes a // multiplication of the frequency domain coefficient endpoints by 0.5 // for a Raised Cosine window (it takes two multiplies at each endpoint // for most other window functions). The remaining 'center' bins are // untouched. This considerably improves band stop rejection. // // If you have two DSK's use one to filter a test tone and the other to // perform an FFT spectrum analysis using DSK3WIN3.EXE. DSK3WIN3.EXE has // the ability to average signals greatly improving side visability in the // presence of white noise. The F11 key increases and F12 decreases averaging // // The inverse Fourier Transform method does not always produce the most // optimal coefficients. Milage will vary! // // See the SDFT documentation in DSK3HELP.HLP for further details // (hit the F1 key if you are using DSK3DW.EXE) //----------------------- extern float D_LEFT []; extern float D_RIGHT[]; asm("_D_LEFT .sect \"FIR_LEFT\" "); // place on 2^N boundary! asm("_D_RIGHT .sect \"FIR_RIGHT\" "); // place on 2^N boundary! float COEF [Max_FIR_size ]; #define DIFFERENTIATOR 0 #define INTEGRATOR 0 #define BANDPASS 1 #define WINCO 0.5 //--------------- // gain is proportional to F //--------------- #if DIFFERENTIATOR float FREQ_RESP[(Max_FIR_size/2)+2] ={ 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, 1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95, 2.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 #if Max_FIR_size=%6 0.00, 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, 0.95, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45, 1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.95, 2.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 #endif ,0,0 }; #endif //--------------------- // Gain is proportional to (N-f) //--------------------- #if INTEGRATOR float FREQ_RESP[(Max_FIR_size/2)+2] ={ 3.15, 3.10, 3.05, 3.00, 2.95, 2.90, 2.85, 2.80, 2.75, 2.70, 2.65, 2.60, 2.55, 2.50, 2.45, 2.40, 2.35, 2.30, 2.25, 2.20, 2.15, 2.10, 2.05, 2.00, 1.95, 1.90, 1.85, 1.80, 1.75, 1.70, 1.65, 1.60, 1.55, 1.50, 1.45, 1.40, 1.35, 1.30, 1.25, 1.20, 1.15, 1.10, 1.05, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00 #if Max_FIR_size=%6 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 #endif ,0,0 }; #endif //---------------------------- // A multi band pass is created with one band increasing in output // with increasing frequency and the other flat, but at -20db //---------------------------- #if BANDPASS float FREQ_RESP[(Max_FIR_size/2)+2] ={ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.60, 0.80, 1.00, 1.20, 1.40, 1.60, 1.80, 2.00, 2.20, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.12, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.12, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 #if Max_FIR_size=%6 ,0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 2.20, 2.00, 1.80, 1.60, 1.40, 1.20, 1.00, 0.80, 0.60, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 #endif ,0,0 }; #endif //--------------------------- // Initialize SDFT twiddles and SDFT impulse // // The starting impulse in each bin can be of any phase angle. An // angle of 90' is interesting since this is the Hilbert transform. // Note that the resulting FIR coefficients are symetric and will // have a group delay of N/2 samples and that any phase can be // achieved using the sin()/cos() phase adjustment // // Another feature that is rolled in is the +/-/+/- phase shifting // that occurs during signal reconstruction. It is instead // incorporated into the input as 'flipsign' to make signal // reconstruction easier // // Two versions of the function are also given. In the first, // a twiddle table is used speeding up the calculation. In the // second the twiddles are computed on the fly, decreasing the // amount of data that must be allocated. //--------------------------- //#define PHASE (pi/2.0) // IMAG 90' phase (HILBERT) #define PHASE 0 // REAL 0' phase #define FR FREQ_RESP #if 0 // fast & compact code //-------------------------- // This version computes the twiddles to an array which is then // accessed by the SDFT. This would result in a faster SDFT, // especially if the twiddle calcalation were seperated from // the main SDFT routine. The cost is a larger data footprint //-------------------------- void sdft_gen_coef(void) { int k, t; float c,f,inv_FIR_size,r; float flipsign=1; float dr[Max_FIR_size/2]; // SDFT bin data float di[Max_FIR_size/2]; float tr[Max_FIR_size/2]; // SDFT bin twiddles float ti[Max_FIR_size/2]; inv_FIR_size = 2.0/FIR_size; for(k=0;k<FIR_size/2;k++) { tr[k] = cos(k*2*pi/FIR_size); ti[k] = sin(k*2*pi/FIR_size); dr[k] = flipsign * cos(PHASE); di[k] = flipsign * sin(PHASE); flipsign*=-1; } // Advance the SDFT filter in time generating an impulse response for(t=0;t<FIR_size;t++) { c = 0; r=0; for(k=0;k<(FIR_size/2);k++) { f = (dr[k]*tr[k])-(di[k]*ti[k]); di[k] = (dr[k]*ti[k])+(di[k]*tr[k]); dr[k] = f; c+= (0.5*FR[k]+FR[k+1]+0.5*FR[k+2])*dr[k]; } COEF[t]= c * inv_FIR_size; } } #else //-------------------------- // This version computes the twiddles as they are needed. // This gives a smaller data footprint at the cost of speed, // which in this case would be OK since this routine is // called only once on startup //-------------------------- void sdft_gen_coef(void) { int k, t; float c,f,inv_FIR_size,r,twr,twi,twr_step,twi_step; float flipsign=1; float dr[Max_FIR_size/2]; // SDFT bin data float di[Max_FIR_size/2]; inv_FIR_size = 2.0/FIR_size; twr_step = cos(2*pi/FIR_size); twi_step = sin(2*pi/FIR_size); for(k=0;k<FIR_size/2;k++) { dr[k] = flipsign * cos(PHASE); di[k] = flipsign * sin(PHASE); flipsign*=-1; } //--------------------------- // Advance the SDFT filter in time generating an impulse response //--------------------------- for(t=0;t<FIR_size;t++) { c = 0; twr = 1; twi = 0; r=0; for(k=0;k<FIR_size/2;k++) { // twr = cos(k*2*pi/FIR_size); // compute twiddle as needed // twi = sin(k*2*pi/FIR_size); // //----- f = (dr[k]*twr)-(di[k]*twi); // compute SDFT bins di[k] = (dr[k]*twi)+(di[k]*twr); // dr[k] = f; // c+= (0.5*FR[k]+FR[k+1]+0.5*FR[k+2])*dr[k]; // and sum the output //--------------------------------- f = (twr*twr_step)-(twi*twi_step); // Or, fast compute next twiddle twi = (twi*twr_step)+(twr*twi_step); // using complex vector multiply twr = f; // ** advance 2*pi/N radians ** f = 1.5 - 0.5*((twr*twr)+(twi*twi)); // correct amplitude errors if needed twi*=f; // (usually not needed) twr*=f; } COEF[t]= c * inv_FIR_size; } } #endif //----------------------------------- // Embed a DSK3DW graphpoint to view the impulse response volatile ETYPES GRAPHTYPE = NT_FLOAT; volatile int GRAPHLENGTH = FIR_size; volatile int GRAPHPOS = (int)COEF; asm(" .global GRAPHPOINT_0"); //----------------------------------- // FIR calculations are embedded into th for(;;) loop in main but can // also be called during the serial port interrupt c_int05(). This is // essentialy a control loop that produces output for the ISR that in // turn simply moves data to and from the serial port. //----------------------------------- #define Begin_Debug_ISR asm(" push IE "); asm(" and 0C4h,IE "); #define End_Debug_ISR asm(" pop IE "); void c_int05(void) { float f; Begin_Debug_ISR // Enable only debug interrupt if(ch) f = LEFT_OUT; else f = RIGHT_OUT; *S0_xdata = f; f = (*S0_rdata >>12); // +/- 2^19 range if(ch) LEFT_IN = f; else RIGHT_IN = f; End_Debug_ISR // restore all interrupts } void main(void) { GRAPHLENGTH = FIR_size; sdft_gen_coef(); asm("GRAPHPOINT_0:"); InitPeripherals(); asm(" ldi 0D4h,IE "); /* Enable XINT/INT2 */ for(;;) { asm(" idle"); // wait for sample ch = *LP_STATUS & 0x10; // ch is used elsewhere in program if(ch) LEFT_OUT = CTRLA_VAL[0]*FIR_L(LEFT_IN) + CTRLA_VAL[1]*FastRand(); else RIGHT_OUT = CTRLA_VAL[0]*FIR_R(RIGHT_IN) + CTRLA_VAL[1]*FastRand(); } } /*=================================================================== Install the XINT/RINT ISR branch vectors ====================================================================*/ asm(" .sect \"BRTBL\""); /* secondary branch table */ asm(" br _c_int05 "); /* XINT0 */ asm(" br _c_int05 "); /* RINT0 */ asm(" .text "); ===================================== At 01:10 PM 10/24/03 -0600, you wrote: Hello! I would like to know, if some one have some information about the company named MultiDSP, I found the reference in the DSP Applications using C a book from R. Chassaing, I am also found their web page (http://users.aol.com/multidsp/), but their email isnt work and either their phone. Well my interest in it, is because they have a digital filter design software (like FIR), and I am working in a speech project with the TMS320C31 DSK. Or may be some one knows some tool to do it the same. Thank you in Advance. Mart Lez CIATEC. Omega 201, Fracc. Delta C.P. 37545 Le, Gto. Mico www.ciatec.mx <http://www.ciatec.mx/> +-----------+ |Keith Larson | |Member Group Technical Staff | |Texas Instruments Incorporated | | | | 281-274-3288 | | | | www.micro.ti.com/~klarson | |-----------+ | TMS320C3x/C4x/VC33 Applications | | | | TMS320VC33 | | The lowest cost and lowest power 500 w/Mflop | | floating point DSP on the planet! | +-----------+ |