/* rockwell.c: */ /*** File: d2.c *** Author: Luigi Rizzo (luigi@iet.unipi.it) *** Date: 28-jul-95 *** *** This program converts from Rockwell 2-bit ADPCM format *** to Ulaw, ascii, signed short, SB (unsigned byte) format. *** *** Input from stdin, output to stdout. The output format *** can be chosen with argv[1] as follows: *** default signed short *** a ascii (useful for debugging) *** u Ulaw (without header) *** S SoundBlaster (unsigned byte) *** *** The program has several limitations, and is only intended to *** help those willing to use Rockwell-based voice modems. *** The decompression/compression code has been derived from *** the sources publicly released from Rockwell itself in 1995 *** The Ulaw conversion routines come from "sox" *** *** BUGS/LIMITATIONS *** + only deals decompression of 2-bit ADPCM format (compression *** should work, modulo some minor bugs, but I haven't tried it); *** + no sampling rate conversion. Use sox for this function. *** + the output quality is not very good. The modem sounds better, *** while running the conversion utility from Rockwell produces *** similar (not very good) results. I am wondering why is this! *** + the program uses FP instead of integer aritmetics. This saved *** me some work in the translation from ASM to C, and *should* *** allow a wider range in computations. However, it slows down the *** program. ***/ #include #include #define WORD short #define FWORD double #define SCALE_FACT 32768.0 /* Maximum Nu for the adaptive quantizer. */ #define QDLMX4 (WORD)0x2ED5 /* I)QDLMX = 1.11V */ #define QDLMX3 (WORD)0x3B7A /* I)QDLMX = 1.41V */ #define QDLMX2 (WORD)0x54C4 /* I)QDLMX = 2.01V */ #define F_QDLMX4 (11989.0/SCALE_FACT) #define F_QDLMX3 (15226.0/SCALE_FACT) /* #define F_QDLMX2 (21700.0/SCALE_FACT) /* original... */ #define F_QDLMX2 (21700.0 / SCALE_FACT ) /* Minimum Nu for the adaptive quantizer. */ #define QDLMN (WORD)0x001F /* QDLMN = IQDLMN = 2.87 mV. */ #define F_QDLMN (31.0 / SCALE_FACT ) /* Emphasis and De-emphasis filter coefficients. */ #define DEMPCF (WORD)0x3333 /* 0.4 */ #define PDEMPCF (WORD)0x3333 /* 0.4 */ #define F_DEMPCF (13107.0 / SCALE_FACT) #define F_PDEMPCF (13107.0 / SCALE_FACT) /* Total delay line length for the pole and zero delay lines. */ #define QORDER 8 /* ; Design Notes: Relation of QDataIndex to position in Qdata Buffer. ; ; Delay data is not rotated by shuffling; the rotation is done by ; the QDataIndex. This variable always points to the data to be ; multiplied by the coefficient a1. ; ; The coefficients, a1..a2 and b1..b6, stay in the same relative ; position in the coefficient array. Updates to these values are ; done in place. ; ; Illustration belows shows the value of QDataIndex and the Delay ; Line data in relation to the coefficient array. ; ; Position ; in Qdata Start 2nd 3rd 4th ; ----------------------------------------- ; 0 a1 Y(n-1) Y(n-2) Q(n-1) Q(n-2) ; 1 a2 Y(n-2) Q(n-1) Q(n-2) Q(n-3) ; 2 b1 Q(n-1) Q(n-2) Q(n-3) Q(n-4) ; 3 b2 Q(n-2) Q(n-3) Q(n-4) Q(n-5) ; 4 b3 Q(n-3) Q(n-4) Q(n-5) Q(n-6) ; 5 b4 Q(n-4) Q(n-5) Q(n-6) Y(n-1) ; 6 b5 Q(n-5) Q(n-6) Y(n-1) Y(n-2) ; 7 b6 Q(n-6) Y(n-1) Y(n-2) Q(n-1) ; ----------------------------------------- ; QDataIndex 0 7 6 5 ; ----------------------------------------- ; */ /* Coefficient Table for the pole and zero linear predictor. */ FWORD pzTable[] = { 0, 0, 0, 0, 0, 0, 0, 0 }; /* a1, a2, b1, b2, b3, b4, b5, b6 */ WORD QDataIndex= 0; /* Delay line pointer to the coefficient a1. */ FWORD Qdata[QORDER]= {0 }; /* Delay line. */ FWORD LastNu= 0; /* Last Nu value. Used in DecomOne, ComOne */ FWORD Dempz= 0; /* De-emphasis filter delay line (one element). */ /* Multiplier table to calculate new Nu for 2 BPS. */ WORD ML2bps[]= { 0x3333, /* ;ML0 = 1.6/4 */ 0x199A, /* ;ML1 = 0.8/4 */ 0x199A, /* ;ML2 = 0.8/4 */ 0x3333 /* ;ML3 = 1.6/4 */ }; FWORD F_ML2bps[]= { /* 1.6, .8, .8, 1.6 /* */ 1.6, .8, .8, 1.6 /* */ }; /* Multiplier table for 2 BPS to calculate inverse ; quantizer output. This number, index by the code ; word, times Nu is the inverse quantizer output. */ WORD Zeta2bps[]= { 0xCFAE, /* ;Zeta0 = -1.510/4 */ 0xF183, /* ;Zeta1 = -0.4528/4 */ 0x0E7D, /* ;Zeta2 = 0.4528/4 */ 0x3052 /* ;Zeta3 = 1.510/4 */ }; FWORD F_Zeta2bps[]= { -1.510, -0.4528, 0.4528, 1.510 /* */ }; /* Array index by BPS for updating QDelayMX */ FWORD QDelayTable[]= { 0, 0, F_QDLMX2, F_QDLMX3, F_QDLMX4 }; FWORD QDelayMX= F_QDLMX2 ; /* Maximum limit for Nu. */ long total_bytes=0; WORD fclip(FWORD x) { extern double rint(double); x *= SCALE_FACT; /* assumes SCALE_FACT-1 < 32767 */ if (x>=SCALE_FACT) return SCALE_FACT-1; else if (x<= -SCALE_FACT) return -SCALE_FACT; else return (short)(rint(x)); } FWORD clip(FWORD x, char *s) { if (x >=1.0) { #if 0 fprintf(stderr,"(%ld)%s: clip %lf to 1.0\n",total_bytes,s,x); #endif x= 1.0; } else if (x < -1.0) { #if 0 fprintf(stderr,"(%ld)%s: clip %lf to 1.0\n",total_bytes,s,x); #endif x= -1.0; } return x; } /*** Routines that both the decompressor and the compressor use. ***/ /*--------------------------------------------------------------------- ; pzPred ; Linear predictor coefficient update routine. Local to this module. ; a_i(n) = .98 a_i(n-1) + .012 sgn( q(n) y(n-i) ) ; b_i(n) = .98 b_i + .006 sgn( q(n) q(n-i) ) ; ; Inputs: Q(n), i.e. NewQdata ;---------------------------------------------------------------------*/ void pzPred(FWORD q_n) { FWORD a, b, c; int i; #define PNT98 .98 #define PNT012 .012 /* .012 */ #define PNT006 .006 /* .006 */ c= PNT012; for (i=0; i<8; i++) { if (i>=2) c=PNT006; a= PNT98 * pzTable[i]; b= q_n * Qdata[ (QDataIndex+i) & 7]; /* y(n-i-1) */ if (b<0) pzTable[i] = a - c; else pzTable[i] = a + c; pzTable[i]= clip(pzTable[i],"pzPred"); } } /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; Design Notes: Sum of Products Multiplications. ; ; lodsw ; a2(n) ; addm di,2 ; Y(n-2) ; imul Qdata[di] ; adjust dx,ax ; decimal point adjustment. ; add bx,ax ; add to 32 bit accumulator. ; adc cx,dx ; clip cx,bx ; check for overflow. ; ; Multiplications are 16-bit signed numbers producing a signed ; 32-bit result. The two operands are usually numbers less than ; one; this requires a 32-bit shift by the macro "adjust" to bring the ; decimal point in line. The 32-bit addition is two 16-bit additions ; with carry. The "clip" macro checks for overflow and limits the ; result of the addition to 0x7fffffff or 0x80000000 (for 32-bit ; results), or 0x7fff or 0x8000 (for 16-bit results). Note that ; the "clip" macro depends on the flags being set because of an ; addition; the 80x86 processor does not update these flags because ; of a move operation. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; XpzCalc ; ; Linear pole and zero predictor calculate. ; ; Inputs: initial value for accumulator ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ FWORD XpzCalc(FWORD acc) { int i; for (i=0; i<8; i++) { acc += pzTable[i]*Qdata[ (QDataIndex+i) & 7 ]; } acc=clip(acc,"XpzCalc"); return acc; } /*----------------------------------------------------------- ; ; Post Filter!! ; ------------------------------------------------------------*/ #ifdef POSTFILTER /* Adaptive post filter number 1 coefficients. */ WORD P8Table[]= { 0x6666, /* ;0.8 */ 0x51eb, /* ;0.8**2 */ 0x4189, /* ;0.8**3 */ 0x346d, /* ;0.8**4 */ 0x29f1, /* ;0.8**5 */ 0x218e /* ;0.8**6 */ }; /*** here I repeat the first two values to simplify the code ***/ FWORD F_P8Table[]= { .8, .64, .8, .64, .512, .4096, .32768, .262144 }; /* Adaptive post filter number 2 coefficients. */ WORD PM5Table[] = { 0xc000, /* ;-0.5 */ 0xe000, /* ;-(0.5**2) */ 0xf000, /* ;-(0.5**3) */ 0xf800, /* ;-(0.5**4) */ 0xfc00, /* ;-(0.5**5) */ 0xfe00 /* ;-(0.5**6) */ }; FWORD F_PM5Table[]= { -.5, -.25, -.5, -.25, -.125, -.0625, -.03125, -.015625 }; /* Delay line for the first stage of the post filter. */ FWORD QPdata[QORDER]= {0}; /* Delay line for the second stage of the post filter. */ FWORD QPPdata[QORDER] = {0 }; /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; App1Calc ; ; Adaptive post filter number 1. Local to this module. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ FWORD App1Calc(FWORD y) { WORD i, di=QDataIndex; FWORD a1, NewAppData=y; for (i=0;i<8; i++) { /* 0-1: a1(n) -> make a1'(n), 2-7: bi(n) --> make bi'(n) */ a1= pzTable[i] * F_P8Table[i]; /* 0-1: Y''(n-1), 2-7: QP(n-i) */ a1 *= QPdata[(di + i) & 7]; y += a1; } y=clip(y, "App1Calc"); /* Update delay line at the a1'(n) table entry position. */ QPdata[di & 7]= y; /* drop b6, now a1 Qdata. */ QPdata[(di+2) & 7]= NewAppData; /* drop a2, now b1 Qdata. */ return y; } /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; App2Calc ; ; Adaptive post filter number 2. Local to this module. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ FWORD App2Calc(FWORD y) { WORD i, di=QDataIndex; FWORD a1, NewAppData=y; for (i=0;i<8; i++) { /* 0-1: a1(n) -> make a1''(n), 2-7: bi(n) --> make bi''(n) */ a1= pzTable[i] * F_PM5Table[i]; /* 0-1: Y''(n-1), 2-7: QP(n-i) */ a1 *= QPPdata[(di + i) & 7]; y += a1; } y=clip(y, "App2Calc"); /* Update delay line at the a1''(n) table entry position. */ QPPdata[di & 7]= y; /* drop b6, now a1 Qdata. */ QPPdata[(di+2) & 7]= NewAppData; /* drop a2, now b1 Qdata. */ return y; } #endif /********** DECOMPRESSION ***************/ /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; RVDecomReset ; Reset local variables for ADPCM decoder. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ void RVDecomReset() { int i; QDataIndex=0; /* clear pointer to a1 */ for (i=0; i QDelayMX) LastNu=QDelayMX; NewQdata= OldNu * zeta; y=XpzCalc(NewQdata); #ifdef POSTFILTER /* Adaptive post filtering. The input to these routines ; are in CX. */ y=App1Calc(y); y=App2Calc(y); #endif pzPred( NewQdata ); /* Update predictor filter coefficients. */ /* Update delay line at the a1(n) table entry position. */ QDataIndex = (QDataIndex + 7) & 7; Qdata[QDataIndex]=y; /* drop b6, now a1 Qdata */ /* Update delay line at the b1(n) table entry position. */ Qdata[(QDataIndex +2) &7]=NewQdata; /* drop a2, now b1 Qdata */ /* Use a de-emphasis filter on Y(n) to remove the effects ; from the emphasis filter used during compression. */ Dempz=clip(F_DEMPCF * Dempz + y,"RVDecomOne"); /* Return Y(n) + de-emphasis filter. */ return fclip(Dempz); } /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; RVDecom2bps ; ; Decode 1 codeword, qd. Store 8 results ; at the location pointed to YnPtr. Use the callback function to ; return the amount of silence to the calling routine if a silence ; codeword is detected. Public routine. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ WORD RVDecom2bps(WORD qd, WORD *YnPtr) { int i; #define SL2CW 0x13ec if (qd==SL2CW) { fprintf(stderr,"Silence codeword!\n"); } for (i=0; i<8; i++) { *YnPtr++ = RVDecomOne(F_ML2bps[ qd & 3 ], F_Zeta2bps[ qd & 3 ]); qd >>= 2; } } /******************** COMPRESSION *********************/ #define CD2X_LENGTH 1 /* 4 level adaptive quantizer 2 bits per sample input constants. */ WORD cd2x= 0x1F69; /* ;Nu3 = 0.9816/4 */ FWORD f_cd2x = .9816; WORD adc16=0; /* Current X(n). */ WORD adc16z1= 0; /* X(n-1) */ WORD NewXpz=0; /* Xp+z(n). */ /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; RVComInit ; ; Reset local variables for the ADPCM encoder. Public routine. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ void RVComInit(WORD bps) { /* The compressor uses part of the decompression routines, ; initialize its local variables too. */ RVDecomInit(bps); /* Initialize the compression routines local variables. */ adc16=adc16z1=0; NewXpz=0; } /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; RVComOne2bps ; ; Code a word into a 2 BPS codeword. Local to this module. ; ; Input: ; ; AX = X(n) ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ WORD RVComOne2bps(WORD x_n) { FWORD delta, OldNu, cx, bx, dx; FWORD myQdata; int index; adc16=x_n; /* Save X(n) for later. */ NewXpz=XpzCalc(0); /* Compute Xp+z(n). */ /* Use the emphasis filter on Xp+z(n). * and compute the difference wrt the predictor */ delta= (adc16 - F_PDEMPCF * adc16z1) - NewXpz; bx= fabs(delta); adc16z1 = adc16; /*** values for index are as follows *** [ 0 ]<-c>[ 1 ]<0>[ 2 ][ 3 ] ***/ if (bx < f_cd2x* LastNu) index=1; else index=2; /* simple for 2 bits... */ if (delta>=0) index = index + 1; else index = 2 - index; /* Update Nu. */ OldNu= LastNu; LastNu = LastNu * ML2bps[index]; if (LastNu <= QDLMN) LastNu=QDLMN; else if (LastNu > QDelayMX) LastNu=QDelayMX; /* Make a new inverse quantizer value. */ myQdata = OldNu * Zeta2bps[index]; /* Update predictor filter coefficients. */ pzPred(myQdata); /* use inverse quantizer output. */ /* Update delay line at the a1(n) table entry position. */ QDataIndex = (QDataIndex + 7) & 7; Qdata[QDataIndex] = NewXpz + myQdata; /* drop b6, now a1 Qdata. */ /* Update delay line at the b1(n) table entry position. */ Qdata[ (QDataIndex + 2) & 7]=myQdata; /* drop a2, now b1 Qdata. */ return index; } /*;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ; RVCom2bps ; ; Encode 8 words pointed to by XnPtr. Returns result. ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;*/ WORD RVCom2bps(WORD *XnPtr) { WORD acc, j; acc=0; for (j=0; j<16;j+=2) { acc += (RVComOne2bps( *XnPtr++) << j); } return acc; } /* ** This routine converts from linear to ulaw. ** ** Craig Reese: IDA/Supercomputing Research Center ** Joe Campbell: Department of Defense ** 29 September 1989 ** ** References: ** 1) CCITT Recommendation G.711 (very difficult to follow) ** 2) "A New Digital Technique for Implementation of Any ** Continuous PCM Companding Law," Villeret, Michel, ** et al. 1973 IEEE Int. Conf. on Communications, Vol 1, ** 1973, pg. 11.12-11.17 ** 3) MIL-STD-188-113,"Interoperability and Performance Standards ** for Analog-to_Digital Conversion Techniques," ** 17 February 1987 ** ** Input: Signed 16 bit linear sample ** Output: 8 bit ulaw sample */ #define ZEROTRAP /* turn on the trap as per the MIL-STD */ #define BIAS 0x84 /* define the add-in bias for 16 bit samples */ #define CLIP 32635 unsigned char st_linear_to_ulaw( WORD sample ) { static int exp_lut[256] = {0,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3, 4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4, 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7}; int sign, exponent, mantissa; unsigned char ulawbyte; /* Get the sample into sign-magnitude. */ sign = (sample >> 8) & 0x80; /* set aside the sign */ if ( sign != 0 ) sample = -sample; /* get magnitude */ if ( sample > CLIP ) sample = CLIP; /* clip the magnitude */ /* Convert from 16 bit linear to ulaw. */ sample = sample + BIAS; exponent = exp_lut[( sample >> 7 ) & 0xFF]; mantissa = ( sample >> ( exponent + 3 ) ) & 0x0F; ulawbyte = ~ ( sign | ( exponent << 4 ) | mantissa ); #ifdef ZEROTRAP if ( ulawbyte == 0 ) ulawbyte = 0x02; /* optional CCITT trap */ #endif return ulawbyte; } /* ** This routine converts from ulaw to 16 bit linear. ** ** Craig Reese: IDA/Supercomputing Research Center ** 29 September 1989 ** ** References: ** 1) CCITT Recommendation G.711 (very difficult to follow) ** 2) MIL-STD-188-113,"Interoperability and Performance Standards ** for Analog-to_Digital Conversion Techniques," ** 17 February 1987 ** ** Input: 8 bit ulaw sample ** Output: signed 16 bit linear sample */ WORD st_ulaw_to_linear_slow(unsigned char ulawbyte ) { static int exp_lut[8] = { 0, 132, 396, 924, 1980, 4092, 8316, 16764 }; int sign, exponent, mantissa, sample; ulawbyte = ~ ulawbyte; sign = ( ulawbyte & 0x80 ); exponent = ( ulawbyte >> 4 ) & 0x07; mantissa = ulawbyte & 0x0F; sample = exp_lut[exponent] + ( mantissa << ( exponent + 3 ) ); if ( sign != 0 ) sample = -sample; return sample; } void MakeSB(WORD *sp, char *dp, int count) { int i; for (i=count; i > 0; --i) *dp++ = (char)(((*sp++ >> 8) + 0x80) & 0xff); } void MakeUlaw(WORD *sp, unsigned char *dp, int count) { int i; for (i=count; i > 0; --i) *dp++ = st_linear_to_ulaw(*sp++); } main(int argc, char *argv[]) { WORD inbuf[1]; WORD outbuf[8]; int i, count=16; RVDecomInit(2); while (read(0, inbuf, 2) == 2) { RVDecom2bps(inbuf[0], outbuf); total_bytes += 8; if (argc>1) { switch (*argv[1]) { case 'a': for (i=0;i<8; i++) fprintf(stdout,"%d\n",outbuf[i]); count=0; /* trick */ break; case 'u': MakeUlaw(outbuf, (char *)outbuf, 8); count=8; break; case 'S': MakeSB(outbuf, (char *)outbuf, 8); count=8; break; default: count=16; } } if (count) write(1, outbuf, count ); } } /*** eof ***/