exile/v1/src/util_xforms.c

   1 /*----------------------------------------------------------------------------
   2     ChucK Concurrent, On-the-fly Audio Programming Language
   3       Compiler and Virtual Machine
   4
   5     Copyright (c) 2004 Ge Wang and Perry R. Cook.  All rights reserved.
   6       http://chuck.cs.princeton.edu/
   7       http://soundlab.cs.princeton.edu/
   8
   9     This program is free software; you can redistribute it and/or modify
  10     it under the terms of the GNU General Public License as published by
  11     the Free Software Foundation; either version 2 of the License, or
  12     (at your option) any later version.
  13
  14     This program is distributed in the hope that it will be useful,
  15     but WITHOUT ANY WARRANTY; without even the implied warranty of
  16     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17     GNU General Public License for more details.
  18
  19     You should have received a copy of the GNU General Public License
  20     along with this program; if not, write to the Free Software
  21     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
  22     U.S.A.
  23 -----------------------------------------------------------------------------*/
  24
  25 //-----------------------------------------------------------------------------
  26 // name: util_xform.c
  27 // desc: transforms
  28 //
  29 // authors: Ge Wang (gewang@cs.princeton.edu)
  30 //          Perry R. Cook (prc@cs.princeton.edu)
  31 //          FFT: CARL music distribution
  32 // date: 11.27.2003
  33 //-----------------------------------------------------------------------------
  34 #include "util_xforms.h"
  35 #include <stdlib.h>
  36 #include <math.h>
  37
  38
  39
  40
  41 //-----------------------------------------------------------------------------
  42 // name: hanning()
  43 // desc: make window
  44 //-----------------------------------------------------------------------------
  45 void hanning( float * window, unsigned long length )
  46 {
  47     unsigned long i;
  48     double pi, phase = 0, delta;
  49
  50     pi = 4.*atan(1.0);
  51     delta = 2 * pi / (double) length;
  52
  53     for( i = 0; i < length; i++ )
  54     {
  55         window[i] = (float)(0.5 * (1.0 - cos(phase)));
  56         phase += delta;
  57     }
  58 }
  59
  60
  61
  62
  63 //-----------------------------------------------------------------------------
  64 // name: hamming()
  65 // desc: make window
  66 //-----------------------------------------------------------------------------
  67 void hamming( float * window, unsigned long length )
  68 {
  69     unsigned long i;
  70     double pi, phase = 0, delta;
  71
  72     pi = 4.*atan(1.0);
  73     delta = 2 * pi / (double) length;
  74
  75     for( i = 0; i < length; i++ )
  76     {
  77         window[i] = (float)(0.54 - .46*cos(phase));
  78         phase += delta;
  79     }
  80 }
  81
  82
  83
  84 //-----------------------------------------------------------------------------
  85 // name: blackman()
  86 // desc: make window
  87 //-----------------------------------------------------------------------------
  88 void blackman( float * window, unsigned long length )
  89 {
  90     unsigned long i;
  91     double pi, phase = 0, delta;
  92
  93     pi = 4.*atan(1.0);
  94     delta = 2 * pi / (double) length;
  95
  96     for( i = 0; i < length; i++ )
  97     {
  98         window[i] = (float)(0.42 - .5*cos(phase) + .08*cos(2*phase));
  99         phase += delta;
 100     }
 101 }
 102
 103
 104
 105
 106 //-----------------------------------------------------------------------------
 107 // name: apply_window()
 108 // desc: apply a window to data
 109 //-----------------------------------------------------------------------------
 110 void apply_window( float * data, float * window, unsigned long length )
 111 {
 112     unsigned long i;
 113
 114     for( i = 0; i < length; i++ )
 115         data[i] *= window[i];
 116 }
 117
 118 static float PI ;
 119 static float TWOPI ;
 120 void bit_reverse( float * x, long N );
 121
 122 //-----------------------------------------------------------------------------
 123 // name: rfft()
 124 // desc: real value fft
 125 //
 126 //   these routines from the CARL software, spect.c
 127 //   check out the CARL CMusic distribution for more source code
 128 //
 129 //   if forward is true, rfft replaces 2*N real data points in x with N complex
 130 //   values representing the positive frequency half of their Fourier spectrum,
 131 //   with x[1] replaced with the real part of the Nyquist frequency value.
 132 //
 133 //   if forward is false, rfft expects x to contain a positive frequency
 134 //   spectrum arranged as before, and replaces it with 2*N real values.
 135 //
 136 //   N MUST be a power of 2.
 137 //
 138 //-----------------------------------------------------------------------------
 139 void rfft( float * x, long N, unsigned int forward )
 140 {
 141     static int first = 1 ;
 142     float c1, c2, h1r, h1i, h2r, h2i, wr, wi, wpr, wpi, temp, theta ;
 143     float xr, xi ;
 144     long i, i1, i2, i3, i4, N2p1 ;
 145
 146     if( first )
 147     {
 148         PI = (float) (4.*atan( 1. )) ;
 149         TWOPI = (float) (8.*atan( 1. )) ;
 150         first = 0 ;
 151     }
 152
 153     theta = PI/N ;
 154     wr = 1. ;
 155     wi = 0. ;
 156     c1 = 0.5 ;
 157
 158     if( forward )
 159     {
 160         c2 = -0.5 ;
 161         cfft( x, N, forward ) ;
 162         xr = x[0] ;
 163         xi = x[1] ;
 164     }
 165     else
 166     {
 167         c2 = 0.5 ;
 168         theta = -theta ;
 169         xr = x[1] ;
 170         xi = 0. ;
 171         x[1] = 0. ;
 172     }
 173
 174     wpr = (float) (-2.*pow( sin( 0.5*theta ), 2. )) ;
 175     wpi = (float) sin( theta ) ;
 176     N2p1 = (N<<1) + 1 ;
 177
 178     for( i = 0 ; i <= N>>1 ; i++ )
 179     {
 180         i1 = i<<1 ;
 181         i2 = i1 + 1 ;
 182         i3 = N2p1 - i2 ;
 183         i4 = i3 + 1 ;
 184         if( i == 0 )
 185         {
 186             h1r =  c1*(x[i1] + xr ) ;
 187             h1i =  c1*(x[i2] - xi ) ;
 188             h2r = -c2*(x[i2] + xi ) ;
 189             h2i =  c2*(x[i1] - xr ) ;
 190             x[i1] =  h1r + wr*h2r - wi*h2i ;
 191             x[i2] =  h1i + wr*h2i + wi*h2r ;
 192             xr =  h1r - wr*h2r + wi*h2i ;
 193             xi = -h1i + wr*h2i + wi*h2r ;
 194         }
 195         else
 196         {
 197             h1r =  c1*(x[i1] + x[i3] ) ;
 198             h1i =  c1*(x[i2] - x[i4] ) ;
 199             h2r = -c2*(x[i2] + x[i4] ) ;
 200             h2i =  c2*(x[i1] - x[i3] ) ;
 201             x[i1] =  h1r + wr*h2r - wi*h2i ;
 202             x[i2] =  h1i + wr*h2i + wi*h2r ;
 203             x[i3] =  h1r - wr*h2r + wi*h2i ;
 204             x[i4] = -h1i + wr*h2i + wi*h2r ;
 205         }
 206
 207         wr = (temp = wr)*wpr - wi*wpi + wr ;
 208         wi = wi*wpr + temp*wpi + wi ;
 209     }
 210
 211     if( forward )
 212         x[1] = xr ;
 213     else
 214         cfft( x, N, forward ) ;
 215 }
 216
 217
 218
 219
 220 //-----------------------------------------------------------------------------
 221 // name: cfft()
 222 // desc: complex value fft
 223 //
 224 //   these routines from CARL software, spect.c
 225 //   check out the CARL CMusic distribution for more software
 226 //
 227 //   cfft replaces float array x containing NC complex values (2*NC float
 228 //   values alternating real, imagininary, etc.) by its Fourier transform
 229 //   if forward is true, or by its inverse Fourier transform ifforward is
 230 //   false, using a recursive Fast Fourier transform method due to
 231 //   Danielson and Lanczos.
 232 //
 233 //   NC MUST be a power of 2.
 234 //
 235 //-----------------------------------------------------------------------------
 236 void cfft( float * x, long NC, unsigned int forward )
 237 {
 238     float wr, wi, wpr, wpi, theta, scale ;
 239     long mmax, ND, m, i, j, delta ;
 240     ND = NC<<1 ;
 241     bit_reverse( x, ND ) ;
 242
 243     for( mmax = 2 ; mmax < ND ; mmax = delta )
 244     {
 245         delta = mmax<<1 ;
 246         theta = TWOPI/( forward? mmax : -mmax ) ;
 247         wpr = (float) (-2.*pow( sin( 0.5*theta ), 2. )) ;
 248         wpi = (float) sin( theta ) ;
 249         wr = 1. ;
 250         wi = 0. ;
 251
 252         for( m = 0 ; m < mmax ; m += 2 )
 253         {
 254             register float rtemp, itemp ;
 255             for( i = m ; i < ND ; i += delta )
 256             {
 257                 j = i + mmax ;
 258                 rtemp = wr*x[j] - wi*x[j+1] ;
 259                 itemp = wr*x[j+1] + wi*x[j] ;
 260                 x[j] = x[i] - rtemp ;
 261                 x[j+1] = x[i+1] - itemp ;
 262                 x[i] += rtemp ;
 263                 x[i+1] += itemp ;
 264             }
 265
 266             wr = (rtemp = wr)*wpr - wi*wpi + wr ;
 267             wi = wi*wpr + rtemp*wpi + wi ;
 268         }
 269     }
 270
 271     // scale output
 272     scale = (float)(forward ? 1./ND : 2.) ;
 273     {
 274         register float *xi=x, *xe=x+ND ;
 275         while( xi < xe )
 276             *xi++ *= scale ;
 277     }
 278 }
 279
 280
 281
 282
 283 //-----------------------------------------------------------------------------
 284 // name: bit_reverse()
 285 // desc: bitreverse places float array x containing N/2 complex values
 286 //       into bit-reversed order
 287 //-----------------------------------------------------------------------------
 288 void bit_reverse( float * x, long N )
 289 {
 290     float rtemp, itemp ;
 291     long i, j, m ;
 292     for( i = j = 0 ; i < N ; i += 2, j += m )
 293     {
 294         if( j > i )
 295         {
 296             rtemp = x[j] ; itemp = x[j+1] ; /* complex exchange */
 297             x[j] = x[i] ; x[j+1] = x[i+1] ;
 298             x[i] = rtemp ; x[i+1] = itemp ;
 299         }
 300
 301         for( m = N>>1 ; m >= 2 && j >= m ; m >>= 1 )
 302             j -= m ;
 303     }
 304 }