ii-sound/sim_pwm.py

141 lines
4.1 KiB
Python

import argparse
import math
import numpy
import librosa
import soundfile as sf
def params(freq, damping, dt):
w = freq * 2 * math.pi * dt
d = damping * dt
e = math.exp(d)
c1 = 2 * e * math.cos(w)
c2 = e * e
t0 = (1 - 2 * e * math.cos(w) + e * e) / (d * d + w * w)
t = d * d + w * w - math.pi * math.pi
t1 = (1 + 2 * e * math.cos(w) + e * e) / math.sqrt(t * t + 4 * d * d *
math.pi * math.pi)
b2 = (t1 - t0) / (t1 + t0)
b1 = b2 * dt * dt * (t0 + t1) / 2
return c1, c2, b1, b2
def filter_audio(
audio: numpy.ndarray, sample_rate: float,
sim_rate: float) -> numpy.ndarray:
freq = 3875
dt = 1 / sim_rate
damping = -1210
c1, c2, b1, b2 = params(freq, damping, dt)
y1 = y2 = 0
x1 = 0
x2 = 0
x1 = 1.0
scale = 1150 # TODO: analytic expression
maxy = 0
cycles_per_sample = int(sim_rate / sample_rate)
audio_idx = 0
sample = audio[0]
duty_cycle = int((sample + 1) * cycles_per_sample / 2)
next_toggle = duty_cycle
cycles = 0
y = 0
bias = 1.0
while audio_idx < len(audio) / 1:
if audio_idx % 10000 == 0:
print(audio_idx, len(audio), y / scale, x1, sample, duty_cycle,
next_toggle)
cycles += 1
# x1 += bias / 4
# if x1 > 1.0:
# x1 = 1.0
# if x1 < -1.0:
# x1 = -1.0
y = (c1 * y1 - c2 * y2 + b1 * x1 + b2 * x2)
x2 = x1
if cycles >= next_toggle:
bias *= -1.0
x1 *= -1.0
# x1 += bias / 4
# if x1 > 1.0:
# x1 = 1.0
# if x1 < -1.0:
# x1 = -1.0
if bias == 1.0:
audio_idx += 1
if audio_idx < len(audio):
sample = audio[audio_idx]
duty_cycle = int((sample + 1) * cycles_per_sample / 2)
else:
duty_cycle = cycles_per_sample - duty_cycle
next_toggle += duty_cycle
y2 = y1
y1 = y
if math.fabs(y) > maxy:
maxy = math.fabs(y)
yield y / scale
print("scale = %f" % maxy)
def preprocess(
filename: str, target_sample_rate: int, normalize: float,
normalization_percentile: int) -> numpy.ndarray:
"""Upscale input audio to target sample rate and normalize signal."""
data, _ = librosa.load(filename, sr=target_sample_rate, mono=True)
max_value = numpy.percentile(data, normalization_percentile)
data /= max_value
data *= normalize
return data
def main():
parser = argparse.ArgumentParser()
parser.add_argument("--clock", choices=['pal', 'ntsc'],
help="Whether target machine clock speed is PAL ("
"1015657Hz) or NTSC (1020484)",
required=True)
parser.add_argument("--normalization", default=1.0, type=float,
help="Overall multiplier to rescale input audio "
"values.")
parser.add_argument("--norm_percentile", default=100,
help="Normalize to specified percentile value of input "
"audio")
parser.add_argument("input", type=str, help="input audio file to convert")
parser.add_argument("output", type=str, help="output audio file")
args = parser.parse_args()
# Effective clock rate, including every-65 cycle "long cycle" that takes
# 16/14 as long.
sim_rate = 1015657 if args.clock == 'pal' else 1020484 # NTSC
input_audio = preprocess(args.input, sim_rate / 46, args.normalization,
args.norm_percentile)
output_audio = numpy.array(
list(filter_audio(input_audio, sim_rate / 46, sim_rate)),
dtype=numpy.float32)
output_rate = 96000
output = librosa.resample(output_audio, orig_sr=sim_rate,
target_sr=output_rate)
with sf.SoundFile(
args.output, "w", output_rate, channels=1, format='WAV') \
as f:
f.write(output)
if __name__ == "__main__":
main()