Resonant sine Teensy Audio

readme.md

Distortion is 95 dB THD or so, has fewer frequency components than linear-interpolated LUT, but overall THD level is similar.

Can only go up to fs/6, though, because signed a can only go to Q31 1.0, not 2.0.

Found in http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.84.1650 but didn't work until I combined with http://www.musicdsp.org/showArchiveComment.php?ArchiveID=10

This is even faster, but distorted at low levels: https://gist.github.com/endolith/14bbb3217f9f58248722

Crude timing results (Should have measured with many instances in parallel instead):

1. Muted: 0.57%
2. arm_sin_q31: 5.77%
3. AudioSynthWaveformSine: 1.85%
4. Resonant sine: 1.67%
5. Resonant with SMMLAR: 1.63%
6. Resonant with SMMLAR/SMMLSR: 1.58% (this)
7. Resonant with SMMLAR/SMMLSR and __attribute__((optimize("unroll-loops"))): 1.53%
8. Resonant with SMMLAR/SMMLSR manually unrolled 8x: 1.41% (this)

Then compared compilers and magnitude-scaled version:

1. Resonant SMMLSR 3 mult in Arduino 26,836 bytes: 1.58% (this)
2. Resonant 2 mult in Arduino 28,340 bytes: 1.4%
3. Resonant SMMLSR 3 mult in UECIDE 21,924 bytes: 1.37%
4. Resonant 2 mult in UECIDE 23,432 bytes: 1.24%
5. Resonant 2 mult in UECIDE but unrolled 8x: 23,504 bytes, 0.94%
6. Quadratic sine first attempt: 1.98%

So best case is (1.85-0.57)/(0.94-0.57) = 2 to 3.4x improvement over linear-interpolated LUT? But doesn't have same frequency range.

resonant.cpp

#include "resonant.h"
#include "utility/dspinst.h"

void AudioSynthWaveformRes::update(void)
{
  audio_block_t *block;
  uint32_t i;

  if (magnitude) {
    block = allocate();
    if (block) {
      for (i = 0; i < AUDIO_BLOCK_SAMPLES; i++) {
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i] = multiply_32x32_rshift32_rounded(s, magnitude);
      }
      transmit(block);
      release(block);
      return;
    }
  }
  // clean up phase/amplitude drift each block? (though it doesn't drift in my tests)
}

resonant.h

#ifndef resonant_h_
#define resonant_h_

#include "AudioStream.h"
#include "arm_math.h"

class AudioSynthWaveformRes : public AudioStream
{
  public:
    AudioSynthWaveformRes() : AudioStream(0, NULL), magnitude(16384) {}
    void frequency(float freq) {
      if (freq < 0.0) freq = 0.0;
      else if (freq > AUDIO_SAMPLE_RATE_EXACT / 2) freq = AUDIO_SAMPLE_RATE_EXACT / 2;
      a = 4.0 * sin(M_PI * freq / AUDIO_SAMPLE_RATE_EXACT) * ((uint)1 << 31);
    }
    void phase(float angle) {
      if (angle < 0.0) angle = 0.0;
      else if (angle > 360.0) {
        angle = angle - 360.0;
        if (angle >= 360.0) return;
      }
      c = cos(angle * M_PI / 180) * (1 << 30);
      s = sin(angle * M_PI / 180) * (1 << 30);
    }
    void amplitude(float n) {
      if (n < 0) n = 0;
      else if (n > 1.0) n = 1.0;
      magnitude = n * 131072;//65536.0;
    }
    virtual void update(void);
  private:
    int32_t a;
    int32_t c = 1 << 30;
    int32_t s = 0;
    int32_t magnitude;
};

#endif

resonant_(unrolled).cpp

#include "resonant.h"
#include "utility/dspinst.h"

void AudioSynthWaveformRes::update(void)
{
  audio_block_t *block;
  uint32_t i;

  if (magnitude) {
    block = allocate();
    if (block) {
      for (i = 0; i < AUDIO_BLOCK_SAMPLES; i+=8) {
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i] = multiply_32x32_rshift32(s, magnitude);
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i+1] = multiply_32x32_rshift32(s, magnitude);
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i+2] = multiply_32x32_rshift32(s, magnitude);
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i+3] = multiply_32x32_rshift32(s, magnitude);        
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i+4] = multiply_32x32_rshift32(s, magnitude);
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i+5] = multiply_32x32_rshift32(s, magnitude);
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i+6] = multiply_32x32_rshift32(s, magnitude);
        c = multiply_subtract_32x32_rshift32_rounded(c, a, s);
        s = multiply_accumulate_32x32_rshift32_rounded(s, a, c);
        block->data[i+7] = multiply_32x32_rshift32(s, magnitude);         
      }
      transmit(block);
      release(block);
      return;
    }
  }
  // clean up phase/amplitude drift each block? (though it doesn't drift in my tests)
}

resonant_sine.ino

#include <Audio.h>
#include <Wire.h>
#include <SD.h>
#include <SPI.h>

#include <math.h>
#include "resonant.h"

AudioSynthWaveformRes    waveform1;      //xy=188,240
AudioOutputI2S           i2s1;           //xy=565,241
AudioConnection          patchCord2(waveform1, 0, i2s1, 0);
AudioConnection          patchCord3(waveform1, 0, i2s1, 1);
AudioControlSGTL5000     audioShield;     //xy=586,175

void setup() {
  Serial.begin(9600);

  AudioMemory(80);
  audioShield.enable();
  audioShield.volume(0.3);

  waveform1.amplitude(0.8);  // 1.26 Vpp for 1.0
  waveform1.frequency(23.45678); // can only go up to 3.5 kHz because of range of a
}

void loop() {
  Serial.print(AudioProcessorUsage());
  Serial.print(",");
  Serial.println(AudioProcessorUsageMax());
  delay(500);
}