v2/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 //-----------------------------------------------------------------------------
  86 // name: blackman()
  87 // desc: make window
  88 //-----------------------------------------------------------------------------
  89 void blackman( FLOAT * window, unsigned long length )
  90 {
  91     unsigned long i;
  92     double pi, phase = 0, delta;
  93
  94     pi = 4.*atan(1.0);
  95     delta = 2 * pi / (double) length;
  96
  97     for( i = 0; i < length; i++ )
  98     {
  99         window[i] = (FLOAT)(0.42 - .5*cos(phase) + .08*cos(2*phase));
 100         phase += delta;
 101     }
 102 }
 103
 104
 105
 106
 107 //-----------------------------------------------------------------------------
 108 // name: bartlett()
 109 // desc: make window
 110 //-----------------------------------------------------------------------------
 111 void bartlett( FLOAT * window, unsigned long length )
 112 {
 113     unsigned long i;
 114     FLOAT half = (FLOAT)length / 2;
 115
 116     for( i = 0; i < length; i++ )
 117     {
 118         if( i < half ) window[i] = i / half;
 119         else window[i] = (length - i) / half;
 120     }
 121 }
 122
 123
 124
 125
 126 //-----------------------------------------------------------------------------
 127 // name: apply_window()
 128 // desc: apply a window to data
 129 //-----------------------------------------------------------------------------
 130 void apply_window( FLOAT * data, FLOAT * window, unsigned long length )
 131 {
 132     unsigned long i;
 133
 134     for( i = 0; i < length; i++ )
 135         data[i] *= window[i];
 136 }
 137
 138 static FLOAT PI ;
 139 static FLOAT TWOPI ;
 140 void bit_reverse( FLOAT * x, long N );
 141
 142 //-----------------------------------------------------------------------------
 143 // name: rfft()
 144 // desc: real value fft
 145 //
 146 //   these routines from the CARL software, spect.c
 147 //   check out the CARL CMusic distribution for more source code
 148 //
 149 //   if forward is true, rfft replaces 2*N real data points in x with N complex
 150 //   values representing the positive frequency half of their Fourier spectrum,
 151 //   with x[1] replaced with the real part of the Nyquist frequency value.
 152 //
 153 //   if forward is false, rfft expects x to contain a positive frequency
 154 //   spectrum arranged as before, and replaces it with 2*N real values.
 155 //
 156 //   N MUST be a power of 2.
 157 //
 158 //-----------------------------------------------------------------------------
 159 void rfft( FLOAT * x, long N, unsigned int forward )
 160 {
 161     static int first = 1 ;
 162     FLOAT c1, c2, h1r, h1i, h2r, h2i, wr, wi, wpr, wpi, temp, theta ;
 163     FLOAT xr, xi ;
 164     long i, i1, i2, i3, i4, N2p1 ;
 165
 166     if( first )
 167     {
 168         PI = (FLOAT) (4.*atan( 1. )) ;
 169         TWOPI = (FLOAT) (8.*atan( 1. )) ;
 170         first = 0 ;
 171     }
 172
 173     theta = PI/N ;
 174     wr = 1. ;
 175     wi = 0. ;
 176     c1 = 0.5 ;
 177
 178     if( forward )
 179     {
 180         c2 = -0.5 ;
 181         cfft( x, N, forward ) ;
 182         xr = x[0] ;
 183         xi = x[1] ;
 184     }
 185     else
 186     {
 187         c2 = 0.5 ;
 188         theta = -theta ;
 189         xr = x[1] ;
 190         xi = 0. ;
 191         x[1] = 0. ;
 192     }
 193
 194     wpr = (FLOAT) (-2.*pow( sin( 0.5*theta ), 2. )) ;
 195     wpi = (FLOAT) sin( theta ) ;
 196     N2p1 = (N<<1) + 1 ;
 197
 198     for( i = 0 ; i <= N>>1 ; i++ )
 199     {
 200         i1 = i<<1 ;
 201         i2 = i1 + 1 ;
 202         i3 = N2p1 - i2 ;
 203         i4 = i3 + 1 ;
 204         if( i == 0 )
 205         {
 206             h1r =  c1*(x[i1] + xr ) ;
 207             h1i =  c1*(x[i2] - xi ) ;
 208             h2r = -c2*(x[i2] + xi ) ;
 209             h2i =  c2*(x[i1] - xr ) ;
 210             x[i1] =  h1r + wr*h2r - wi*h2i ;
 211             x[i2] =  h1i + wr*h2i + wi*h2r ;
 212             xr =  h1r - wr*h2r + wi*h2i ;
 213             xi = -h1i + wr*h2i + wi*h2r ;
 214         }
 215         else
 216         {
 217             h1r =  c1*(x[i1] + x[i3] ) ;
 218             h1i =  c1*(x[i2] - x[i4] ) ;
 219             h2r = -c2*(x[i2] + x[i4] ) ;
 220             h2i =  c2*(x[i1] - x[i3] ) ;
 221             x[i1] =  h1r + wr*h2r - wi*h2i ;
 222             x[i2] =  h1i + wr*h2i + wi*h2r ;
 223             x[i3] =  h1r - wr*h2r + wi*h2i ;
 224             x[i4] = -h1i + wr*h2i + wi*h2r ;
 225         }
 226
 227         wr = (temp = wr)*wpr - wi*wpi + wr ;
 228         wi = wi*wpr + temp*wpi + wi ;
 229     }
 230
 231     if( forward )
 232         x[1] = xr ;
 233     else
 234         cfft( x, N, forward ) ;
 235 }
 236
 237
 238
 239
 240 //-----------------------------------------------------------------------------
 241 // name: cfft()
 242 // desc: complex value fft
 243 //
 244 //   these routines from CARL software, spect.c
 245 //   check out the CARL CMusic distribution for more software
 246 //
 247 //   cfft replaces FLOAT array x containing NC complex values (2*NC FLOAT
 248 //   values alternating real, imagininary, etc.) by its Fourier transform
 249 //   if forward is true, or by its inverse Fourier transform ifforward is
 250 //   false, using a recursive Fast Fourier transform method due to
 251 //   Danielson and Lanczos.
 252 //
 253 //   NC MUST be a power of 2.
 254 //
 255 //-----------------------------------------------------------------------------
 256 void cfft( FLOAT * x, long NC, unsigned int forward )
 257 {
 258     FLOAT wr, wi, wpr, wpi, theta, scale ;
 259     long mmax, ND, m, i, j, delta ;
 260     ND = NC<<1 ;
 261     bit_reverse( x, ND ) ;
 262
 263     for( mmax = 2 ; mmax < ND ; mmax = delta )
 264     {
 265         delta = mmax<<1 ;
 266         theta = TWOPI/( forward? mmax : -mmax ) ;
 267         wpr = (FLOAT) (-2.*pow( sin( 0.5*theta ), 2. )) ;
 268         wpi = (FLOAT) sin( theta ) ;
 269         wr = 1. ;
 270         wi = 0. ;
 271
 272         for( m = 0 ; m < mmax ; m += 2 )
 273         {
 274             register FLOAT rtemp, itemp ;
 275             for( i = m ; i < ND ; i += delta )
 276             {
 277                 j = i + mmax ;
 278                 rtemp = wr*x[j] - wi*x[j+1] ;
 279                 itemp = wr*x[j+1] + wi*x[j] ;
 280                 x[j] = x[i] - rtemp ;
 281                 x[j+1] = x[i+1] - itemp ;
 282                 x[i] += rtemp ;
 283                 x[i+1] += itemp ;
 284             }
 285
 286             wr = (rtemp = wr)*wpr - wi*wpi + wr ;
 287             wi = wi*wpr + rtemp*wpi + wi ;
 288         }
 289     }
 290
 291     // scale output
 292     scale = (FLOAT)(forward ? 1./ND : 2.) ;
 293     {
 294         register FLOAT *xi=x, *xe=x+ND ;
 295         while( xi < xe )
 296             *xi++ *= scale ;
 297     }
 298 }
 299
 300
 301
 302
 303 //-----------------------------------------------------------------------------
 304 // name: bit_reverse()
 305 // desc: bitreverse places FLOAT array x containing N/2 complex values
 306 //       into bit-reversed order
 307 //-----------------------------------------------------------------------------
 308 void bit_reverse( FLOAT * x, long N )
 309 {
 310     FLOAT rtemp, itemp ;
 311     long i, j, m ;
 312     for( i = j = 0 ; i < N ; i += 2, j += m )
 313     {
 314         if( j > i )
 315         {
 316             rtemp = x[j] ; itemp = x[j+1] ; /* complex exchange */
 317             x[j] = x[i] ; x[j+1] = x[i+1] ;
 318             x[i] = rtemp ; x[i+1] = itemp ;
 319         }
 320
 321         for( m = N>>1 ; m >= 2 && j >= m ; m >>= 1 )
 322             j -= m ;
 323     }
 324 }
 325
 326
 327
 328
 329 //-----------------------------------------------------------------------------
 330 // name: the_dct()
 331 // desc: type ii dct on N reals
 332 //-----------------------------------------------------------------------------
 333 void the_dct( FLOAT * x, unsigned long N, FLOAT * out, unsigned long Nout )
 334 {
 335     unsigned long k, n;
 336
 337     // sanity check
 338     assert( Nout <= N );
 339     // zero out
 340     memset( out, 0, sizeof(FLOAT)*Nout );
 341
 342     // go for it
 343     for( k = 0; k < Nout; k++ )
 344     {
 345         for( n = 0; n < N; n++ )
 346         {
 347             out[k] += x[n] * cos( ONE_PI / N * k * (n + .5) );
 348         }
 349     }
 350 }
 351
 352
 353
 354
 355
 356 //-----------------------------------------------------------------------------
 357 // name: the_dct_matrix()
 358 // desc: type ii dct matrx on N reals
 359 //-----------------------------------------------------------------------------
 360 void the_dct_matrix( FLOAT ** out, unsigned long N )
 361 {
 362     unsigned long k, n;
 363
 364     // go for it
 365     for( k = 0; k < N; k++ )
 366     {
 367         for( n = 0; n < N; n++ )
 368         {
 369             out[k][n] = cos( ONE_PI / N * k * (n + .5) );
 370         }
 371     }
 372 }
 373
 374
 375 //-----------------------------------------------------------------------------
 376 // name: the_dct_now()
 377 // desc: apply dct matrx on N reals
 378 //-----------------------------------------------------------------------------
 379 void the_dct_now( FLOAT * x, FLOAT ** matrix, unsigned long N,
 380                   FLOAT * out, unsigned long Nout )
 381 {
 382     unsigned long k, n;
 383
 384     // sanity check
 385     assert( Nout <= N );
 386     // zero out
 387     memset( out, 0, sizeof(FLOAT)*Nout );
 388
 389     // go for it
 390     for( k = 0; k < Nout; k++ )
 391     {
 392         for( n = 0; n < N; n++ )
 393         {
 394             out[k] = x[n] * matrix[k][n];
 395         }
 396     }
 397 }
 398
 399
 400
 401
 402 //-----------------------------------------------------------------------------
 403 // name: the_inverse_dct()
 404 // desc: type iii dct, or inverse dct
 405 //-----------------------------------------------------------------------------
 406 void the_inverse_dct( FLOAT * x, unsigned long N, FLOAT * out, unsigned long Nout )
 407 {
 408 }
 409
 410
 411 //-----------------------------------------------------------------------------
 412 // name: the_inverse_dct_matrix()
 413 // desc: generates the NxN DCT matrix
 414 //-----------------------------------------------------------------------------
 415 void the_inverse_dct_matrix( FLOAT ** out, unsigned long N )
 416 {
 417     unsigned long k, n;
 418
 419     // go for it
 420     for( k = 0; k < N; k++ )
 421     {
 422         for( n = 0; n < N; n++ )
 423         {
 424             out[k][n] = cos( ONE_PI / N * n * (k + .5) );
 425         }
 426     }
 427 }
 428
 429
 430
 431
 432 //-----------------------------------------------------------------------------
 433 // name: the_inverse_dct_now()
 434 // desc: apply inverse dct matrx on N reals
 435 //-----------------------------------------------------------------------------
 436 void the_inverse_dct_now( FLOAT * x, FLOAT ** matrix, unsigned long N,
 437                   FLOAT * out, unsigned long Nout )
 438 {
 439     unsigned long k, n;
 440
 441     // sanity check
 442     assert( Nout <= N );
 443     // zero out
 444     memset( out, 0, sizeof(FLOAT)*Nout );
 445
 446     // go for it
 447     for( k = 0; k < Nout; k++ )
 448     {
 449         for( n = 0; n < N; n++ )
 450         {
 451             out[k] = x[0] / 2 + x[n] * matrix[k][n];
 452         }
 453     }
 454 }