include/plugin_interface/SC_Complex.h

   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 #ifndef _SC_Complex_
  23 #define _SC_Complex_
  24
  25 #include <cmath>
  26
  27 #include "SC_Types.h"
  28 #include "SC_Constants.h"
  29 #include "float.h"
  30
  31 #ifdef _MSC_VER
  32 // hypotf is c99, but not c++
  33 #define hypotf _hypotf
  34 #endif
  35
  36 ////////////////////////////////////////////////////////////////////////////////
  37
  38 namespace detail {
  39
  40 const int kSineSize = 8192;
  41 const int kSineMask = kSineSize - 1;
  42 const double kSinePhaseScale = kSineSize / twopi;
  43 const int32 kPolarLUTSize = 2049;
  44 const int32 kPolarLUTSize2 = kPolarLUTSize >> 1;
  45
  46
  47 /* each object file that is including this header will have separate lookup tables */
  48 namespace {
  49
  50 float gMagLUT[kPolarLUTSize];
  51 float gPhaseLUT[kPolarLUTSize];
  52 float gSine[kSineSize+1];
  53
  54 static bool initTables(void)
  55 {
  56         double sineIndexToPhase = twopi / kSineSize;
  57         double pmf = (1L << 29) / twopi;
  58         for (int i=0; i <= kSineSize; ++i) {
  59                 double phase = i * sineIndexToPhase;
  60                 float32 d = sin(phase);
  61                 gSine[i] = d;
  62         }
  63
  64         double rPolarLUTSize2 = 1. / kPolarLUTSize2;
  65         for (int i=0; i < kPolarLUTSize; ++i) {
  66                 double slope = (i - kPolarLUTSize2) * rPolarLUTSize2;
  67                 double angle = atan(slope);
  68                 gPhaseLUT[i] = (float)angle;
  69                 gMagLUT[i] = (float)(1.f / cos(angle));
  70         }
  71
  72         return true;
  73 }
  74
  75 bool dummy = initTables();
  76
  77 }
  78
  79 struct Polar;
  80
  81
  82 struct Complex
  83 {
  84         Complex() {}
  85         Complex(float r, float i) : real(r), imag(i) {}
  86         void Set(float r, float i) { real = r; imag = i; }
  87
  88         Complex& operator=(Complex b) { real = b.real; imag = b.imag; return *this; }
  89         Complex& operator=(float b) { real = b; imag = 0.; return *this; }
  90
  91         Polar ToPolar();
  92
  93         /**
  94         * Converts cartesian to polar representation, using lookup tables.
  95         * Note: in this implementation the phase values returned lie in the range [-pi/4, 7pi/4]
  96         * rather than the more conventional [0, 2pi] or [-pi, pi].
  97         */
  98         Polar ToPolarApx();
  99
 100         void ToPolarInPlace();
 101
 102         void ToPolarApxInPlace();
 103
 104         float real, imag;
 105 };
 106
 107 struct Polar
 108 {
 109         Polar() {}
 110         Polar(float m, float p) : mag(m), phase(p) {}
 111         void Set(float m, float p) { mag = m; phase = p; }
 112
 113         Complex ToComplex()
 114     {
 115         return Complex(mag * std::cos(phase), mag * std::sin(phase));
 116     }
 117
 118         Complex ToComplexApx()
 119         {
 120                 uint32 sinindex = (int32)(kSinePhaseScale * phase) & kSineMask;
 121                 uint32 cosindex = (sinindex + (kSineSize>>2)) & kSineMask;
 122                 return Complex(mag * gSine[cosindex], mag * gSine[sinindex]);
 123         }
 124
 125         void ToComplexInPlace()
 126         {
 127                 Complex complx = ToComplex();
 128                 mag = complx.real;
 129                 phase = complx.imag;
 130         }
 131
 132         void ToComplexApxInPlace()
 133         {
 134                 Complex complx = ToComplexApx();
 135                 mag = complx.real;
 136                 phase = complx.imag;
 137         }
 138
 139         float mag, phase;
 140 };
 141
 142 inline Polar Complex::ToPolar()
 143 {
 144         return Polar(hypotf(imag, real), std::atan2(imag, real));
 145 }
 146
 147 inline Polar Complex::ToPolarApx()
 148 {
 149         int32 index;
 150         float absreal = fabs(real);
 151         float absimag = fabs(imag);
 152         float mag, phase, slope;
 153         if (absreal > absimag) {
 154                 slope = imag/real;
 155                 index = (int32)(kPolarLUTSize2 + kPolarLUTSize2 * slope);
 156                 mag = gMagLUT[index] * absreal;
 157                 phase = gPhaseLUT[index];
 158                 if (real > 0) {
 159                         return Polar(mag, phase);
 160                 } else {
 161                         return Polar(mag, (float)(pi + phase));
 162                 }
 163         } else if (absimag > 0) {
 164                 slope = real/imag;
 165                 index = (int32)(kPolarLUTSize2 + kPolarLUTSize2 * slope);
 166                 mag = gMagLUT[index] * absimag;
 167                 phase = gPhaseLUT[index];
 168                 if (imag > 0) {
 169                         return Polar(mag, (float)(pi2 - phase));
 170                 } else {
 171                         return Polar(mag, (float)(pi32 - phase));
 172                 }
 173         } else return Polar(0, 0);
 174 }
 175
 176 inline void Complex::ToPolarInPlace()
 177 {
 178         Polar polar = ToPolar();
 179         real = polar.mag;
 180         imag = polar.phase;
 181 }
 182
 183 inline void Complex::ToPolarApxInPlace()
 184 {
 185         Polar polar = ToPolarApx();
 186         real = polar.mag;
 187         imag = polar.phase;
 188 }
 189
 190 }
 191
 192 using detail::Polar;
 193 using detail::Complex;
 194
 195 struct ComplexFT
 196 {
 197         float dc, nyq;
 198         Complex complex[1];
 199 };
 200
 201 struct PolarFT
 202 {
 203         float dc, nyq;
 204         Polar polar[1];
 205 };
 206
 207 void ToComplex(Polar in, Complex& out);
 208
 209 inline Complex operator+(Complex a, Complex b) { return Complex(a.real + b.real, a.imag + b.imag); }
 210 inline Complex operator+(Complex a, float b) { return Complex(a.real + b, a.imag); }
 211 inline Complex operator+(float a, Complex b) { return Complex(a + b.real, b.imag); }
 212
 213 inline Complex& operator+=(Complex& a, const Complex& b) { a.real += b.real, a.imag += b.imag; return a; }
 214 inline Complex& operator+=(Complex& a, float b) { a.real += b; return a; }
 215
 216 inline Complex operator-(Complex a, Complex b) { return Complex(a.real - b.real, a.imag - b.imag); }
 217 inline Complex operator-(Complex a, float b) { return Complex(a.real - b, a.imag); }
 218 inline Complex operator-(float a, Complex b) { return Complex(a - b.real, b.imag); }
 219
 220 inline Complex operator-=(Complex a, Complex b) { a.real -= b.real, a.imag -= b.imag; return a; }
 221 inline Complex operator-=(Complex a, float b) { a.real -= b; return a; }
 222
 223 inline Complex operator*(Complex a, Complex b)
 224 {
 225         return Complex(a.real * b.real - a.imag * b.imag, a.real * b.imag + a.imag * b.real);
 226 }
 227
 228 inline Complex operator*(Complex a, float b)
 229 {
 230         return Complex(a.real * b, a.imag * b);
 231 }
 232
 233 inline Complex operator*(float a, Complex b)
 234 {
 235         return Complex(b.real * a, b.imag * a);
 236 }
 237
 238 inline Complex operator*=(Complex a, Complex b)
 239 {
 240         a.Set(
 241                 a.real * b.real - a.imag * b.imag,
 242                 a.real * b.imag + a.imag * b.real
 243         );
 244         return a;
 245 }
 246
 247 inline Complex operator*=(Complex a, float b)
 248 {
 249         a.real *= b;
 250         a.imag *= b;
 251         return a;
 252 }
 253
 254
 255 inline Polar operator*(Polar a, float b)
 256 {
 257         return Polar(a.mag * b, a.phase);
 258 }
 259
 260 inline Polar operator*(float a, Polar b)
 261 {
 262         return Polar(a * b.mag, b.phase);
 263 }
 264
 265 inline Polar operator*=(Polar a, float b)
 266 {
 267         a.mag *= b;
 268         return a;
 269 }
 270
 271 #endif