home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Multimed / Multimed.zip / rxwavsrc.zip / RxWavAnalyze.c < prev    next >
C/C++ Source or Header  |  2000-03-06  |  11KB  |  429 lines
  1. /*  RxWav
  2.    Copyright (C) 1999  Giorgio Vicario
  3.  
  4.    This program is free software; you can redistribute it and/or modify
  5.    it under the terms of the GNU General Public License as published by
  6.    the Free Software Foundation; either version 2 of the License, or
  7.    (at your option) any later version.
  8.  
  9.    This program is distributed in the hope that it will be useful,
  10.    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11.    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12.    GNU General Public License for more details.
  13.  
  14.    You should have received a copy of the GNU General Public License
  15.    along with this program; if not, write to the Free Software
  16.    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA     */
  17.  
  18. #define INCL_REXXSAA
  19. #include <os2emx.h>
  20. #include <stdlib.h>
  21. #include <stdio.h>
  22. #include <string.h>
  23. #include <regexp.h>
  24. #include <math.h>
  25. #include <float.h>
  26. #include "RxWav.h"
  27.  
  28. int *BitReversed;
  29.  
  30. typedef double t_fft;
  31.  
  32. /***********************************************************************
  33. Analyze
  34. Vuole il puntatore ad una traccia ed un parametro ('OTTAVE', 'TERZI', 'LINEARE')
  35. indicante il tipo di analisi (ottava, terzi di ottava o lineare) e
  36. restituisce uno stem con i valori richiesti
  37. ***********************************************************************/
  38. ULONG
  39. WavAnalyze (PCSZ name, LONG argc, const RXSTRING * argv,
  40.         PCSZ queuename, PRXSTRING retstr)
  41. {
  42.   PSHORT pCh = NULL;
  43.   ULONG nCamp, TipoOut, i, j, cband;
  44.   INT PuntiFFT;
  45.   SHORT c;
  46.   t_fft *data, *window, *output, re, im, outbuf, elementi;
  47.   double banda, k;
  48.   APIRET rc;
  49.   RXSTEMDATA spettro;
  50.  
  51.   if (argc != 3)
  52.     {
  53.       SendMsg (FUNC_ANALYZE, ERR_NUMERO_PARAMETRI);
  54.       return INVALID_ROUTINE;
  55.     }
  56.  
  57.   if (!sscanf (argv[0].strptr, "%d", &pCh))
  58.     {
  59.       SendMsg (FUNC_ANALYZE, ERR_PUNTATORE_ERRATO);
  60.       return INVALID_ROUTINE;
  61.     }
  62.  
  63.   TipoOut = BANDA_OTTAVE;
  64.   if (!strncmp (argv[1].strptr, "OTT", 3))
  65.     TipoOut = BANDA_OTTAVE;
  66.   if (!strncmp (argv[1].strptr, "TER", 3))
  67.     TipoOut = BANDA_TERZE;
  68.   if (!strncmp (argv[1].strptr, "CRO", 3))
  69.     TipoOut = BANDA_CROMATICA;
  70.   if (!strncmp (argv[1].strptr, "LOG", 3))
  71.     TipoOut = BANDA_LOGARITMICA;
  72.   if (!strncmp (argv[1].strptr, "LIN", 3))
  73.     TipoOut = BANDA_LINEARE;
  74.  
  75.   spettro.count = 0;
  76.   strcpy (spettro.varname, argv[2].strptr);
  77.   spettro.stemlen = argv[2].strlength;
  78.   strupr (spettro.varname);
  79.  
  80.   if (spettro.varname[spettro.stemlen - 1] != '.')
  81.     spettro.varname[spettro.stemlen++] = '.';
  82.  
  83.   switch (TipoOut)
  84.     {
  85.     case BANDA_LINEARE:
  86.       PuntiFFT = PUNTI_FFT_LINEARE;
  87.       banda = 1;
  88.       break;
  89.     case BANDA_OTTAVE:
  90.       PuntiFFT = PUNTI_FFT_OTTAVE;
  91.       banda = 20;
  92.       break;
  93.     case BANDA_TERZE:
  94.       PuntiFFT = PUNTI_FFT_TERZE;
  95.       banda = 20;
  96.       break;
  97.     case BANDA_CROMATICA:
  98.       PuntiFFT = PUNTI_FFT_CROMATICA;
  99.       banda = 20;
  100.       break;
  101.     case BANDA_LOGARITMICA:
  102.       PuntiFFT = PUNTI_FFT_LOGARITMICA;
  103.       banda = 20;
  104.       break;
  105.     default:
  106.       SendMsg (FUNC_ANALYZE, ERR_BANDA_ANALISI);
  107.       sprintf (retstr->strptr, "%i", ERR_BANDA_ANALISI);
  108.       retstr->strlength = strlen (retstr->strptr);
  109.       return INVALID_ROUTINE;
  110.       break;
  111.     }
  112.  
  113.   pCh = AllineaCh (pCh, (ULONG) PuntiFFT, FUNC_ANALYZE);
  114.   if (!pCh)
  115.     return INVALID_ROUTINE;
  116.  
  117.   InitFFT (PuntiFFT);
  118.  
  119.   data = (t_fft *) malloc (PuntiFFT * sizeof (t_fft));
  120.   window = (t_fft *) malloc (PuntiFFT * sizeof (t_fft));
  121.   output = (t_fft *) malloc (PuntiFFT * sizeof (t_fft));
  122.   if ((window == NULL) || (data == NULL) || (output == NULL))
  123.     {
  124.       SendMsg (FUNC_ANALYZE, ERR_ALLOCMEM);
  125.       return INVALID_ROUTINE;
  126.     }
  127.  
  128. /* windowing ???
  129.    for(i=0;i<PuntiFFT;i++)
  130.    window[i]=0.42-0.5*cos(2*M_PI*i/PuntiFFT)+0.08*cos(4*M_PI*i/PuntiFFT);
  131.  
  132.    for(i=0;i<PuntiFFT;i++)
  133.    data[i]=(t_fft)(*pCh++) * window[i];
  134.  */
  135.   for (i = 0; i < PuntiFFT; i++)
  136.     data[i] = (t_fft) (*pCh++);
  137.  
  138.   RealFFT (data);
  139.  
  140.   k = (double) FreqCamp / (double) PuntiFFT;
  141.   cband = 1;
  142.  
  143.   for (i = 0; i < (PuntiFFT / 2); i++)
  144.     {
  145.       re = data[BitReversed[i]];
  146.       im = data[BitReversed[i] + 1];
  147.       output[i] = sqrt (re * re + im * im);
  148.  
  149.       /*
  150.          if(output[i]>1) printf("freq %f pow---> %f\n", i*k, output[i]);
  151.        */
  152.  
  153.       switch (TipoOut)
  154.     {
  155.     case BANDA_LINEARE:
  156.       spettro.count++;
  157.       if (impostaOut (spettro, pow (output[i], 2), i))
  158.         {
  159.           sprintf (retstr->strptr, "%i", ERR_REXXPOOL);
  160.           retstr->strlength = strlen (retstr->strptr);
  161.           return INVALID_ROUTINE;
  162.         }
  163.       break;
  164.     case BANDA_OTTAVE:
  165.       if ((i * k) < banda)
  166.         {
  167.           cband++;
  168.           outbuf = outbuf + output[i];
  169.         }
  170.       else
  171.         {
  172.           spettro.count++;
  173.           if (impostaOut (spettro, (pow (outbuf, 2) / cband), banda))
  174.         {
  175.           sprintf (retstr->strptr, "%i", ERR_REXXPOOL);
  176.           retstr->strlength = strlen (retstr->strptr);
  177.           return INVALID_ROUTINE;
  178.         }
  179.           outbuf = output[i];
  180.           cband = 1;
  181.           banda = banda * 2;
  182.         }
  183.       break;
  184.     case BANDA_TERZE:
  185.       if ((i * k) < banda)
  186.         {
  187.           cband++;
  188.           outbuf = outbuf + output[i];
  189.         }
  190.       else
  191.         {
  192.           spettro.count++;
  193.           if (impostaOut (spettro, (pow (outbuf, 2) / cband), banda))
  194.         {
  195.           sprintf (retstr->strptr, "%i", ERR_REXXPOOL);
  196.           retstr->strlength = strlen (retstr->strptr);
  197.           return INVALID_ROUTINE;
  198.         }
  199.           outbuf = output[i];
  200.           cband = 1;
  201.           banda = banda * pow ((double) 2, (double) 1 / (double) 3);
  202.         }
  203.       break;
  204.     case BANDA_CROMATICA:
  205.       if ((i * k) < banda)
  206.         {
  207.           cband++;
  208.           outbuf = outbuf + output[i];
  209.         }
  210.       else
  211.         {
  212.           spettro.count++;
  213.           if (impostaOut (spettro, (pow (outbuf, 2) / cband), banda))
  214.         {
  215.           sprintf (retstr->strptr, "%i", ERR_REXXPOOL);
  216.           retstr->strlength = strlen (retstr->strptr);
  217.           return INVALID_ROUTINE;
  218.         }
  219.           outbuf = output[i];
  220.           cband = 1;
  221.           banda = banda * pow ((double) 2, (double) 1 / (double) 12);
  222.         }
  223.       break;
  224.     case BANDA_LOGARITMICA:
  225.       if ((i * k) < banda)
  226.         {
  227.           cband++;
  228.           outbuf = outbuf + output[i];
  229.         }
  230.       else
  231.         {
  232.           spettro.count++;
  233.           if (impostaOut (spettro, (pow (outbuf, 2) / cband), banda))
  234.         {
  235.           sprintf (retstr->strptr, "%i", ERR_REXXPOOL);
  236.           retstr->strlength = strlen (retstr->strptr);
  237.           return INVALID_ROUTINE;
  238.         }
  239.           outbuf = output[i];
  240.           cband = 1;
  241.           banda = banda * pow ((double) 2, (double) 1 / (double) 32);
  242.         }
  243.       break;
  244.     default:
  245.       SendMsg (FUNC_ANALYZE, ERR_BANDA_ANALISI);
  246.       sprintf (retstr->strptr, "%i", ERR_BANDA_ANALISI);
  247.       retstr->strlength = strlen (retstr->strptr);
  248.       return INVALID_ROUTINE;
  249.       break;
  250.     }
  251.     }
  252.  
  253.   elementi = spettro.count;
  254.   spettro.count = 0;
  255.  
  256.   if (impostaOut (spettro, elementi, (double) 0))
  257.     {
  258.       sprintf (retstr->strptr, "%i", ERR_REXXPOOL);
  259.       retstr->strlength = strlen (retstr->strptr);
  260.       return INVALID_ROUTINE;
  261.     }
  262.  
  263.   EndFFT ();
  264.   free (data);
  265.   free (window);
  266.   free (output);
  267.  
  268.   sprintf (retstr->strptr, "%f", ERR_OK);
  269.   retstr->strlength = strlen (retstr->strptr);
  270.   return VALID_ROUTINE;
  271. }
  272.  
  273.  
  274.  
  275. /***********************************************************************
  276. Inizializzazione e terminazione FFT (allocazione memoria...)
  277. ***********************************************************************/
  278. int *BitReversed;
  279. t_fft *SinTable;
  280. int punti = 0;
  281.  
  282. void
  283. InitFFT (int fftlen)
  284. {
  285.   int i;
  286.   int temp;
  287.   int mask;
  288.  
  289.   punti = fftlen / 2;
  290.  
  291.   if ((SinTable = (t_fft *) malloc (2 * punti * sizeof (t_fft))) == NULL)
  292.     {
  293.       puts ("Error allocating memory for Sine table.");
  294.       exit (1);
  295.     }
  296.  
  297.   if ((BitReversed = (int *) malloc (punti * sizeof (int))) == NULL)
  298.     {
  299.       puts ("Error allocating memory for BitReversed.");
  300.       exit (1);
  301.     }
  302.  
  303.   for (i = 0; i < punti; i++)
  304.     {
  305.       temp = 0;
  306.       for (mask = punti / 2; mask > 0; mask >>= 1)
  307.     temp = (temp >> 1) + (i & mask ? punti : 0);
  308.       BitReversed[i] = temp;
  309.     }
  310.  
  311.   for (i = 0; i < punti; i++)
  312.     {
  313.       SinTable[BitReversed[i]] = -sin (2 * M_PI * i / (2 * punti));
  314.       SinTable[BitReversed[i] + 1] = -cos (2 * M_PI * i / (2 * punti));
  315.     }
  316. }
  317.  
  318. void
  319. EndFFT ()
  320. {
  321.   free (BitReversed);
  322.   free (SinTable);
  323.   punti = 0;
  324. }
  325.  
  326.  
  327. RXSTEMDATA spettro;
  328.  
  329. int
  330. impostaOut (RXSTEMDATA spettro, t_fft valore, t_fft banda)
  331. {
  332.   /*
  333.      if (valore < 1)
  334.      valore = 0;
  335.    */
  336.   sprintf (spettro.varname + spettro.stemlen, "%d", spettro.count);
  337.   if (banda > 0)
  338.     sprintf (spettro.ibuf, "%f", sqrt (valore));
  339.   else
  340.     sprintf (spettro.ibuf, "%f", valore);
  341.  
  342.   spettro.shvb.shvnext = NULL;
  343.   spettro.shvb.shvname.strptr = spettro.varname;
  344.   spettro.shvb.shvname.strlength = strlen (spettro.varname);
  345.   spettro.shvb.shvvalue.strptr = spettro.ibuf;
  346.   spettro.shvb.shvvalue.strlength = strlen (spettro.ibuf);
  347.   spettro.shvb.shvcode = RXSHV_SET;
  348.   spettro.shvb.shvret = 0;
  349.   if (RexxVariablePool (&spettro.shvb) == RXSHV_BADN)
  350.     {
  351.       SendMsg (FUNC_ANALYZE, ERR_REXXPOOL);
  352.       printf ("       RexxVariablePool rc:%i\n", spettro.shvb.shvret);
  353.       return spettro.shvb.shvret;
  354.     }
  355.   else
  356.     return 0;
  357. }
  358.  
  359.  
  360.  
  361.  
  362. /***********************************************************************
  363. Esecuzione FFT reale
  364. ***********************************************************************/
  365. t_fft *A, *B;
  366. t_fft *sptr;
  367. t_fft *endptr1, *endptr2;
  368. int *br1, *br2;
  369. t_fft HRpiu, HRmeno, HIpiu, HImeno;
  370.  
  371. void
  372. RealFFT (t_fft * buffer)
  373. {
  374.   int ButterfliesPerGroup = punti / 2;
  375.  
  376.   endptr1 = buffer + punti * 2;
  377.  
  378.  
  379.   while (ButterfliesPerGroup > 0)
  380.     {
  381.       A = buffer;
  382.       B = buffer + ButterfliesPerGroup * 2;
  383.       sptr = SinTable;
  384.  
  385.       while (A < endptr1)
  386.     {
  387.       register t_fft sin = *sptr;
  388.       register t_fft cos = *(sptr + 1);
  389.       endptr2 = B;
  390.       while (A < endptr2)
  391.         {
  392.           t_fft v1 = *B * cos + *(B + 1) * sin;
  393.           t_fft v2 = *B * sin - *(B + 1) * cos;
  394.           *B = (*A + v1) * 0.5;
  395.           *(A++) = *(B++) - v1;
  396.           *B = (*A - v2) * 0.5;
  397.           *(A++) = *(B++) + v2;
  398.         }
  399.       A = B;
  400.       B += ButterfliesPerGroup * 2;
  401.       sptr += 2;
  402.     }
  403.       ButterfliesPerGroup >>= 1;
  404.     }
  405.   br1 = BitReversed + 1;
  406.   br2 = BitReversed + punti - 1;
  407.  
  408.   while (br1 <= br2)
  409.     {
  410.       register t_fft temp1;
  411.       register t_fft temp2;
  412.       t_fft sin = SinTable[*br1];
  413.       t_fft cos = SinTable[*br1 + 1];
  414.       A = buffer + *br1;
  415.       B = buffer + *br2;
  416.       HRpiu = (HRmeno = *A - *B) + (*B * 2);
  417.       HIpiu = (HImeno = *(A + 1) - *(B + 1)) + (*(B + 1) * 2);
  418.       temp1 = (sin * HRmeno - cos * HIpiu);
  419.       temp2 = (cos * HRmeno + sin * HIpiu);
  420.       *B = (*A = (HRpiu + temp1) * 0.5) - temp1;
  421.       *(B + 1) = (*(A + 1) = (HImeno + temp2) * 0.5) - HImeno;
  422.  
  423.       br1++;
  424.       br2--;
  425.     }
  426.   buffer[0] += buffer[1];
  427.   buffer[1] = 0;
  428. }
  429.