sclass ContinuousOscillators {
  macro dprint { if (debug) printVars } bool debug;
    
  int channels = 1;
  double sampleRate = 48000;
  double currentSample;
  IQuerySound sound;
  IAudioSample audio;
  SlidingWindow window;
  bool renderUsingSine; // or square wave
  bool realOnly, imagOnly; // for testing
  double maxRelFreq = 0.25; // maximum oscillator frequency relative to sample rate

  long haarFeaturesChecked;

  persistable class Oscillator extends HasKey<Frequency> {
    double currentPeriodStart;
    // all recorded intensities
    new TreeMap<DoubleRange, Channels<Complex>> intensities;
    double interval; // length of period in samples
    
    *(Frequency f) {
      super(f);
      interval = frequency().interval()*sampleRate;
    }
    
    Frequency frequency() { ret key(); }
    double interval() { ret interval; }
    
    void record_impl(DoubleRange r, Channels<Complex> c) {
      intensities.put(r, c);
    }
    
    void measure() {
      DoubleRange p1 = doubleRangeWithLength(currentPeriodStart, interval);
      AudioHaarFeature haar = new(audio, p1);
      record(this, p1, haar.getComplex());
      haarFeaturesChecked += 2;
      schuleNextOscillation();
    }
    
    void initial() {
      completionsForSample.put(currentPeriodStart, this);
    }
    
    void schuleNextOscillation() {
      completionsForSample.put(currentPeriodStart = phase360(), this);
    }
    
    double phase360() { ret currentPeriodStart + interval(); }
  } // end of Oscillator
  
  TreeHasKeyMap<Frequency, Oscillator> oscillators;
  new TreeMultiMap<Double, Oscillator> completionsForSample;
  
  *() {}
  *(IAudioSample audio, double *currentSample) { setAudio(audio); }
  *(IAudioSample audio) { setAudio(audio); } // start at time 0

  void setFrequencies(Iterable<Frequency> frequencies) {
    oscillators = new TreeHasKeyMap<Frequency, Oscillator>(map(f -> new Oscillator(f), filter(frequencies,
      f -> between(f!, 1, sampleRate*maxRelFreq))));
  }
  
  void setAudio(IAudioSample audio) {
    this.audio = audio;
    window = optCast SlidingWindow(audio);
  }
  
  int nFrequencies() { ret oscillators.size(); }
  
  void startOscillators {
    for (o : oscillators) o.initial();
  }
  
  void stepTo(DoubleRange timeRange) {
    while ping (true) {
      Double t = firstKey(completionsForSample);
      if (t != null && t <= timeRange.end) {
        var list = completionsForSample.getAndClear(t);
        dprint(+t, +list);
        for (Oscillator o : list)
          o.measure();
      } else
        break;
    }
  }
  
  void record(Oscillator o, DoubleRange r, Channels<Complex> c) {
    if (o == null) ret;
    o.record_impl(r, c);
  }
  
  Frequency lowestFrequency() { ret oscillators.firstKey(); }
  Frequency highestFrequency() { ret oscillators.lastKey(); }
  
  int minWindowSize() {
    ret iceil(sampleRate*lowestFrequency().interval());
  }
  
  void makeSmallestWindow(IQuerySound sound) {
    this.sound = sound;
    setAudio(new SlidingWindow(channels, sampleRate, sound, 0, minWindowSize());
  }
  
  void makeFullWindow(IQuerySound sound, IntRange r) {
    this.sound = sound;
    setAudio(new SlidingWindow(channels, sampleRate, sound,
      r.start, r.length()));
  }
  
  double windowSize() { ret window.length(); }
  
  S stats() { ret renderVars(+audio, +haarFeaturesChecked); }
  
  // how much horizontal area have we covered
  L<Double> coverageByFreq() {
    ret mapReversed(oscillators, o -> totalLengthOfDoubleRanges(keys(o.intensities)));
  }
  
  DoubleRange windowBounds() { ret window.bounds(); }
  
  IIntegralImage imageColumn(DoubleRange timeRange) {
    int n = nFrequencies();
    double[] col = rawImageColumn(timeRange);
    double[] col2 = normalizeDoubles(col, 255);
    ret bwIntegralImageFromFunction(1, n, (x, y) -> ifloor(col2[y]));
  }
  
  double[] rawImageColumn(DoubleRange timeRange) {
    int n = nFrequencies();
    var l = reversedList(oscillators);
    double[] col = new[n];
    for y to n: {
      TreeMap<DoubleRange, Channels<Complex>> l2 = l.get(y).intensities;
      
      // find right period
      DoubleRange period = l2.floorKey(timeRange);
      if (period == null) continue;
      
      double amplitude = l2.get(period).first().abs();
      col[y] = amplitude;
    }
    
    ret col;
  }
  
  MakesBufferedImage imageForWindow(double pixelsPerSecond) {
    new L<IIntegralImage> images;
    double len = windowBounds().length()/sampleRate;
    int pixels = iround(len*pixelsPerSecond);
    for i to pixels: {
      DoubleRange range = doubleRange(
        i/pixelsPerSecond*sampleRate+windowBounds().start(),
        (i+1)/pixelsPerSecond*sampleRate+windowBounds().start());
      //print(+range);
      images.add(imageColumn(range));
    }
    
    ret mergeBWImagesHorizontally(map toBWImage(images), spacing := 0);
  }
  
  double[] toAudio(IntRange sampleRange, Iterable<Oscillator> oscillatorsToRender default oscillators) {
    double[] samples = new[l(sampleRange)];
    
    class PerFreq {
      Oscillator o;
      DoubleRange period;
      Channels<Complex> intensity;

      *(Oscillator *o) {
        setPeriod(o.intensities.floorKey(toDoubleRange(sampleRange)));
        //printVars("PerFreq", f := o.frequency(), +sampleRange, +period, firstPeriod := firstKey(o.intensities));
      }
      
      void setPeriod(DoubleRange period) {
        this.period = period;
        intensity = period == null ?: o.intensities.get(period);
      }
      
      void nextPeriod(double t) {
        while (period != null && t >= period.end)
          setPeriod(o.intensities.higherKey(period));
      }
    }
    
    L<PerFreq> perFreqs = map(oscillatorsToRender, o -> new PerFreq(o));
    
    for i over samples: {
      double t = sampleRange.start+i; // t is time in samples
      double sum = 0;
      for (pf : perFreqs) {
        pf.nextPeriod(t);
        bool shouldPrint = toAudio_shouldPrint(i, l(samples));
        
        if (pf.intensity != null) {
          Complex c = div(pf.intensity.first(), l(pf.period));
          double t2 = t-pf.period.start;
          double frac = t2/pf.o.interval;
          //double phiRE = frac*pi()*2;
          //double phiIM = phiRE+pi()/2;
          if (shouldPrint)
            printVars(f := pf.o.frequency(), +i, interval := formatDouble1(pf.o.interval)/*, interval2 := l(pf.period)*/, +t, +t2, frac := formatDouble2(frac),
              c := renderComplexWithAngle(div(c, 32768)));
          if (realOnly) c = complex(c.re());
          else if (imagOnly) c = complex(0, c.im());
          sum += renderFrequencySample(c, frac);
        } else
          if (shouldPrint)
            printVars(f := pf.o.frequency(), +i, interval := pf.o.interval, +t, c := 0);
        
      }
      
      samples[i] = sum;
    }
  
    ret samples;
  }

  double renderFrequencySample(Complex c, double frac) {
    ret (renderUsingSine
      ? new RenderFrequencySample_Sine
      : new RenderFrequencySample_SquareWave).get(c, frac);
  }
  
  swappable bool toAudio_shouldPrint(int i, int n) { false; }
}

Began life as a copy of #1032970