<reponame>drepetto/alles<filename>analysis/2021_04_29_blosca_simulation.py<gh_stars>10-100
# -*- coding: utf-8 -*-
"""2021-04-29-blosca_simulation.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1PSBC98l3rYV6GIQkevNhWRUU4SLY7CtA

# 2021-04-29-blosca_simulation.ipynb

Simulate the process of synthesizing waveforms by lookup of bandlimited impulse tables.  This is to prototype the C/C++ implementation in http://github.com/bwhitman/alles.
"""

import collections

import numpy as np
import matplotlib.pyplot as plt

def plotspec(waveform, samplerate=44100, window_fn=np.hanning):
  fftsize = len(waveform)
  windowlength = fftsize
  window = window_fn(windowlength)
  wavepart = waveform[:len(window)]
  logspecmag = 20 * np.log10(
      np.maximum(1e-10, 
                 np.abs(np.fft.fft(wavepart * window)))[:(fftsize // 2 + 1)])
  freqs = samplerate * np.arange(len(logspecmag)) / fftsize
  plt.subplot(211)
  times = np.arange(len(wavepart)) / samplerate
  plt.plot(times, wavepart, '.')
  plt.subplot(212)
  plt.plot(freqs, logspecmag, '.-')
  plt.ylim(np.array([-100, 0]) + np.max(logspecmag))

x = np.zeros((1024,))
x[512] = 1
x[519] = 1
plotspec(x)

def interp(table, fractions, order=1):
  table_len = len(table)
  samples = fractions * table_len
  sample_base = np.floor(samples).astype(int)
  sample_frac = samples - sample_base
  a = table[(sample_base - 1) % table_len]
  b = table[(sample_base + 0) % table_len]
  c = table[(sample_base + 1) % table_len]
  d = table[(sample_base + 2) % table_len]
  # Linear interp
  if order == 1:
    vals = b + sample_frac * (c - b)
  elif order == 3:
    vals = (- sample_frac * (sample_frac - 1) * (sample_frac - 2) / 6.0 * a
            + (sample_frac + 1) * (sample_frac - 1) * (sample_frac - 2) / 2.0 * b
            - (sample_frac + 1) * sample_frac * (sample_frac - 2) / 2.0 * c
            + (sample_frac + 1) * sample_frac * (sample_frac - 1) / 6.0 * d)
  else:
    raise ValueError("order of {} is not 1 or 3".format(order))
  return vals

def cos_lut(table_size, harmonics_weights, harmonics_phases=None):
  if harmonics_phases is None:
    harmonics_phases = np.zeros(len(harmonics_weights))
  table = np.zeros(table_size)
  phases = np.arange(table_size) * 2 * np.pi / table_size
  for harmonic_number, harmonic_weight in enumerate(harmonics_weights):
    table += harmonic_weight * np.cos(phases * harmonic_number + harmonics_phases[harmonic_number])
  return table

cos_32 = cos_lut(32, [0, 0, 0, 0, 1])
plotspec(interp(cos_32, np.linspace(0, 1, num=1024, endpoint=False)))

plotspec(interp(cos_32, np.linspace(0, 1, num=1024, endpoint=False), order=3))

plotspec(interp(cos_lut(256, [0, 0, 1, 0, 0]), np.linspace(0, 1, num=1024, endpoint=False), order=1))

plotspec(interp(cos_lut(256, [0, 0, 1, 0, 0], -np.pi/2*np.ones(5)), np.linspace(0, 1, num=1024, endpoint=False), order=1))

# Centered
#imp16 = cos_lut(256, np.ones(16), np.arange(16)*-np.pi) / 256
# Left-edge
imp16 = cos_lut(256, np.ones(16)) / 256
plotspec(imp16, window_fn=np.ones)

print(np.sum(imp16))
print(np.max(imp16)*256)
print(np.sum(imp16 * np.cos(np.arange(256)*2*np.pi/256)))

def integrate(signal, alpha, state=0):
  out = np.zeros(signal.shape)
  for i in range(len(out)):
    out[i] = signal[i] + alpha * state
    state = out[i]
  return out

def square_wave(period, length, initial_phase=0, initial_state=0):
  """Synthesize <length> samples of a square wave of <period> by integrating the bandlimited impulse."""
  out = np.zeros(length)
  # Phase is in cycles (1.0 = 1 cycle), i.e. radians/(2pi)
  normphase = initial_phase + (np.arange(length) / period)
  # Actual step-size within LUT.  If <period> samples in output corresponds to len(LUT), 
  # then each sample of output corresponds to a step of len(LUT)/period LUT samples.
  step = len(imp16) / period
  print("step=",step)
  # Because stepping by <step> reduces the effective total number of samples in the LUT
  # by 1/step, and thus the total integral is reduced by 1/step, we scale the impulses 
  # by step so their integral stays constant.
  impscale = 2 * step
  impulses = impscale * (interp(imp16, normphase) - interp(imp16, normphase + 0.5))
  # Smoothing time is 100x period for little droop.
  alpha = 1 - 1 / (1000 * period)
  # To get the exact balance from the first sample, we have to discount half the first value of the
  # bandlimited impulse, because to start perfectly in the middle of the interpolate impulse, we 
  # would start "half-way through" the zero'th sample.
  # The summation is a little off because it isn't exactly intg(imp), but has the decay forgetting.
  # Making the forgetting large makes it very close to the correct preset.
  initial_state -= 0.5 * (impscale * imp16[0])
  return integrate(impulses, alpha=alpha, state=initial_state)

sq = square_wave(32, 4096) #, initial_phase=0.25, initial_state=0.5)

plotspec(sq[:64])

print(32*np.sum(sq[:32]))
print(np.sum(sq[32:64]))
print(np.sum(sq[512:768]))
print(np.sum(sq[768:1024]))

def saw_wave(period, length, initial_phase=0, initial_state=0):
  """Synthesize <length> samples of a sawtooth of <period> by integrating the bandlimited impulse."""
  out = np.zeros(length)
  # Phase is in cycles (1.0 = 1 cycle), i.e. radians/(2pi)
  normphase = initial_phase + (np.arange(length) / period)
  # Actual step-size within LUT.  If <period> samples in output corresponds to len(LUT), 
  # then each sample of output corresponds to a step of len(LUT)/period LUT samples.
  step = len(imp16) / period
  print("step=",step)
  # Because stepping by <step> reduces the effective total number of samples in the LUT
  # by 1/step, and thus the total integral is reduced by 1/step, we scale the impulses 
  # by step so their integral stays constant.
  impscale = 2 * step
  impulses = impscale * interp(imp16, normphase)
  # Smoothing time is 100x period for little droop.
  alpha = 1 - 1 / (10 * period)
  # To get the exact balance from the first sample, we have to discount half the first value of the
  # bandlimited impulse, because to start perfectly in the middle of the interpolate impulse, we 
  # would start "half-way through" the zero'th sample.
  # The summation is a little off because it isn't exactly intg(imp), but has the decay forgetting.
  # Making the forgetting large makes it very close to the correct preset.
  initial_state = -0.5 * (impscale * imp16[0])
  #initial_state = 0
  #return impulses
  return integrate(impulses - np.mean(impscale * imp16), alpha=alpha, state=initial_state)

print(sum(imp16))

saw = saw_wave(200, 4096) #, initial_phase=0.25, initial_state=0.5)

plotspec(saw)

def triangle_wave(period, length):
  # The square wave is at +0.5 for period/2, so its integral would be period/4.
  # But we start it 1/4 through its cycle, so we can start the triangle at 0 
  # (midway through the upswing).
  sq_wave = square_wave(period, length, initial_phase=0.25, initial_state=0.5)
  # 500; + 0.012; 400; + 0.021; 300: + 0.034; 200: + 0.061; 100 + 0.143; 50: 0.307
  #return sq_wave
  alpha = 1 - 1 / (10 * period)
  # We start at zero.
  return integrate((4/period) * sq_wave, alpha)

tri = triangle_wave(400, 4096)

plotspec(tri)#[:112])

print(tri[:6])
print(tri[512:518])
print(tri[2048:2054])

print(np.sum(tri[:256]))
plt.plot(np.arange(256), tri[:256], np.arange(256), -tri[256:512])

clut = cos_lut(256, np.ones(16))   # 256 * imp16
period = 256
clint = interp(clut, np.arange(period)/period)
sqint = clint - np.hstack([clint[len(clint)//2:], clint[:len(clint)//2]])
plt.plot(sqint)
print(len(clut), len(clint), np.sum(clut), np.sum(clint))
print(np.sum(sqint[:(period//4)]), np.sum(sqint[:(period//2)])) #, np.sum(sqint[:(3*period//4))]), np.sum(sqint[:period]))
print(np.sum(sqint[1:129]), np.sum(sqint[1:257]), np.sum(sqint[1:385]), np.sum(sqint[:512]))

(131 - 128)/256

# Make a separate LUT for triangle.
n_harms = 64
coefs = (np.arange(n_harms) % 2)*(np.maximum(1, np.arange(n_harms, dtype=float))**(-2))
#print(coefs)
triLUT = cos_lut(256, coefs, np.arange(n_harms)*-np.pi/2)
plotspec(triLUT, window_fn=np.ones)

plt.plot(triLUT[64 - 16 + np.arange(32)], '.')

# Implement the multiple lookup tables.
# A LUT is stored as an array of values (table) and the harmonic number of the
# highest harmonic they contain (i.e., the number of cycles it completes in the
# entire table, so must be <= len(table)/2.)
LUTentry = collections.namedtuple('LUTentry', ['table', 'highest_harmonic'])

# A LUTset is a list of LUTentries describing downsampled versions of the same
# basic waveform, sorted with the longest (highest-bandwidth) first.

def create_lutset(harmonic_weights, harmonic_phases=None, 
                  length_factor=8, bandwidth_factor=np.sqrt(0.5)):
  """Create an ordered list of LUTs with decreasing harmonic content.

  These can then be used in interp_from_lutset to make an adaptive-bandwidth
  interpolation.

  Args:
    harmonic_weights: vector of amplitudes for cosine harmonic components.
    harmonic_phases: initial phases for each harmonic, in radians.  Zero 
      (default) indicates cosine phase.
    length_factor: Each table's length is at least this factor times the order
      of the highest harmonic it contains.  Thus, this is a lower bound on the
      number of samples per cycle for the highest harmonic.  Higher factors make
      the interpolation easier.
    bandwidth_factor: Target ratio between the highest harmonics in successive
      table entries.  Default is sqrt(0.5), so after two tables, bandwidth is
      reduced by 1/2 (and length with follow).

  Returns:
    A list of LUTentry objects, sorted in decreasing order of the highest 
    harmonic they contain.  Each LUT's length is a power of 2, and as small as
    possible while respecting the length_factor for the highest contained 
    harmonic.
  """
  if harmonic_phases is None:
    harmonic_phases = np.zeros(len(harmonic_weights))
  # Calculate the length of the longest LUT we need.  Must be a power of 2, 
  # must have at least length_factor * highest_harmonic samples.
  # Harmonic 0 (dc) doesn't count.
  float_num_harmonics = float(len(harmonic_weights))
  lutsets = []
  done = False
  # harmonic 0 is DC; there's no point in generating that table.
  while float_num_harmonics >= 2:
    num_harmonics = int(round(float_num_harmonics))
    highest_harmonic = num_harmonics - 1  # because zero doesn't count.
    lut_size = int(2 ** np.ceil(np.log(length_factor * highest_harmonic) /
                                np.log(2)))
    print(float_num_harmonics, num_harmonics, lut_size)
    lutsets.append(LUTentry(
        table=cos_lut(lut_size, 
                      harmonic_weights[:num_harmonics], 
                      harmonic_phases[:num_harmonics]),  # / lut_size,
        highest_harmonic=highest_harmonic))
    float_num_harmonics = bandwidth_factor * float_num_harmonics
  return lutsets

# Centered
#lutset = create_lutset(np.ones(128), np.arange(128) * -np.pi)
# Left-aligned
lutset = create_lutset(np.ones(128))
print([(len(s.table), s.highest_harmonic) for s in lutset])

plt.figure(figsize=(12, 6))
lut = lutset[0]
plotspec(lut.table, window_fn=np.ones)
print(np.sum(lutset[0].table))
print(np.sum(lutset[1].table))
print(np.sum(lutset[2].table))

def choose_from_lutset(period, lutset):
  """Select the best entry from a lutset for a given period.

  Args:
    period: (float) Target period of waveform in fractional samples.
    lutset: Sorted list of LUTs, as generated by create_lutset().
  
  Returns:
    One of the LUTs from the lutset, best suited to interpolating to generate
    a waveform of the desired period.
  """
  # Use the earliest (i.e., longest, most harmonics) LUT that works
  # (i.e., will not actually cause aliasing).
  # So start with the highest-bandwidth (and longest) LUTs, but skip them
  # if they result in aliasing.
  for lut_index, lutentry in enumerate(lutset):
    lut_table = lutentry.table
    lut_size = len(lut_table)
    # What proportion of nyquist does the highest harmonic in this table occupy?
    lut_bandwidth = 2 * lutentry.highest_harmonic / lut_size
    # To complete one cycle of <lut_size> points in <period> steps, each step
    # will need to be this many samples:
    lut_hop = lut_size / period
    # If we have a signal with a given bandwidth, but then speed it up by 
    # skipping lut_hop samples per sample, its bandwidth will increase 
    # proportionately.
    interp_bandwidth = lut_bandwidth * lut_hop
    if interp_bandwidth < 1.0:
      # No aliasing
      break
  # At this point, we either got to the end of the LUT table, or we found a
  # table we could interpolate without aliasing.
  print("lut_size=", lut_size, "lut_hop=", lut_hop, "lut_ix=", lut_index, "lut_bw=", lut_bandwidth)
  return lut_table

def interp_from_lutset(period, num_samples, lutset):
  """Generate a stretch of interpolated waveform from a set of LUTs.

  Args:
    period: (float) Target period of waveform in fractional samples.
    num_samples: (int) Total number of samples to generate at that period.
    lutset: Sorted list of LUTs, as generated by create_lutset().
  
  Returns:
    Vector of num_samples samples containing interpolated waveform with
    requested period.
  """
  lut_table = choose_from_lutset(period, lutset)
  lut_size = len(lut_table)
  # To complete one cycle of <lut_size> points in <period> steps, each step
  # will need to be this many samples:
  lut_hop = lut_size / period
  return interp(lut_table, np.arange(num_samples) * lut_hop / lut_size, order=3)

plotspec(interp_from_lutset(70, 210, lutset), window_fn=np.ones)

plotspec(interp_from_lutset(45, 135, lutset), window_fn=np.ones)

plotspec(interp_from_lutset(44, 132, lutset), window_fn=np.ones)

# Saving out the lutset as a C-compatible header file.

def write_lutset_to_h(filename, variable_base, lutset):
  num_luts = len(lutset)
  with open(filename, "w") as f:
    f.write("// Automatically-generated LUTset\n")
    f.write("#ifndef LUTSET_{:s}_DEFINED\n".format(variable_base.upper()))
    f.write("#define LUTSET_{:s}_DEFINED\n".format(variable_base.upper()))
    f.write("\n")
    # Define the structure.
    f.write("#ifndef LUTENTRY_DEFINED\n")
    f.write("#define LUTENTRY_DEFINED\n")
    f.write("typedef struct {\n")
    f.write("  const float *table;\n")
    f.write("  int table_size;\n")
    f.write("  int highest_harmonic;\n")
    f.write("} lut_entry;\n")
    f.write("#endif // LUTENTRY_DEFINED\n")
    f.write("\n")
    # Define the content of the individual tables.
    samples_per_row = 8
    for i in range(num_luts):
      table_size = len(lutset[i].table)
      f.write("const float {:s}_lutable_{:d}[{:d}] = {{\n".format(
          variable_base, i, table_size))
      for row_start in range(0, table_size, samples_per_row):
        for sample_index in range(row_start, 
                                  min(row_start + samples_per_row, table_size)):
          f.write("{:f},".format(lutset[i].table[sample_index]))
        f.write("\n")
      f.write("};\n")
      f.write("\n")
    # Define the table of LUTs.
    f.write("lut_entry {:s}_lutset[{:d}] = {{\n".format(
        variable_base, num_luts + 1))
    for i in range(num_luts):
      f.write("  {{{:s}_lutable_{:d}, {:d}, {:d}}},\n".format(
          variable_base, i, len(lutset[i].table), 
          lutset[i].highest_harmonic))
    # Final entry is null to indicate end of table.
    f.write("  {NULL, 0, 0},\n")
    f.write("};\n")
    f.write("\n")
    f.write("#endif // LUTSET_x_DEFINED\n")
  print("wrote", filename)

write_lutset_to_h('impulse_lutset.h', 'impulse', lutset)

!cat impulse_lutset.h

# Triangle wave lutset
n_harms = 64
coefs = (np.arange(n_harms) % 2)*(np.maximum(1, np.arange(n_harms, dtype=float))**(-2))
triangle_lutset = create_lutset(coefs, np.arange(len(coefs)) * -np.pi/2)

print([(len(s.table), s.highest_harmonic) for s in triangle_lutset])
plotspec(triangle_lutset[-3].table, window_fn=np.ones)

write_lutset_to_h('triangle_lutset.h', 'triangle', triangle_lutset)
!cat triangle_lutset.h

# Sinusoid "lutset" (only one table)
sine_lutset = create_lutset(np.array([0, 1]), length_factor=256)

print([(len(s.table), s.highest_harmonic) for s in sine_lutset])
plotspec(sine_lutset[0].table, window_fn=np.ones)

write_lutset_to_h('sine_lutset.h', 'sine', sine_lutset)
!cat sine_lutset.h

print(0.9*32768, 0.15 * 32768)

# Soft clipping lookup table scratchpad.

INT16_MAX = 32767
LINMAX = 29491  #// int(round(0.9 * 32768))
NONLIN_RANGE = 4915  # // size of nonlinearity lookup table = round(1.5 * (INT16_MAX - LINMAX))

clipping_lookup_table = np.arange(LINMAX + NONLIN_RANGE)
print(len(clipping_lookup_table))

for x in range(NONLIN_RANGE):
  x_dash = float(x) / NONLIN_RANGE
  clipping_lookup_table[x + LINMAX] = LINMAX + int(np.floor(NONLIN_RANGE * (x_dash - x_dash * x_dash * x_dash / 3.0)))

plt.plot(clipping_lookup_table[25000:]) # - np.arange(25000, len(clipping_lookup_table)))
print(max(clipping_lookup_table))

filename = "clipping_lookup_table.h"

with open(filename, "w") as f:
  f.write("// Automatically generated.\n// Clipping lookup table\n")
  f.write("const uint16 clipping_lookup_table[%d] = {\n" % NONLIN_RANGE)
  samples_per_row = 8
  for row_start in range(0, NONLIN_RANGE, samples_per_row):
    for sample in range(row_start, min(NONLIN_RANGE, row_start + samples_per_row)):
      f.write("%d," % clipping_lookup_table[LINMAX + sample])
    f.write("\n")
  f.write("};\n")
print("wrote", filename)

!cat clipping_lookup_table.h


# -*- coding: utf-8 -*-
"""2021-04-29-blosca_simulation.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1PSBC98l3rYV6GIQkevNhWRUU4SLY7CtA

# 2021-04-29-blosca_simulation.ipynb

Simulate the process of synthesizing waveforms by lookup of bandlimited impulse tables.  This is to prototype the C/C++ implementation in http://github.com/bwhitman/alles.
"""

import collections

import numpy as np
import matplotlib.pyplot as plt

def plotspec(waveform, samplerate=44100, window_fn=np.hanning):
  fftsize = len(waveform)
  windowlength = fftsize
  window = window_fn(windowlength)
  wavepart = waveform[:len(window)]
  logspecmag = 20 * np.log10(
      np.maximum(1e-10, 
                 np.abs(np.fft.fft(wavepart * window)))[:(fftsize // 2 + 1)])
  freqs = samplerate * np.arange(len(logspecmag)) / fftsize
  plt.subplot(211)
  times = np.arange(len(wavepart)) / samplerate
  plt.plot(times, wavepart, '.')
  plt.subplot(212)
  plt.plot(freqs, logspecmag, '.-')
  plt.ylim(np.array([-100, 0]) + np.max(logspecmag))

x = np.zeros((1024,))
x[512] = 1
x[519] = 1
plotspec(x)

def interp(table, fractions, order=1):
  table_len = len(table)
  samples = fractions * table_len
  sample_base = np.floor(samples).astype(int)
  sample_frac = samples - sample_base
  a = table[(sample_base - 1) % table_len]
  b = table[(sample_base + 0) % table_len]
  c = table[(sample_base + 1) % table_len]
  d = table[(sample_base + 2) % table_len]
  # Linear interp
  if order == 1:
    vals = b + sample_frac * (c - b)
  elif order == 3:
    vals = (- sample_frac * (sample_frac - 1) * (sample_frac - 2) / 6.0 * a
            + (sample_frac + 1) * (sample_frac - 1) * (sample_frac - 2) / 2.0 * b
            - (sample_frac + 1) * sample_frac * (sample_frac - 2) / 2.0 * c
            + (sample_frac + 1) * sample_frac * (sample_frac - 1) / 6.0 * d)
  else:
    raise ValueError("order of {} is not 1 or 3".format(order))
  return vals

def cos_lut(table_size, harmonics_weights, harmonics_phases=None):
  if harmonics_phases is None:
    harmonics_phases = np.zeros(len(harmonics_weights))
  table = np.zeros(table_size)
  phases = np.arange(table_size) * 2 * np.pi / table_size
  for harmonic_number, harmonic_weight in enumerate(harmonics_weights):
    table += harmonic_weight * np.cos(phases * harmonic_number + harmonics_phases[harmonic_number])
  return table

cos_32 = cos_lut(32, [0, 0, 0, 0, 1])
plotspec(interp(cos_32, np.linspace(0, 1, num=1024, endpoint=False)))

plotspec(interp(cos_32, np.linspace(0, 1, num=1024, endpoint=False), order=3))

plotspec(interp(cos_lut(256, [0, 0, 1, 0, 0]), np.linspace(0, 1, num=1024, endpoint=False), order=1))

plotspec(interp(cos_lut(256, [0, 0, 1, 0, 0], -np.pi/2*np.ones(5)), np.linspace(0, 1, num=1024, endpoint=False), order=1))

# Centered
#imp16 = cos_lut(256, np.ones(16), np.arange(16)*-np.pi) / 256
# Left-edge
imp16 = cos_lut(256, np.ones(16)) / 256
plotspec(imp16, window_fn=np.ones)

print(np.sum(imp16))
print(np.max(imp16)*256)
print(np.sum(imp16 * np.cos(np.arange(256)*2*np.pi/256)))

def integrate(signal, alpha, state=0):
  out = np.zeros(signal.shape)
  for i in range(len(out)):
    out[i] = signal[i] + alpha * state
    state = out[i]
  return out

def square_wave(period, length, initial_phase=0, initial_state=0):
  """Synthesize <length> samples of a square wave of <period> by integrating the bandlimited impulse."""
  out = np.zeros(length)
  # Phase is in cycles (1.0 = 1 cycle), i.e. radians/(2pi)
  normphase = initial_phase + (np.arange(length) / period)
  # Actual step-size within LUT.  If <period> samples in output corresponds to len(LUT), 
  # then each sample of output corresponds to a step of len(LUT)/period LUT samples.
  step = len(imp16) / period
  print("step=",step)
  # Because stepping by <step> reduces the effective total number of samples in the LUT
  # by 1/step, and thus the total integral is reduced by 1/step, we scale the impulses 
  # by step so their integral stays constant.
  impscale = 2 * step
  impulses = impscale * (interp(imp16, normphase) - interp(imp16, normphase + 0.5))
  # Smoothing time is 100x period for little droop.
  alpha = 1 - 1 / (1000 * period)
  # To get the exact balance from the first sample, we have to discount half the first value of the
  # bandlimited impulse, because to start perfectly in the middle of the interpolate impulse, we 
  # would start "half-way through" the zero'th sample.
  # The summation is a little off because it isn't exactly intg(imp), but has the decay forgetting.
  # Making the forgetting large makes it very close to the correct preset.
  initial_state -= 0.5 * (impscale * imp16[0])
  return integrate(impulses, alpha=alpha, state=initial_state)

sq = square_wave(32, 4096) #, initial_phase=0.25, initial_state=0.5)

plotspec(sq[:64])

print(32*np.sum(sq[:32]))
print(np.sum(sq[32:64]))
print(np.sum(sq[512:768]))
print(np.sum(sq[768:1024]))

def saw_wave(period, length, initial_phase=0, initial_state=0):
  """Synthesize <length> samples of a sawtooth of <period> by integrating the bandlimited impulse."""
  out = np.zeros(length)
  # Phase is in cycles (1.0 = 1 cycle), i.e. radians/(2pi)
  normphase = initial_phase + (np.arange(length) / period)
  # Actual step-size within LUT.  If <period> samples in output corresponds to len(LUT), 
  # then each sample of output corresponds to a step of len(LUT)/period LUT samples.
  step = len(imp16) / period
  print("step=",step)
  # Because stepping by <step> reduces the effective total number of samples in the LUT
  # by 1/step, and thus the total integral is reduced by 1/step, we scale the impulses 
  # by step so their integral stays constant.
  impscale = 2 * step
  impulses = impscale * interp(imp16, normphase)
  # Smoothing time is 100x period for little droop.
  alpha = 1 - 1 / (10 * period)
  # To get the exact balance from the first sample, we have to discount half the first value of the
  # bandlimited impulse, because to start perfectly in the middle of the interpolate impulse, we 
  # would start "half-way through" the zero'th sample.
  # The summation is a little off because it isn't exactly intg(imp), but has the decay forgetting.
  # Making the forgetting large makes it very close to the correct preset.
  initial_state = -0.5 * (impscale * imp16[0])
  #initial_state = 0
  #return impulses
  return integrate(impulses - np.mean(impscale * imp16), alpha=alpha, state=initial_state)

print(sum(imp16))

saw = saw_wave(200, 4096) #, initial_phase=0.25, initial_state=0.5)

plotspec(saw)

def triangle_wave(period, length):
  # The square wave is at +0.5 for period/2, so its integral would be period/4.
  # But we start it 1/4 through its cycle, so we can start the triangle at 0 
  # (midway through the upswing).
  sq_wave = square_wave(period, length, initial_phase=0.25, initial_state=0.5)
  # 500; + 0.012; 400; + 0.021; 300: + 0.034; 200: + 0.061; 100 + 0.143; 50: 0.307
  #return sq_wave
  alpha = 1 - 1 / (10 * period)
  # We start at zero.
  return integrate((4/period) * sq_wave, alpha)

tri = triangle_wave(400, 4096)

plotspec(tri)#[:112])

print(tri[:6])
print(tri[512:518])
print(tri[2048:2054])

print(np.sum(tri[:256]))
plt.plot(np.arange(256), tri[:256], np.arange(256), -tri[256:512])

clut = cos_lut(256, np.ones(16))   # 256 * imp16
period = 256
clint = interp(clut, np.arange(period)/period)
sqint = clint - np.hstack([clint[len(clint)//2:], clint[:len(clint)//2]])
plt.plot(sqint)
print(len(clut), len(clint), np.sum(clut), np.sum(clint))
print(np.sum(sqint[:(period//4)]), np.sum(sqint[:(period//2)])) #, np.sum(sqint[:(3*period//4))]), np.sum(sqint[:period]))
print(np.sum(sqint[1:129]), np.sum(sqint[1:257]), np.sum(sqint[1:385]), np.sum(sqint[:512]))

(131 - 128)/256

# Make a separate LUT for triangle.
n_harms = 64
coefs = (np.arange(n_harms) % 2)*(np.maximum(1, np.arange(n_harms, dtype=float))**(-2))
#print(coefs)
triLUT = cos_lut(256, coefs, np.arange(n_harms)*-np.pi/2)
plotspec(triLUT, window_fn=np.ones)

plt.plot(triLUT[64 - 16 + np.arange(32)], '.')

# Implement the multiple lookup tables.
# A LUT is stored as an array of values (table) and the harmonic number of the
# highest harmonic they contain (i.e., the number of cycles it completes in the
# entire table, so must be <= len(table)/2.)
LUTentry = collections.namedtuple('LUTentry', ['table', 'highest_harmonic'])

# A LUTset is a list of LUTentries describing downsampled versions of the same
# basic waveform, sorted with the longest (highest-bandwidth) first.

def create_lutset(harmonic_weights, harmonic_phases=None, 
                  length_factor=8, bandwidth_factor=np.sqrt(0.5)):
  """Create an ordered list of LUTs with decreasing harmonic content.

  These can then be used in interp_from_lutset to make an adaptive-bandwidth
  interpolation.

  Args:
    harmonic_weights: vector of amplitudes for cosine harmonic components.
    harmonic_phases: initial phases for each harmonic, in radians.  Zero 
      (default) indicates cosine phase.
    length_factor: Each table's length is at least this factor times the order
      of the highest harmonic it contains.  Thus, this is a lower bound on the
      number of samples per cycle for the highest harmonic.  Higher factors make
      the interpolation easier.
    bandwidth_factor: Target ratio between the highest harmonics in successive
      table entries.  Default is sqrt(0.5), so after two tables, bandwidth is
      reduced by 1/2 (and length with follow).

  Returns:
    A list of LUTentry objects, sorted in decreasing order of the highest 
    harmonic they contain.  Each LUT's length is a power of 2, and as small as
    possible while respecting the length_factor for the highest contained 
    harmonic.
  """
  if harmonic_phases is None:
    harmonic_phases = np.zeros(len(harmonic_weights))
  # Calculate the length of the longest LUT we need.  Must be a power of 2, 
  # must have at least length_factor * highest_harmonic samples.
  # Harmonic 0 (dc) doesn't count.
  float_num_harmonics = float(len(harmonic_weights))
  lutsets = []
  done = False
  # harmonic 0 is DC; there's no point in generating that table.
  while float_num_harmonics >= 2:
    num_harmonics = int(round(float_num_harmonics))
    highest_harmonic = num_harmonics - 1  # because zero doesn't count.
    lut_size = int(2 ** np.ceil(np.log(length_factor * highest_harmonic) /
                                np.log(2)))
    print(float_num_harmonics, num_harmonics, lut_size)
    lutsets.append(LUTentry(
        table=cos_lut(lut_size, 
                      harmonic_weights[:num_harmonics], 
                      harmonic_phases[:num_harmonics]),  # / lut_size,
        highest_harmonic=highest_harmonic))
    float_num_harmonics = bandwidth_factor * float_num_harmonics
  return lutsets

# Centered
#lutset = create_lutset(np.ones(128), np.arange(128) * -np.pi)
# Left-aligned
lutset = create_lutset(np.ones(128))
print([(len(s.table), s.highest_harmonic) for s in lutset])

plt.figure(figsize=(12, 6))
lut = lutset[0]
plotspec(lut.table, window_fn=np.ones)
print(np.sum(lutset[0].table))
print(np.sum(lutset[1].table))
print(np.sum(lutset[2].table))

def choose_from_lutset(period, lutset):
  """Select the best entry from a lutset for a given period.

  Args:
    period: (float) Target period of waveform in fractional samples.
    lutset: Sorted list of LUTs, as generated by create_lutset().
  
  Returns:
    One of the LUTs from the lutset, best suited to interpolating to generate
    a waveform of the desired period.
  """
  # Use the earliest (i.e., longest, most harmonics) LUT that works
  # (i.e., will not actually cause aliasing).
  # So start with the highest-bandwidth (and longest) LUTs, but skip them
  # if they result in aliasing.
  for lut_index, lutentry in enumerate(lutset):
    lut_table = lutentry.table
    lut_size = len(lut_table)
    # What proportion of nyquist does the highest harmonic in this table occupy?
    lut_bandwidth = 2 * lutentry.highest_harmonic / lut_size
    # To complete one cycle of <lut_size> points in <period> steps, each step
    # will need to be this many samples:
    lut_hop = lut_size / period
    # If we have a signal with a given bandwidth, but then speed it up by 
    # skipping lut_hop samples per sample, its bandwidth will increase 
    # proportionately.
    interp_bandwidth = lut_bandwidth * lut_hop
    if interp_bandwidth < 1.0:
      # No aliasing
      break
  # At this point, we either got to the end of the LUT table, or we found a
  # table we could interpolate without aliasing.
  print("lut_size=", lut_size, "lut_hop=", lut_hop, "lut_ix=", lut_index, "lut_bw=", lut_bandwidth)
  return lut_table

def interp_from_lutset(period, num_samples, lutset):
  """Generate a stretch of interpolated waveform from a set of LUTs.

  Args:
    period: (float) Target period of waveform in fractional samples.
    num_samples: (int) Total number of samples to generate at that period.
    lutset: Sorted list of LUTs, as generated by create_lutset().
  
  Returns:
    Vector of num_samples samples containing interpolated waveform with
    requested period.
  """
  lut_table = choose_from_lutset(period, lutset)
  lut_size = len(lut_table)
  # To complete one cycle of <lut_size> points in <period> steps, each step
  # will need to be this many samples:
  lut_hop = lut_size / period
  return interp(lut_table, np.arange(num_samples) * lut_hop / lut_size, order=3)

plotspec(interp_from_lutset(70, 210, lutset), window_fn=np.ones)

plotspec(interp_from_lutset(45, 135, lutset), window_fn=np.ones)

plotspec(interp_from_lutset(44, 132, lutset), window_fn=np.ones)

# Saving out the lutset as a C-compatible header file.

def write_lutset_to_h(filename, variable_base, lutset):
  num_luts = len(lutset)
  with open(filename, "w") as f:
    f.write("// Automatically-generated LUTset\n")
    f.write("#ifndef LUTSET_{:s}_DEFINED\n".format(variable_base.upper()))
    f.write("#define LUTSET_{:s}_DEFINED\n".format(variable_base.upper()))
    f.write("\n")
    # Define the structure.
    f.write("#ifndef LUTENTRY_DEFINED\n")
    f.write("#define LUTENTRY_DEFINED\n")
    f.write("typedef struct {\n")
    f.write("  const float *table;\n")
    f.write("  int table_size;\n")
    f.write("  int highest_harmonic;\n")
    f.write("} lut_entry;\n")
    f.write("#endif // LUTENTRY_DEFINED\n")
    f.write("\n")
    # Define the content of the individual tables.
    samples_per_row = 8
    for i in range(num_luts):
      table_size = len(lutset[i].table)
      f.write("const float {:s}_lutable_{:d}[{:d}] = {{\n".format(
          variable_base, i, table_size))
      for row_start in range(0, table_size, samples_per_row):
        for sample_index in range(row_start, 
                                  min(row_start + samples_per_row, table_size)):
          f.write("{:f},".format(lutset[i].table[sample_index]))
        f.write("\n")
      f.write("};\n")
      f.write("\n")
    # Define the table of LUTs.
    f.write("lut_entry {:s}_lutset[{:d}] = {{\n".format(
        variable_base, num_luts + 1))
    for i in range(num_luts):
      f.write("  {{{:s}_lutable_{:d}, {:d}, {:d}}},\n".format(
          variable_base, i, len(lutset[i].table), 
          lutset[i].highest_harmonic))
    # Final entry is null to indicate end of table.
    f.write("  {NULL, 0, 0},\n")
    f.write("};\n")
    f.write("\n")
    f.write("#endif // LUTSET_x_DEFINED\n")
  print("wrote", filename)

write_lutset_to_h('impulse_lutset.h', 'impulse', lutset)

!cat impulse_lutset.h

# Triangle wave lutset
n_harms = 64
coefs = (np.arange(n_harms) % 2)*(np.maximum(1, np.arange(n_harms, dtype=float))**(-2))
triangle_lutset = create_lutset(coefs, np.arange(len(coefs)) * -np.pi/2)

print([(len(s.table), s.highest_harmonic) for s in triangle_lutset])
plotspec(triangle_lutset[-3].table, window_fn=np.ones)

write_lutset_to_h('triangle_lutset.h', 'triangle', triangle_lutset)
!cat triangle_lutset.h

# Sinusoid "lutset" (only one table)
sine_lutset = create_lutset(np.array([0, 1]), length_factor=256)

print([(len(s.table), s.highest_harmonic) for s in sine_lutset])
plotspec(sine_lutset[0].table, window_fn=np.ones)

write_lutset_to_h('sine_lutset.h', 'sine', sine_lutset)
!cat sine_lutset.h

print(0.9*32768, 0.15 * 32768)

# Soft clipping lookup table scratchpad.

INT16_MAX = 32767
LINMAX = 29491  #// int(round(0.9 * 32768))
NONLIN_RANGE = 4915  # // size of nonlinearity lookup table = round(1.5 * (INT16_MAX - LINMAX))

clipping_lookup_table = np.arange(LINMAX + NONLIN_RANGE)
print(len(clipping_lookup_table))

for x in range(NONLIN_RANGE):
  x_dash = float(x) / NONLIN_RANGE
  clipping_lookup_table[x + LINMAX] = LINMAX + int(np.floor(NONLIN_RANGE * (x_dash - x_dash * x_dash * x_dash / 3.0)))

plt.plot(clipping_lookup_table[25000:]) # - np.arange(25000, len(clipping_lookup_table)))
print(max(clipping_lookup_table))

filename = "clipping_lookup_table.h"

with open(filename, "w") as f:
  f.write("// Automatically generated.\n// Clipping lookup table\n")
  f.write("const uint16 clipping_lookup_table[%d] = {\n" % NONLIN_RANGE)
  samples_per_row = 8
  for row_start in range(0, NONLIN_RANGE, samples_per_row):
    for sample in range(row_start, min(NONLIN_RANGE, row_start + samples_per_row)):
      f.write("%d," % clipping_lookup_table[LINMAX + sample])
    f.write("\n")
  f.write("};\n")
print("wrote", filename)

!cat clipping_lookup_table.h
# -*- coding: utf-8 -*- 2021-04-29-blosca_simulation.ipynb  Automatically generated by Colaboratory.  Original file is located at     https://colab.research.google.com/drive/1PSBC98l3rYV6GIQkevNhWRUU4SLY7CtA  # 2021-04-29-blosca_simulation.ipynb  Simulate the process of synthesizing waveforms by lookup of bandlimited impulse tables.  This is to prototype the C/C++ implementation in http://github.com/bwhitman/alles. # Linear interp # Centered #imp16 = cos_lut(256, np.ones(16), np.arange(16)*-np.pi) / 256 # Left-edge Synthesize <length> samples of a square wave of <period> by integrating the bandlimited impulse. # Phase is in cycles (1.0 = 1 cycle), i.e. radians/(2pi) # Actual step-size within LUT.  If <period> samples in output corresponds to len(LUT), # then each sample of output corresponds to a step of len(LUT)/period LUT samples. # Because stepping by <step> reduces the effective total number of samples in the LUT # by 1/step, and thus the total integral is reduced by 1/step, we scale the impulses # by step so their integral stays constant. # Smoothing time is 100x period for little droop. # To get the exact balance from the first sample, we have to discount half the first value of the # bandlimited impulse, because to start perfectly in the middle of the interpolate impulse, we # would start "half-way through" the zero'th sample. # The summation is a little off because it isn't exactly intg(imp), but has the decay forgetting. # Making the forgetting large makes it very close to the correct preset. #, initial_phase=0.25, initial_state=0.5) Synthesize <length> samples of a sawtooth of <period> by integrating the bandlimited impulse. # Phase is in cycles (1.0 = 1 cycle), i.e. radians/(2pi) # Actual step-size within LUT.  If <period> samples in output corresponds to len(LUT), # then each sample of output corresponds to a step of len(LUT)/period LUT samples. # Because stepping by <step> reduces the effective total number of samples in the LUT # by 1/step, and thus the total integral is reduced by 1/step, we scale the impulses # by step so their integral stays constant. # Smoothing time is 100x period for little droop. # To get the exact balance from the first sample, we have to discount half the first value of the # bandlimited impulse, because to start perfectly in the middle of the interpolate impulse, we # would start "half-way through" the zero'th sample. # The summation is a little off because it isn't exactly intg(imp), but has the decay forgetting. # Making the forgetting large makes it very close to the correct preset. #initial_state = 0 #return impulses #, initial_phase=0.25, initial_state=0.5) # The square wave is at +0.5 for period/2, so its integral would be period/4. # But we start it 1/4 through its cycle, so we can start the triangle at 0 # (midway through the upswing). # 500; + 0.012; 400; + 0.021; 300: + 0.034; 200: + 0.061; 100 + 0.143; 50: 0.307 #return sq_wave # We start at zero. #[:112]) # 256 * imp16 #, np.sum(sqint[:(3*period//4))]), np.sum(sqint[:period])) # Make a separate LUT for triangle. #print(coefs) # Implement the multiple lookup tables. # A LUT is stored as an array of values (table) and the harmonic number of the # highest harmonic they contain (i.e., the number of cycles it completes in the # entire table, so must be <= len(table)/2.) # A LUTset is a list of LUTentries describing downsampled versions of the same # basic waveform, sorted with the longest (highest-bandwidth) first. Create an ordered list of LUTs with decreasing harmonic content.    These can then be used in interp_from_lutset to make an adaptive-bandwidth   interpolation.    Args:     harmonic_weights: vector of amplitudes for cosine harmonic components.     harmonic_phases: initial phases for each harmonic, in radians.  Zero        (default) indicates cosine phase.     length_factor: Each table's length is at least this factor times the order       of the highest harmonic it contains.  Thus, this is a lower bound on the       number of samples per cycle for the highest harmonic.  Higher factors make       the interpolation easier.     bandwidth_factor: Target ratio between the highest harmonics in successive       table entries.  Default is sqrt(0.5), so after two tables, bandwidth is       reduced by 1/2 (and length with follow).    Returns:     A list of LUTentry objects, sorted in decreasing order of the highest      harmonic they contain.  Each LUT's length is a power of 2, and as small as     possible while respecting the length_factor for the highest contained      harmonic. # Calculate the length of the longest LUT we need.  Must be a power of 2, # must have at least length_factor * highest_harmonic samples. # Harmonic 0 (dc) doesn't count. # harmonic 0 is DC; there's no point in generating that table. # because zero doesn't count. # / lut_size, # Centered #lutset = create_lutset(np.ones(128), np.arange(128) * -np.pi) # Left-aligned Select the best entry from a lutset for a given period.    Args:     period: (float) Target period of waveform in fractional samples.     lutset: Sorted list of LUTs, as generated by create_lutset().      Returns:     One of the LUTs from the lutset, best suited to interpolating to generate     a waveform of the desired period. # Use the earliest (i.e., longest, most harmonics) LUT that works # (i.e., will not actually cause aliasing). # So start with the highest-bandwidth (and longest) LUTs, but skip them # if they result in aliasing. # What proportion of nyquist does the highest harmonic in this table occupy? # To complete one cycle of <lut_size> points in <period> steps, each step # will need to be this many samples: # If we have a signal with a given bandwidth, but then speed it up by # skipping lut_hop samples per sample, its bandwidth will increase # proportionately. # No aliasing # At this point, we either got to the end of the LUT table, or we found a # table we could interpolate without aliasing. Generate a stretch of interpolated waveform from a set of LUTs.    Args:     period: (float) Target period of waveform in fractional samples.     num_samples: (int) Total number of samples to generate at that period.     lutset: Sorted list of LUTs, as generated by create_lutset().      Returns:     Vector of num_samples samples containing interpolated waveform with     requested period. # To complete one cycle of <lut_size> points in <period> steps, each step # will need to be this many samples: # Saving out the lutset as a C-compatible header file. # Define the structure. # Define the content of the individual tables. # Define the table of LUTs. # Final entry is null to indicate end of table. # Triangle wave lutset # Sinusoid "lutset" (only one table) # Soft clipping lookup table scratchpad. #// int(round(0.9 * 32768)) # // size of nonlinearity lookup table = round(1.5 * (INT16_MAX - LINMAX)) # - np.arange(25000, len(clipping_lookup_table)))