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