home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / demos / audio / amesh / Power.c++ < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  3.3 KB  |  147 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17.  
  18. #include <stdio.h>
  19. // #include "Sound.h"
  20. // #include "Audio.h"
  21. extern "C" {
  22. #include        <math.h>
  23. }
  24. typedef struct { float  r; float  i; } complex;
  25.  
  26. //  Object
  27. // To read in a sampled sound file 
  28. // compute an estimate of the power Spectrum of a Delta t
  29. // write out the pse for each time interval
  30. // pse should be normalized to 0 to 1.0
  31. // Question should we plot log(power) or power???
  32.  
  33.  
  34. void NRICfft( float *, unsigned int, int);
  35.  
  36.  
  37.  
  38. #define FILTERSIZE  16384
  39. #define FILTERSIZE2 32768
  40. #define LOGFILTERSIZE2 15
  41.  
  42. #define PI2 6.283185307179586476925287
  43. int nfreqs, nintervals;
  44.  
  45. #define MAXFREQS 256
  46. #define MAXINTERVALS 2000
  47. #define MAXSAMPS 8192
  48.  
  49. complex  tmpc[8192];
  50. void NRICfft(float *data,unsigned int nn,int isign);
  51. typedef signed short SAMPLE;
  52.  
  53.  
  54. float spectrum( SAMPLE * samples, int nsamples, float * ps, int nfreqs, float maxpwr )
  55.     {
  56.     int nfft;
  57.     int i;
  58.     SAMPLE * pb;
  59.     float pwr,nfftsqd;
  60.     complex * pc;
  61.  
  62.  
  63.     for(pb = samples, pc = tmpc, i = 0; i < nsamples ; i+=2)
  64.         { pc->r  = (float) *pb++;
  65.           pc->r  += (float) *pb++;
  66.           pc->i = 0.0; pc++;
  67.         } 
  68.     nsamples >>= 1;
  69.     nfft = 1;
  70.     while (nfft < nsamples) nfft <<= 1;
  71.     if (nfft > nsamples )   // zero fill
  72.        for(i = 0; i < nfft-nsamples ; i++)
  73.         { pc->r = 0.0; pc->i = 0.0 ; pc++; }
  74.     pc = tmpc;
  75.  
  76.     NRICfft( (float *) pc,nfft,1);
  77.      nfftsqd = float(nfft)*float(nfft);
  78.  
  79.     for(i=1; i<nfreqs; i++)
  80.         {
  81.         
  82.         pwr = tmpc[i].r * tmpc[i].r;
  83.         pwr += tmpc[i].i * tmpc[i].i;
  84.         pwr = fsqrt(pwr/nfftsqd);
  85.         pwr *= 1.0 + (( (float)(i)/(float)nfreqs)*.3);
  86.         ps[i-1] = pwr;
  87.         if (pwr > maxpwr) maxpwr = pwr;
  88.         }
  89.     return(maxpwr);
  90.     }
  91.  
  92.  
  93.  
  94.  
  95.  
  96. /* FFT from Numerical recipes in C page 411 */
  97. #define SWAP(a,b) tempr=(a); (a)=(b); (b) = tempr
  98.  
  99. void NRICfft(float *data,unsigned int nn,int isign) 
  100.     {
  101.     int istep,m,mmax,i,j;
  102.     unsigned int n;
  103.     double wtemp,wr,wpr,wpi,wi,theta;
  104.     float tempr,tempi;
  105.  
  106.         data -= 1; ;   /* Because the stupid routine is written for data[1...2*nn] */
  107.     n = nn << 1;
  108.     j=1;
  109.     for(i=1;i<n;i+=2) {
  110.         if( j>i) {
  111.             SWAP(data[j],data[i]);
  112.             SWAP(data[j+1],data[i+1]);
  113.             }
  114.         m = n >> 1;
  115.         while (m >= 2 && j > m) {
  116.             j -=m;
  117.             m >>= 1;
  118.             }
  119.         j +=m;
  120.         }
  121.     mmax =2;
  122.     while (n > mmax) {
  123.         istep = 2*mmax;
  124.         theta = PI2 / (isign*mmax);
  125.         wtemp = sin(0.5*theta);
  126.         wpr = -2.0*wtemp*wtemp;
  127.         wpi = sin(theta);
  128.         wr = 1.0;
  129.         wi = 0.0;
  130.         for ( m=1; m<mmax;m+=2) {
  131.             for(i=m;i<n;i+=istep) {
  132.                 j = i+mmax;
  133.                 tempr = wr*data[j]-wi*data[j+1];
  134.                 tempi = wr*data[j+1]+wi*data[j];
  135.                 data[j] = data[i]-tempr;
  136.                 data[j+1] = data[i+1]-tempi;
  137.                 data[i] += tempr;
  138.                 data[i+1] += tempi;
  139.                 }
  140.             wr = (wtemp=wr)*wpr-wi*wpi+wr;
  141.             wi = wi*wpr+wtemp*wpi+wi;
  142.             }
  143.         mmax = istep;
  144.         }
  145.  
  146.     }
  147.