home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / magazine / drdobbs / 1991 / 05 / co_dsp.asc < prev    next >
Text File  |  1991-03-15  |  21KB  |  529 lines
  1. _ADDING THE POWER OF DSP TO YOUR APPLICATIONS_
  2. by Jim Bittman
  3.  
  4. [LISTING ONE]
  5. /*-------------------- DDJ.C-----------------------------------------
  6. |  PC-Resident program for FFT and for controlling DSP co-processor
  7. |  The Makefile for this program is:
  8. |   .c.obj: 
  9. |         cl /AL -c -DANSI $<
  10. |   ddj.exe: ddj.obj dspif.obj dspif.h four1.obj realft.obj
  11. |         cl /AL ddj.obj dspif.obj four1.obj realft.obj \
  12. |             -link dsputill.lib hqmcil.lib
  13. |  This program assumes the DspHq environment and requires additional
  14. |  support files (menu definition and function specification files)
  15. |  not listed in the Makefile. However, it should not be difficult to
  16. |  excerpt the relevant portions of this code for standalone execution.
  17. +------------------------------------------------------------------*/
  18.  
  19. #include <stdio.h>              /* include for NULL define         */
  20. #include <math.h>               /* include for math functions      */
  21. #include "hq_ci_.h"             /* include for dsphq interface     */
  22. #include "nr.h"                 /* include nr function prototypes  */
  23. #include "nrutil.h"             /* include nr utility prototypes   */
  24. #include "dsputil.h"            /* include CAC prototypes          */
  25. #include "dspif.h"              /* dsp interface function header   */
  26.  
  27. /*---------- define the menu parameter types ----------------------*/
  28. typedef unsigned long parm1_type;       /* Buffer 1                */
  29. typedef unsigned long parm2_type;       /* Buffer 3                */
  30. typedef menu_float    parm3_type;       /* Cosine Amplitude        */
  31. typedef menu_float    parm4_type;       /* Cosine DC-offset        */
  32. typedef menu_float    parm5_type;       /* Cosine Frequency        */
  33. typedef menu_float    parm6_type;       /* Cosine Sammple Rate     */
  34. /*------------ constant definitions -------------------------------*/
  35. #define SQR(a) ((a)*(a))
  36. #define PI               (float) 3.141592654
  37. #define BAD_FUNC_NUM     (int) -1       /* user error code */
  38. /*----- Function Constants For Dsp32 'C' & Assembly Code ------*/
  39. #define DSP32_X     1      /* converts IEEE ==> DSP numbers */
  40. #define IEEE32_X    2      /* converts DSP  ==> IEEE numbers */
  41. /*----- Function Constants For Dsp32 'C' Code ----------------*/
  42. #define GENCOS_C    3      /* generate test cosine */
  43. #define REALFT_C    4      /* performs realft, from numerical recipes */
  44. #define IREALFT_C   5      /* performs inv_realft, from numerical recipes */
  45. #define LOGMAG_C    6      /* calculates log-magnitude */
  46. #define RFFTA_C     7      /* calculates real fft using app-lib */
  47. /*----- Function Constants For Dsp32 Assembly Code ------------*/
  48. #define RFFT_S      3      /* performs rfft, from AT&T application library */
  49. #define MAG_S       4      /* calculates magnitude-squared */
  50. #define LOG_S       5      /* calculates log */
  51. /*----- Function Prototypes ----------------------------------*/
  52. int     ilog2      (unsigned int);
  53. void    log_mag    (float far * i1, float far * o1, long bs);
  54. void    scale_data (float far *output1, float scale_val, long np);
  55. void    synth_cos  (float far *data1, long np, float a, float d, float f, 
  56.                                                                      float s);
  57. main(int argc, char *argv[])
  58. {  long                 indx;           /* used for loop index */
  59.    int                  func_num;       /* for function number from dsphq */
  60.    long                 np;             /* for block size from dsphq */
  61.    float far *          input1;         /* array address */
  62.    float far *          output1;        /* array address */
  63.    long                 bs;             /* address of DSP blocksize var */
  64.    long                 flag;           /* address of DSP funcnum flag */
  65.    parm1_type           b1;             /* DSP Buffer #1 */
  66.    parm2_type           b2;             /* DSP Buffer #2 */
  67.    parm3_type           amp;            /* cosine amplitude */
  68.    parm4_type           dco;            /* cosine DC offset */
  69.    parm5_type           freq;           /* cosine frequency */
  70.    parm6_type           samprate;       /* cosine sample rate */
  71.    init_intfc(argc, argv);              /* init dsphq interface */
  72.    func_num = get_func_num();           /* get the function number */
  73.    np = get_block_size();               /* get the block size */
  74.    /* get menu parameters */
  75.    b1   = get_parm(1);                  /* DSP Buffer #1 */
  76.    b2   = get_parm(2);                  /* DSP Buffer #2 */
  77.    amp  = get_parm(3);                  /* cosine amplitude */
  78.    dco  = get_parm(4);                  /* cosine DC offset */
  79.    freq = get_parm(5);                  /* cosine frequency */
  80.    samprate = get_parm(6);              /* cosine sample rate */
  81.    /* get array addresses */
  82.    input1  = get_data_in_ptr(1);        /* base address of input #1 */
  83.    output1 = get_data_out_ptr(1);       /* base address of output #1 */
  84.    /* perform selected function */
  85.    switch (func_num)
  86.    {
  87.       case 1 :       /*--- Synthesize Cosine Using PC ---------------------*/
  88.          synth_cos(output1, np, amp, dco, freq, samprate);
  89.          break;
  90.       case 2 :      /*---- Numerical Recipes Forward Real FFT Using PC ----*/
  91.          output1--;                             /* NR funcs index at 1     */
  92.          realft(output1, np>>1, 1);             /* forward real fft        */
  93.          break;
  94.       case 3 :       /*---- Numerical Recipes Inverse Real FFT Using PC ----*/
  95.          output1--;                             /* NR funcs index at 1      */
  96.          realft(output1, np>>1, -1);            /* inverse real fft         */
  97.          output1++;                             /* reallign address         */
  98.          scale_data(output1,1.0/(np >> 1),np);  /* restore original ampl.   */
  99.          break;
  100.       case 4 :     /*------ Calculate LOG(10)-[MAGNITUDE] using PC --------*/
  101.          if (input1)
  102.             log_mag(input1, output1, np);       /* take logmag of input     */
  103.          else
  104.             log_mag(output1, output1, np);      /* perform logmag in-place  */
  105.          break;
  106.       case 5 :     /*--------- Synthesize Cosine Using DSP-32C -------------*/
  107.          init_dsp("ddj_32c.32c",&flag,&bs,np);  /* download dsp code & init */
  108.          dsp_dl_fp(get_addr("amp"),amp);        /* download floats          */
  109.          dsp_dl_fp(get_addr("dco"),dco);
  110.          dsp_dl_fp(get_addr("freq"),freq);
  111.          dsp_dl_fp(get_addr("samprate"),samprate);
  112.          set_dspbuf("o1", b1);                  /* set dsp buffer address   */
  113.          dsp_dl_int(flag,GENCOS_C);             /* invoke function on dsp   */
  114.          wait4dsp(flag);                        /* wait for dsp to finish   */
  115.          ul_float(output1,np,flag,b1);          /* upload results           */
  116.          break;
  117.       case 6 :  /*---- Numerical Recipes Forward Real FFT Using DSP-32C ----*/
  118.          init_dsp("ddj_32c.32c",&flag,&bs,np);  /* download dsp code & init */
  119.          set_dspbuf("o1", b1);                  /* set dsp buffer address   */
  120.          dl_float(input1,np,flag,b1);           /* download float array     */
  121.          dsp_dl_int(flag,REALFT_C);             /* execute "realft" on dsp  */
  122.          wait4dsp(flag);                        /* wait for dsp to finish   */
  123.          ul_float(output1,np,flag,b1);          /* upload results           */
  124.          break;
  125.       case 7 : /*---- Numerical Recipes Inverse Real FFT Using DSP-32C -----*/
  126.          init_dsp("ddj_32c.32c",&flag,&bs,np);  /* download dsp code & init */
  127.          dl_float(input1,np,flag,b1);           /* download float array     */
  128.          set_dspbuf("o1", b1);                  /* set dsp buffer address   */
  129.          dsp_dl_int(flag,IREALFT_C);         /* execute "inv_realft" on dsp */
  130.          wait4dsp(flag);                        /* wait for dsp to finish   */
  131.          ul_float(output1,np,flag,b1);          /* upload results           */
  132.          break;
  133.       case 8 :      /*----- Calculate LOG(10)-[MAGNITUDE] using DSP-32C ----*/
  134.          init_dsp("ddj_32c.32c",&flag,&bs,np);  /* download dsp code & init */
  135.          dl_float(input1,np,flag,b1);           /* download float array     */
  136.          set_dspbuf("i1", b1);                  /* set dsp buffer address   */
  137.          set_dspbuf("o1", b2);                  /* set dsp buffer address   */
  138.          dsp_dl_int(flag,LOGMAG_C);          /* execute "inv_realft" on dsp */
  139.          wait4dsp(flag);                        /* wait for dsp to finish   */
  140.          ul_float(output1,np,flag,b2);          /* upload results           */
  141.          break;
  142.       case 9 :     /*------- Forward Real FFT Using DSP-32C App-Lib --------*/
  143.          init_dsp("ddj_32s.32c",&flag,&bs,np);  /* download dsp code & init */
  144.          dl_float(input1,np,flag,b1);           /* download float array     */
  145.          set_dspbuf("o1", b1);                  /* set dsp buffer address   */
  146.          dsp_dl_int(get_addr("stages"),ilog2(np));
  147.                                                 /* download int             */
  148.          dsp_dl_int(flag,RFFT_S);               /* execute "rfft" on dsp    */
  149.          wait4dsp(flag);                        /* wait for dsp to finish   */
  150.          ul_float(output1,np,flag,b1);          /* upload results           */
  151.          break;
  152.       case 10:     /*-------- Download Data To DSP-32C ---------------------*/
  153.          init_dsp("ddj_32s.32c",&flag,&bs,np);  /* download dsp code & init */
  154.          dl_float(input1,np,flag,b1);       /* download data from pc to dsp */
  155.          break;
  156.       case 11:    /*--------- Upload Data From DSP-32C ---------------------*/
  157.          init_dsp("ddj_32s.32c",&flag,&bs,np);  /* download dsp code & init */
  158.          ul_float(output1,np,flag,b1);          /* upload results           */
  159.          break;
  160.       default : set_err_return(BAD_FUNC_NUM); break;
  161.    }
  162. }
  163. /*--- This function returns the integer part of the log base 2 of x. ---*/
  164. int ilog2(unsigned int x)
  165. {
  166.    return( x >> 1 ? 1 + ilog2(x >> 1) : 0);
  167. }
  168. /*--- This function scales np elements of data1[] by scale_val. ---*/
  169. void scale_data(float far * data1, float scale_val, long np)
  170. {  long i;
  171.    for (i = 0; i < np; i++)   {
  172.       data1[i] *= scale_val;
  173.    }
  174. }
  175. /*--- Function generates cosine data: data[i] = A cos((2 pi f/s) i) + d ---*/
  176. void synth_cos(float far * data1, long np, float a, float d, float f, float s)
  177. {  long i;
  178.    float theta, angle_step;
  179.    angle_step = 2.0 * PI * f / s;
  180.    theta = 0.0;
  181.    for (i = 0; i < np; i++) {
  182.       data1[i] = (a * cos(theta)) + d;
  183.       theta += angle_step;
  184.    }
  185. }
  186. /*--- log_mag ---*/
  187. void log_mag(float far * i1, float far * o1, long bs)
  188. {  long i;
  189.    long n;
  190.    n = bs >> 1;
  191.    o1[0] = log10(SQR(i1[0]));
  192.    for (i = 1; i < n; i++)   {
  193.       o1[i] = log10(SQR(i1[2*i]) + SQR(i1[2*i+1]));
  194.    }
  195.    for (i = n; i < bs; i++)   {
  196.       o1[i] = 0.0;
  197.    }
  198. }
  199.  
  200.  
  201. [LISTING TWO]
  202.  
  203. /*-----------------DDJ_MP87.C-------------------------------------------
  204. |  MathPak 87 FFT Function Group Execution Source File 
  205. |  The makefile is: ddj_mp87.exe: ddj_mp87.c
  206. |                      cl /AL ddj_mp87.c -link hqmcil.lib mpak87.lib
  207. +----------------------------------------------------------------------*/
  208. #include <hq_ci_.h>                     /* include function prototypes */
  209. #include <mpak87.h>                     /* include math pak header */
  210. #define BAD_BLOCK_SIZE        -1        /* user defined error codes */
  211. #define BAD_FUNC_NUM          -2
  212.  
  213. int fft_stages(long fft_size); /*function prototype*/
  214.  
  215. void main(int argc, char *argv[])
  216. {
  217.    int         m;       /* number of fft stages */
  218.    far_array_of_double o1;  /* data pointer */
  219.    init_intfc(argc, argv); /* MUST do this before other interface functions */
  220.    o1 = get_data_out_ptr(1);   /* get address of data */
  221.    if ((m = fft_stages(get_block_size())) == 0)
  222.       set_err_return(BAD_BLOCK_SIZE);  /* won't happen if .fnc file correct */
  223.    else
  224.    switch(get_func_num())   {
  225.       case 1: rvfft(o1, m); break;    /* real fft from MathPak library     */
  226.       case 2: irvfft(o1, m); break;   /* inverse real fft from MathPak     */
  227.       default:                           
  228.         set_err_return(BAD_FUNC_NUM); /* won't happen if .fnc file correct */
  229.       break;
  230.    }
  231. } 
  232. /* return the log(base 2) of the input, or 0 if input is not a power of 2 */
  233. int fft_stages(long fft_size)
  234. {  int rtn;
  235.    int sw_fft;
  236.    sw_fft = fft_size;
  237.    switch(sw_fft)   {
  238.       case    8 : rtn =  3; break;
  239.       case   16 : rtn =  4; break;
  240.       case   32 : rtn =  5; break;
  241.       case   64 : rtn =  6; break;
  242.       case  128 : rtn =  7; break;
  243.       case  256 : rtn =  8; break;
  244.       case  512 : rtn =  9; break;
  245.       case 1024 : rtn = 10; break;
  246.       case 2048 : rtn = 11; break;
  247.       case 4096 : rtn = 12; break;
  248.       case 8192 : rtn = 13; break;
  249.       default   : rtn =  0; break;
  250.    }
  251.    return(rtn);
  252. }
  253.  
  254.  
  255.  
  256. [LISTING THREE]
  257.  
  258. /*--------------------DDJ_32C.C--------------------------------------
  259. | This file is will run on the DSP add-in board after compilation
  260. | by the AT&T C compiler d3cc. The makefile is:
  261. |  .c.o: 
  262. |       d3cc -c $<
  263. |  .s.o: 
  264. |       d3as -W -Q $<
  265. |   ddj_32c.32c: ddj_32c.o startup.o memory.map four1.o realft.o
  266. |       d3cc -lm -lap -o ddj_32c.32c ddj_32c.o four1.o realft.o \
  267. |        -s startup.o -m memory.map
  268. +--------------------------------------------------------------------*/
  269. #include <stdio.h>
  270. #include <math.h>
  271. #include <nr.h>
  272.  
  273. #define PI 3.1415926
  274. #define SQR(a) ((a)*(a))
  275.  
  276. asm(".global i1, i2, o1, o2");
  277. asm(".global funcnum, bs");
  278. asm(".global amp, dco, freq, samprate");
  279.  
  280. short funcnum;
  281. short bs;
  282. float *i1, *i2, *o1, *o2;
  283. float amp, dco, freq, samprate;
  284. /*------------------------------------------------------------*/
  285. short fft_stages(fft_size)   
  286. short fft_size;
  287. {
  288.    short rtn;
  289.    switch(fft_size)
  290.    {
  291.       case   32 :   rtn =  5;    break;
  292.       case   64 :   rtn =  6;    break;
  293.       case  128 :   rtn =  7;    break;
  294.       case  256 :   rtn =  8;    break;
  295.       case  512 :   rtn =  9;    break;
  296.       case 1024 :   rtn = 10;    break;
  297.       case 2048 :   rtn = 11;    break;
  298.       case 4096 :   rtn = 12;    break;
  299.       case 8192 :   rtn = 13;    break;
  300.       default   :   rtn =  0;    break;
  301.    }
  302.    return(rtn);
  303. }
  304. /*------------------------------------------------------------*/
  305. main()
  306. {    register float scal;
  307.      short n;
  308.      float *data1, *data2, temp;
  309.      register short i,j;
  310.  
  311.      while (1) {
  312.        funcnum=0;
  313.        while (!(funcnum))
  314.        ;       /* Wait for PC to Download Function Number */
  315.        n = bs >> 1;
  316.        switch(funcnum) {
  317.           case 1: /*--------- Convert to DSP Format ------------*/
  318.                     dsp32(bs,o1); /*IEEE-->DSP */
  319.                     break;
  320.           case 2: /*--------- Convert to IEEE format -----------*/
  321.                     ieee32(bs,o1); /*DSP-->IEEE*/
  322.                     break;
  323.           case 3: /*--------- Synthesize Cosine ----------------*/
  324.                     scal = 2.0 * PI * freq / samprate;
  325.                     j = 0;
  326.                     data1 = o1;
  327.                     for (i=bs; i-- > 0; j++) {
  328.                       *data1++ = (amp * cos(scal * j)) + dco;
  329.                     }
  330.                     break;
  331.           case 4: /*------- Forward FFT ------------------------*/
  332.                     data1 = o1;
  333.                     data1--;
  334.                     realft(data1,n,1); /*from numerical recipes*/
  335.                     break;
  336.           case 5: /*------- Inverse FFT ------------------------*/
  337.                     data1 = o1;
  338.                     data1--;
  339.                     realft(data1,n,-1); /*from numerical recipes*/
  340.                     /* scale by 1/n to retain original amplitude*/
  341.                     data1 = o1;
  342.                     scal = 1.0 / n;
  343.                     for (i=bs; i-- > 0; data1++) {
  344.                       *data1 = *data1 * scal;
  345.                     }
  346.                     break;
  347.           case 6: /*----- Calc LOG-MAGNITUDE (data output from NR-realft)--*/
  348.                     o1[0] = log10(SQR(i1[0]));
  349.                     temp  = log10(SQR(i1[1]));
  350.                     for (i=1;i<n;i++) {
  351.                       o1[i]=log10(SQR(i1[2*i])+SQR(i1[2*i+1]));
  352.                     }
  353.                     o1[n] = temp;
  354.                     for (i=n+1; i<bs; i++) {
  355.                       o1[i] = 0.0;
  356.                     }
  357.                     break;
  358.           case 7: /*------ Forward FFT ----------------------------------*/
  359.                     data1 = o1;
  360.                     rffta(bs,fft_stages(bs),data1); /*uses AT&T app lib*/
  361.                     break;
  362.           default : break;
  363.        }
  364.      }
  365. }
  366.  
  367.  
  368. [LISTING FOUR]
  369.  
  370. /*------ file DDJ_32S.S-----------------------------------------------
  371. | Assembly language version of FFT test. The makefile is as follows:
  372. |   ddj_32s.32c: ddj_32s.s
  373. |        d3make -o ddj_32s.32c -M6 -Q -W ddj_32s.s
  374. +---------------------------------------------------------------------*/
  375. #include <dspregs.h>
  376. /*--------------------------------------------------------------------*/
  377. .global i1,o1,o2
  378. .global funcnum, bs, stages
  379. .global endofcode
  380. /*--------------------------------------- initialization -------------*/
  381.        r5 = 0x0003
  382.        pcw = r5
  383.        ioc  = 0x30CC0
  384. /*--------------------------------------- wait until funcnum != 0 ----*/
  385. begin:
  386.        r5 = *funcnum
  387.        nop
  388.        if (eq) goto begin
  389.        nop
  390. /*--------------- switch on funcnum, which is set to function number --*/
  391. /* func 1:  IEEE-->DSP */
  392.      r5 = r5 - 1
  393.      if (eq) goto do_dsp32
  394. /* func 2: DSP-->IEEE */
  395.      r5 = r5 - 1
  396.      if (eq) goto do_ieee32
  397. /* func 3: invoke rffta, from AT&T app library */
  398.      r5 = r5 - 1
  399.      if (eq) goto rffta
  400. /* func 4: calc magnitude squared */
  401.      r5 = r5 - 1
  402.      if (eq) goto do_mag
  403.      nop
  404. /* func 5: calc log */
  405.      r5 = r5 - 1
  406.      if (eq) goto do_log10
  407.      nop
  408. /* illegal function number */
  409.      goto finished
  410.      nop
  411. /*---------------------------------call to _rffta -----------------*/
  412. rffta:
  413.        r2e = *o1
  414.        r1  = *bs
  415.        r3  = *stages
  416.        *fft_b = r2e
  417.        *fft_n = r1e
  418.        *fft_m = r3e
  419.        call _rffta (r14)
  420.        nop
  421. fft_lv: int24    localV
  422. fft_n:  int24    0
  423. fft_m:  int24    0
  424. fft_b:  int24    0
  425. .align 4
  426.        goto finished
  427.        nop
  428. /*---------------------------------calc magnitude------------------*/
  429. do_mag:
  430.        r8   = *bs
  431.        r10e = *i1
  432.        r8   = r8 - 2
  433.        r11e = *o1
  434.        nop
  435.        a0 = *r10++                 /* DC value */
  436.        nop
  437.        a2 = *r10++                 /* Nyquist  */
  438.        *r11++ = a0 = a0*a0         /* save DC mag */
  439. magloop:
  440.        a0 = *r10++
  441.        nop
  442.        a1 = *r10++
  443.        a0 = a0*a0
  444.        nop
  445.        a1 = a1*a1
  446.        nop
  447.        nop
  448.        *r11++ = a1 = a0 + a1
  449.        if (r8-->=0) goto magloop
  450.        nop
  451.        *r11++ = a0 = a2*a2         /* save Nyquist mag */
  452.        goto finished
  453.        nop
  454. /*---------------------------------calc log10------------------*/
  455. do_log10:
  456.        r12e = *i1
  457.        r11e = *o1
  458.        r13  = *bs
  459.        r10e = in
  460.        r9e  = out
  461.        r13  = r13 - 2
  462. logloop:
  463.        *r10 = a3 = *r12++
  464.        call _log10 (r14)
  465.        nop
  466.        int24     localV
  467.        int24     in, out
  468. .align 4
  469.        *r11++ = a3 = *r9
  470.        if (r13-->=0) goto logloop
  471.        nop
  472.        goto finished
  473.        nop
  474. /*--------------------------------------call _ieee32 converter-----*/
  475. do_ieee32:
  476.        r1e = *o1
  477.        r2 = *bs
  478.        *o1_ieee32 = r1e
  479.        *bs_ieee32 = r2e
  480.        call _ieee32 (r14)
  481.        nop
  482.             int24       localV
  483. bs_ieee32:  int24       0
  484. o1_ieee32:  int24       0
  485. .align 4
  486.        goto finished
  487.        nop
  488. /*--------------------------------------call _dsp32 converter-----*/
  489. do_dsp32:
  490.        r1e = *o1
  491.        r2 = *bs
  492.        *o1_dsp32 = r1e
  493.        *bs_dsp32 = r2e
  494.        call _dsp32 (r14)
  495.        nop
  496. bs_dsp32:  int24       0
  497. o1_dsp32:  int24       0
  498. .align 4
  499.        goto finished
  500.        nop
  501. /*-------------------------------finished, set funcnum=0 -----------*/
  502. finished: 
  503.        r1=0
  504.        goto begin
  505.        *funcnum=r1
  506. bs:         int 256
  507. stages:     int 8
  508. funcnum:    int 0
  509. .align 4
  510. i1: int24 0x2000
  511. o1: int24 0x2000
  512. o2: int24 0x3000
  513.  
  514. .align 4
  515. localV:   2*float   0.0
  516. in:       float     0.0
  517. out:      float     0.0
  518. max:      float     0.0
  519. scalefac: float     0.0
  520. #include <_rffta.asm>
  521. #include <_log10.asm>
  522. #include <_ieee32.asm>
  523. #include <_dsp32.asm>
  524. /*------------------------------------- mark end of code------------*/
  525. endofcode: int 0xDEAD, 0xC0DE
  526.  
  527.  
  528.  
  529.