// Sweetdreams: a toy softsynth in C++.

/* To run with sox:
make -k sweetdreams-c++ LDFLAGS=-lm && ./sweetdreams-c++ | play -t raw -r 48000 -s -b 16 -L -c 2 -
*/

// This program is in the public domain.  To the extent possible under
// law, Kragen Javier Sitaker has waived all copyright and related or
// neighboring rights to Sweetdreams softsynth. This work is published
// from Argentina; see
// <http://creativecommons.org/publicdomain/zero/1.0/> for details.

// That awesome synth sound in Eurythmics’ “Sweet Dreams” — is it just
// a flanged triangle wave?  Or does it maybe have some AM and/or FM
// in it?  I wanted to find out, so I wrote this stateless modular
// synth library in C++.  It turns out to be a bit more complicated
// than that; basically you take some harmonic-rich carrier wave,
// such as a sawtooth (not a triangle wave!),
// flange it so the harmonics change over time, and amplitude-modulate
// it at a lower frequency.  That isn’t the whole story, but that’s
// the basic timbre.

// The signal flow graph here to generate the Eurythmics synth
// consists of 29 nodes.  Executing it interpretively takes about 669
// instructions per output sample, about 23 instructions per node.

// Emits 48ksps (16-bit little-endian 2-channel signed, like DAT).

// (I started writing this in C, but man, you know what?  Fuck C.  I’m
// going to use C++.)

// Currently 43,854,340 instrutions and 32,093,296 data refs to
// produce only 65536 samples.  For shame!

// After batching up the output into 64-sample chunks: 35,662,215 and
// 28,883,826.

// After making Compose batch: 35,254,223 and 28,448,004 (1.2%
// faster).

// After making Amplify batch: 35,362,767 (!) and 28,085,508.  No
// improvement, at least.

// After making Constant batch: 34,826,191 and 27,684,100.  1.2%
// faster than after Compose.

// After making Identity batch: 34,355,151 and 27,282,692.  1.3%
// faster.

// After fixing the test in Compose to allow rounding errors:
// 34,371,368 and 27,269,996.

// After making Sum batch: 34,290,603 and 27,016,662.  However, it
// only batched 96 times.

// After loosening the rounding error test to 0.5%: 31,448,455 and
// 15,648,626.  9% better!

// After making Repeat batch: 25,137,245 and 10,311,366.  Awesome!
// 20% better!

// After making Triangle batch: 23,778,067 and 9,529,450.  6.4%
// better.

// Changing buffer_size from 64 to 32 or 128 didn’t help; it slightly
// worsened performance, by like 3% or so.

// After manually strength-reducing the base-class get loop:
// 23,776,979 and 9,873,352.  This is a tiny win, but I think it shows
// that the compiler isn’t doing the strength-reduction automatically,
// maybe because float addition isn’t a safe substitute for
// multiplication in general (though in this case it is).

// After eliminating the second buffer in Sum, 23,771,191.
// After eliminating the second buffer in Amplify, 23,761,975.
// These are not improvements in performance but they do simplify the
// code.

#include <cmath>
#include <cstdio>
#include <cstdint>
#include <map>
#include <memory>
#include <set>
#include <string>

using std::enable_shared_from_this;
using std::map;
using std::make_pair;
using std::make_shared;
using std::pair;
using std::set;
using std::shared_ptr;
using std::string;
using std::to_string;

class StdioGraph {
  FILE *out;
  int next_node_id;
  map<uintptr_t, int> node_ids;
  map<int, string> node_names;
  set<pair<int, int> > edges;
  map<pair<int, int>, string> edge_labels;

  int get_id(void *node) {
    auto node_int = uintptr_t(node);
    auto p = node_ids.find(node_int);
    if (p != node_ids.end()) return p->second;
    int id = next_node_id;
    next_node_id++;
    node_ids.insert(make_pair(node_int, id));
    return id;
  }

public:
  StdioGraph(FILE *out_) : out(out_), next_node_id(0) {}
  bool start() {return true;}
  bool finish() {
    if (0 > fprintf(out, "digraph sweetdreams {\n")) return false;

    for (auto it = node_names.begin(); it != node_names.end(); it++) {
      if (0 > fprintf(out, "  %d [label=\"%s\"];\n",
                      it->first, it->second.c_str())) return false;
    }

    for (auto it = edges.begin(); it != edges.end(); it++) {
      auto p = edge_labels.find(*it);
      if (p == edge_labels.end()) {
        if (0 > fprintf(out, "  %d -> %d;\n",
                        it->first, it->second)) return false;
      } else {
        if (0 > fprintf(out, "  %d -> %d [label=\"%s\"];\n",
                        it->first, it->second, p->second.c_str())) return false;
      }
    }

    if (0 > fprintf(out, "}\n")) return false;
    return true;
  }

  bool node_name(void *node, string name) {
    int id = get_id(node);
    node_names.insert(make_pair(id, name));
    return true;
  }

  bool node_edge(void *from, void *to, string label = "") {
    int id_from = get_id(from), id_to = get_id(to);
    auto pair = make_pair(id_from, id_to);
    edges.insert(pair);
    if (label != "") {
      edge_labels.insert(make_pair(pair, label));
    }
    return true;
  }
};

class Signal;

// Previously I just had `typedef shared_ptr<Signal> sig;` here.
// Unfortunately, some version of clang on MacOS was willing to
// implicitly convert a shared_ptr to an int, resulting in fatal
// ambiguity in the operator overloading.  So here we define a dummy
// class that just wraps a shared_ptr but doesn’t have a conversion to
// int.
class sig {
  shared_ptr<Signal> s;
public:
  sig(shared_ptr<Signal> s_) : s(s_) {}
  sig(Signal *s_) : s(s_) {}
  Signal &operator*() { return *s; }
  Signal *operator->() { return &*s; }
};

const int buffer_size = 64;

class Signal : public enable_shared_from_this<Signal> {
public:
  virtual float get(float time) = 0;
  virtual void get(float time, float step, float *buf) {
    float t = time;
    for (int i = 0; i != buffer_size; i++) {
      buf[i] = get(t);
      t += step;
    }
  }
  virtual bool output_graphviz_digraph(StdioGraph &out) = 0;
  virtual ~Signal() = default;

  // Arguably these methods really belong on `sig`.
  sig repeat();
  sig compose(sig timebase);
  sig speed(sig hz);
  sig speed(float hz);
  sig amplify(sig multiplier);
  sig sum(sig addend);
  sig fm(sig modulator);
  sig flange(sig modulator);
};

// Allows a reference from a `sig` to a non-dynamically-allocated
// signal.
class Leak : public Signal {
  Signal &s;
public:
  Leak(Signal &s_) : s(s_) {}
  float get(float time) {
    return s.get(time);
  }
  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, "Leak")
      && out.node_edge(this, &s)
      && s.output_graphviz_digraph(out);
  }
};

// Fundamental signal: always return a constant.
class Constant : public Signal {
  float value;
public:
  Constant(float value_) : value(value_) {}

  float get(float _) {
    return value;
  }

  void get(float time, float step, float *buf) {
    for (int i = 0; i != buffer_size; i++) {
      buf[i] = value;
    }
  }

  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, to_string(value));
  }
};

inline sig constant(float value) {
  return sig(make_shared<Constant>(value));
}

// Fundamental signal: equal to the time base.  If you consider it
// over the interval [0, 1.0] it’s one cycle of a sawtooth wave.
class Identity : public Signal {
  float get(float t) {
    return t;
  }

  void get(float t, float dt, float *buf) {
    for (int i = 0; i != buffer_size; i++) {
      buf[i] = t;
      t += dt;
    }
  }

  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, "identity");
  }
};
sig identity = sig(make_shared<Identity>());

// Generate a single triangle wave with peak-to-peak amplitude of 1.0.
// XXX this is wrong.  Its amplitude is 0.5.
class Triangle : public Signal {
  float get(float phase_fraction) {
    return phase_fraction < 0.5
      ? phase_fraction - 0.25
      : 0.75 - phase_fraction;
  }

  void get(float t, float dt, float *buf) {
    float end = t * dt * buffer_size;
    if (t < 0.5 && end < 0.5) {
      float v = t - 0.25;
      for (int i = 0; i < buffer_size; i++) {
        buf[i] = v;
        v += dt;
      }
    } else if (t >= 0.5 && end >= 0.5) {
      float v = 0.75 - t;
      for (int i = 0; i < buffer_size; i++) {
        buf[i] = v;
        v -= dt;
      }
    } else {
      Signal::get(t, dt, buf);
    }
  }

  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, "triangle");
  }
};

// Repeat a single wave.
class Repeat : public Signal {
  sig s;
public:
  Repeat(sig s_) : s(s_) {}

  float get(float time) {
    return s->get(time - floor(time));
  }

  void get(float t, float dt, float *buf) {
    float t0 = t - floor(t);
    if (t0 + dt * (buffer_size-1) > 1) {
      Signal::get(t, dt, buf);
      return;
    }
    s->get(t0, dt, buf);
  }

  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, "repeat")
      && out.node_edge(this, &*s)
      && s->output_graphviz_digraph(out);
  }
};

inline sig Signal::repeat() {
  return sig(make_shared<Repeat>(shared_from_this()));
}

sig triangle = sig(make_shared<Triangle>())->repeat();

// Add two signals.
class Sum : public Signal {
  sig a, b;
public:
  Sum(sig a_, sig b_) : a(a_), b(b_) {}

  float get(float t) {
    return a->get(t) + b->get(t);
  }

  void get(float t, float dt, float *buf) {
    // fprintf(stderr, "hello from sum\n");
    float bs[buffer_size];
    a->get(t, dt, buf);
    b->get(t, dt, bs);
    for (int i = 0; i != buffer_size; i++) {
      buf[i] += bs[i];
    }
  }

  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, "+")
      && out.node_edge(this, &*a)
      && a->output_graphviz_digraph(out)
      && out.node_edge(this, &*b)
      && b->output_graphviz_digraph(out);
    // XXX we should really avoid retraversing the same graph nodes
    // again...
  }
};

inline sig Signal::sum(sig addend) {
  return sig(make_shared<Sum>(shared_from_this(), addend));
}

static inline sig operator+(sig a, sig b) {
  return a->sum(b);
}

static inline sig operator+(sig a, float b) {
  return a + constant(b);
}

// Generate a single sawtooth with peak-to-peak amplitude of 1.0 and no DC.
sig sawtooth = (identity + -0.5)->repeat();

// Amplify a waveform to a given amplitude.  Also does AM modulation
// as a side effect.
class Amplify : public Signal {
  sig a, b;
public:
  Amplify(sig a_, sig b_) : a(a_), b(b_) {}

  float get(float t) {
    return a->get(t) * b->get(t);
  }

  void get(float t, float dt, float *buf) {
    float bs[buffer_size];
    a->get(t, dt, buf);
    b->get(t, dt, bs);
    for (int i = 0; i != buffer_size; i++) {
      buf[i] *= bs[i];
    }
  }

  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, "×")
      && out.node_edge(this, &*a)
      && a->output_graphviz_digraph(out)
      && out.node_edge(this, &*b)
      && b->output_graphviz_digraph(out);
  }
};

inline sig Signal::amplify(sig multiplier) {
  return sig(make_shared<Amplify>(shared_from_this(), multiplier));
}

static inline sig operator*(sig a, sig b) {
  return a->amplify(b);
}

static inline sig operator*(sig a, float b) {
  return a * constant(b);
}

// Use one signal as the time base for another signal; the right-hand
// side is applied to the original timebase to get the timebase for
// the left-hand side.
class Compose : public Signal {
  sig a, b;
public:
  Compose(sig a_, sig b_) :  a(a_), b(b_) {}
  float get(float t) {
    return a->get(b->get(t));
  }

  void get(float t, float dt, float *buf) {
    float modulator[buffer_size];
    b->get(t, dt, modulator);

    // Is time linear?
    float dm = modulator[1] - modulator[0];
    for (int i = 2; i != buffer_size; i++) {
      if (fabsf(modulator[i] - modulator[i-1] - dm) > fabsf(dm * 5e-2)) {
        // Time is uneven; fall back on more or less brute force.
        for (int j = 0; j != buffer_size; j++) {
          buf[j] = a->get(modulator[j]);
        }
        return;
      }
    }

    // Time is linear.
    a->get(modulator[0], dm, buf);
  }

  bool output_graphviz_digraph(StdioGraph &out) {
    return out.node_name(this, "∘")
      && out.node_edge(this, &*a, "carrier")
      && a->output_graphviz_digraph(out)
      && out.node_edge(this, &*b, "modulator")
      && b->output_graphviz_digraph(out);
  }
};

inline sig Signal::compose(sig timebase) {
  return sig(make_shared<Compose>(shared_from_this(), timebase));
}

// Adjust the time of a wave to a given number of cycles per second.
// In the version that takes a signal as an argument, note that
// although discrete changes in frequency will work correctly,
// continuous ones will also alter it with their derivative.
inline sig Signal::speed(float hz) {
  return speed(constant(hz));
}
inline sig Signal::speed(sig hz) {
  return compose(identity * hz);
}

// Modulate the time base of one signal with another.  Exactly as with
// flange below and compose above,
// the second signal must produce values that can be
// interpreted as time lags for the first signal.
inline sig Signal::fm(sig modulator) {
  return compose(identity + modulator);
}

// Flange one signal with another.  Exactly as with fm and compose, the second
// signal must produce values that can be interpreted as time lags for
// the first signal.
inline sig Signal::flange(sig modulator) {
  return shared_from_this() + fm(modulator);
}

static inline bool output_raw(sig samples) {
  float buf[buffer_size];
  for (int t = 0; t < 65536; t += buffer_size) {
    samples->get(t, 1, buf);
    char cbuf[4 * buffer_size];
    for (int i = 0; i != buffer_size; i++) {
      short sample = buf[i] + 0.5;
      cbuf[4*i    ] = cbuf[4*i + 2] = sample & 0xff;
      cbuf[4*i + 1] = cbuf[4*i + 3] = sample >> 8;
    }
    if (1 != fwrite(cbuf, sizeof(cbuf), 1, stdout)) return false;
  }
  return true;
}

int main(int argc, char **argv) {
  int sampling_rate = 48000;
  // Part of the first note’s base frequency is C2
  // I use that for the amplitude modulation, too, but I’m not
  // sure that’s right.
  float cfreq = 65.4064;
  sig c2 = triangle->speed(cfreq)
    , flanged = sawtooth->speed(2*cfreq)->flange((triangle * .002)->speed(1.5))
    , flangeam = flanged * (c2 + 0.5)
    , s2 = flangeam + c2
    , sampled = (s2 * 20000)->speed(1.0/sampling_rate)
    ;
  if (argc == 1) {
    return !output_raw(sampled);
  } else {
    StdioGraph out(stdout);
    if (!out.start()) return 1;
    if (!sampled->output_graphviz_digraph(out)) return 1;
    return !out.finish();
  }
}