server/scsynth/SC_Complex.cpp

   1 /*
   2         SuperCollider real time audio synthesis system
   3     Copyright (c) 2002 James McCartney. All rights reserved.
   4         http://www.audiosynth.com
   5
   6     This program is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     This program is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with this program; if not, write to the Free Software
  18     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301  USA
  19 */
  20
  21
  22 #include "SC_Complex.h"
  23 #include "SC_Constants.h"
  24 #include "SC_Samp.h"
  25 #include <math.h>
  26
  27 ////////////////////////////////////////////////////////////////////////////////
  28
  29 const int32 kPolarLUTSize = 2049;
  30 const int32 kPolarLUTSize2 = kPolarLUTSize >> 1;
  31 const double kSinePhaseScale = (double)kSineSize / twopi;
  32 float gMagLUT[kPolarLUTSize];
  33 float gPhaseLUT[kPolarLUTSize];
  34
  35 void BuildPolarLUT();
  36 void BuildPolarLUT()
  37 {
  38         double rPolarLUTSize2 = 1. / kPolarLUTSize2;
  39         for (int i=0; i < kPolarLUTSize; ++i) {
  40                 double slope = (i - kPolarLUTSize2) * rPolarLUTSize2;
  41                 double angle = atan(slope);
  42                 gPhaseLUT[i] = (float)angle;
  43                 gMagLUT[i] = 1.f / (float)cos(angle);
  44         }
  45 }
  46
  47 Polar Complex::ToPolar()
  48 {
  49         return Polar(hypot(imag, real), atan2(imag, real));
  50 }
  51
  52 Complex Polar::ToComplex()
  53 {
  54         return Complex(mag * cos(phase), mag * sin(phase));
  55 }
  56
  57 /*
  58 float fhypotx(float real, float imag);
  59 float fhypotx(float real, float imag)
  60 {
  61
  62         int32 index;
  63         float absreal = fabs(real);
  64         float absimag = fabs(imag);
  65         float slope;
  66         if (absreal > absimag) {
  67                 slope = imag/real;
  68                 index = kPolarLUTSize2 + kPolarLUTSize2 * slope;
  69                 return gMagLUT[index] * absreal;
  70         } else {
  71                 slope = real/imag;
  72                 index = kPolarLUTSize2 + kPolarLUTSize2 * slope;
  73                 return gMagLUT[index] * absimag;
  74         }
  75 }
  76 */
  77
  78 /**
  79 * Converts cartesian to polar representation, using lookup tables.
  80 * Note: in this implementation the phase values returned lie in the range [-pi/4, 7pi/4]
  81 * rather than the more conventional [0, 2pi] or [-pi, pi].
  82 */
  83 Polar Complex::ToPolarApx()
  84 {
  85         int32 index;
  86         float absreal = fabs(real);
  87         float absimag = fabs(imag);
  88         float mag, phase, slope;
  89         if (absreal > absimag) {
  90                 slope = imag/real;
  91                 index = (int32)(kPolarLUTSize2 + kPolarLUTSize2 * slope);
  92                 mag = gMagLUT[index] * absreal;
  93                 phase = gPhaseLUT[index];
  94                 if (real > 0) {
  95                         return Polar(mag, phase);
  96                 } else {
  97                         return Polar(mag, (float)(pi + phase));
  98                 }
  99         } else {
 100                 slope = real/imag;
 101                 index = (int32)(kPolarLUTSize2 + kPolarLUTSize2 * slope);
 102                 mag = gMagLUT[index] * absimag;
 103                 phase = gPhaseLUT[index];
 104                 if (imag > 0) {
 105                         return Polar(mag, (float)(pi2 - phase));
 106                 } else {
 107                         return Polar(mag, (float)(pi32 - phase));
 108                 }
 109         }
 110 }
 111
 112 void Complex::ToPolarInPlace()
 113 {
 114         Polar polar = ToPolar();
 115         real = polar.mag;
 116         imag = polar.phase;
 117 }
 118
 119 void Polar::ToComplexInPlace()
 120 {
 121         Complex complx = ToComplex();
 122         mag = complx.real;
 123         phase = complx.imag;
 124 }
 125
 126 Complex Polar::ToComplexApx()
 127 {
 128         uint32 sinindex = (int32)(kSinePhaseScale * phase) & kSineMask;
 129         uint32 cosindex = (sinindex + (kSineSize>>2)) & kSineMask;
 130         return Complex(mag * gSine[cosindex], mag * gSine[sinindex]);
 131 }
 132
 133 void Complex::ToPolarApxInPlace()
 134 {
 135         Polar polar = ToPolarApx();
 136         real = polar.mag;
 137         imag = polar.phase;
 138 }
 139
 140 void Polar::ToComplexApxInPlace()
 141 {
 142         Complex complx = ToComplexApx();
 143         mag = complx.real;
 144         phase = complx.imag;
 145 }
 146
 147 ////////////////////////////////////////////////////////////////////////////////
 148