diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 0000000..8a63c4a --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,19 @@ +{ + "configurations": [ + { + "name": "Mac", + "includePath": [ + "${workspaceFolder}/**" + ], + "defines": [], + "macFrameworkPath": [ + "/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/System/Library/Frameworks" + ], + "compilerPath": "/usr/bin/clang", + "cStandard": "c11", + "cppStandard": "c++98", + "intelliSenseMode": "macos-clang-x64" + } + ], + "version": 4 +} \ No newline at end of file diff --git a/playground_v2/exp_fab.py b/playground_v2/exp_fab.py index 8c1f97a..82a4689 100644 --- a/playground_v2/exp_fab.py +++ b/playground_v2/exp_fab.py @@ -9,6 +9,16 @@ import mz +class Test(mz.Module): + + def setup(self): + src = mz.SkewedTriangleSource(alpha=mz.Constant(0.99)) + f = mz.Butterworth(src, cutoffs=mz.LFO(lo=100, hi=500), + mode="lp") + self.out = f + + + class SignalWithEnvelope(mz.BaseModule): src: mz.BaseModule @@ -81,6 +91,14 @@ def setup(self): class Testing(mz.Module): def setup(self): + + src = mz.SkewedTriangleSource(frequency=mz.Parameter(220, knob="l_mixer_lo")) + f = mz.Butterworth(src, cutoffs=mz.LFO(lo=100, hi=500), + mode="lp") + self.out = f + return + + src = mz.SkewedTriangleSource( alpha=mz.Parameter(0.5, lo=0., hi=1., knob="r_mixer_mi"), frequency=mz.Parameter(300., key="f", @@ -154,45 +172,46 @@ def out_given_inputs(self, clock_signal: mz.ClockSignal, src): class Nick(mz.Module): - bpm: int = 130 - def setup(self): - bpm = mz.Constant(self.bpm) + bpm = mz.Constant(100) length = mz.Constant(5000.) - base_freq = 420. + base_freq = 320. freq = base_freq * PiecewiseLinearEnvelope( xs=[0., 0.2, 0.4, 0.7], - ys=[1., 0.7, 0.4, 0.1], - length=length) + ys=[0.9, 0.5, 0.3, 0.1], + length=length) #+ mz.LFO(mz.Constant(0.1), lo=100, hi=400) src = SignalWithEnvelope( src=mz.TimeIndependentSineSource(frequency=freq), env=mz.ADSREnvelopeGenerator(total_length=length*2)) trigger = mz.PeriodicTriggerSource(bpm) kick_trigger = trigger kick = mz.TriggerModulator(src, triggers=trigger) + de = mz.DelayElement(time=2, feedback=0.9, hi_cut=0.7) - kick = mz.SimpleDelay(kick, de, mix=0.4) - de = mz.DelayElement(time=5, feedback=0.6, hi_cut=0.7) + kick = mz.SimpleDelay(kick, de, mix=0.3) + de = mz.DelayElement(time=10, feedback=0.6, hi_cut=0.4) kick = mz.SimpleDelay(kick, de, mix=0.2) + # sNaRe - freq = 320 * PiecewiseLinearEnvelope( - xs=[0., 0.1, 0.3, 0.5, 0.9], - ys=[0.3, 0.4, 0.9, 0.4, 0.3], - length=length) + #freq = 320 * PiecewiseLinearEnvelope( + # xs=[0., 0.1, 0.3, 0.5, 0.9], + # ys=[0.1, 0.5, 1.0, 0.5, 0.1], + # #ys=[0.1, 0.1, 0.2, 0.1, 0.1], + # length=length) - src = mz.TimeIndependentSineSource(frequency=freq) + #src = mz.TimeIndependentSineSource(frequency=freq) src = Noise(mz.Constant(0.), eps=0.3) src = Overdrive(src, 0.8) src = Noise(src, eps=0.1) src = Overdrive(src, 0.8) + src = mz.ButterworthFilter(src, f_low=mz.Parameter(13000, key="f"), mode="lp") - src = SignalWithEnvelope( src=src, env=mz.ADSREnvelopeGenerator( @@ -206,13 +225,17 @@ def setup(self): snare = mz.TriggerModulator(src, triggers=trigger) + + #self.out = snare + kick * mz.Parameter(1., lo=0., hi=1., knob="r_mixer_hi") + #return + # hI hAt - freq = 320 * PiecewiseLinearEnvelope( - xs=[0., 0.1, 0.3, 0.5, 0.9], - ys=[0.3, 0.4, 0.9, 0.4, 0.3], - length=length) - src = mz.TimeIndependentSineSource(frequency=freq) +# freq = 320 * PiecewiseLinearEnvelope( +# xs=[0., 0.1, 0.3, 0.5, 0.9], +# ys=[0.1, 0.4, 1., 0.4, 0.1], +# length=length) + src = mz.TimeIndependentSineSource(frequency=mz.Constant(32)) src = NicStein(src, eps = 0.1) src = NicStein(src) @@ -231,29 +254,32 @@ def setup(self): total_length=mz.Constant(2000)) ) - trigger = mz.PeriodicTriggerSource(4 * bpm) + trigger = mz.PeriodicTriggerSource(bpm * 2) hi_hat = mz.TriggerModulator(src, triggers=trigger) + melody_cycler = 120. * mz.FreqFactors.STEP.value ** mz.Cycler( #(1, 1, 2, 2, 8, 8,) + - (1, 1, 0, 0, 0, 0, 1, 1, - 6, 6, 6, 8, 8, 9, 12, 12) + (1, 1, 1, 1, 0, 0, 0, 0, + 6, 6, 6, 6, 8, 8, 8, 8,) + #(1, 1, 0, 0, 0, 0, 1, 1, + #6, 6, 6, 8, 8, 9, 12, 12) ) melody = mz.TriggerModulator(melody_cycler, trigger) freq = mz.Hold(melody) + #freq = melody voice = mz.SkewedTriangleSource( alpha=mz.Parameter(0.9, lo=0., hi=1., knob="r_mixer_hi"), frequency=freq, ) melody_trigger = trigger - src = mz.ADSREnvelopeGenerator(attack=mz.Cycler((3, 8)),#mz.Constant(3), - hold=mz.Cycler((1., 4., 8.)), + src = mz.ADSREnvelopeGenerator(attack=mz.Constant(0.01),#mz.Cycler((3, 8)),#mz.Constant(3), + hold=mz.Constant(4.),#mz.Cycler((1., 4., 8.)), release=mz.Constant(5), - total_length=mz.LFO( - frequency=mz.Constant(0.01), - lo=1000,hi=20000), - ) + total_length=mz.Parameter(5000, 100, 10000, knob="r_mixer_lo")#mz.LFO( + #frequency=mz.Constant(0.01), + #lo=5000,hi=5000), #total_length=mz.Cycler(( # 4000, # 10000, @@ -265,7 +291,8 @@ def setup(self): # 6000, # 10000, # 4000, - # ))) + # )) + ) voice_env = mz.TriggerModulator(src, triggers=melody_trigger) voice_with_env = voice * voice_env voice_with_env = mz.SimpleDelay(voice_with_env, de, mix=0.2) @@ -273,9 +300,13 @@ def setup(self): voice = mz.ButterworthFilter( voice, - f_low=mz.Parameter(5000, lo=1., hi=24000, knob="r_mixer_lo", clip=True), mode="lp") + f_low=mz.Parameter(5000, lo=1., hi=24000, knob="r_mixer_mi", clip=True), mode="lp") + self.out = voice + kick + + return - self.out = hi_hat + kick + snare# + voice + + self.out = snare #+ kick + hi_hat# + voice class RandomDropper(mz.BaseModule): @@ -313,7 +344,7 @@ def setup(self): # kick = ) melody_cycler = 150. * mz.FreqFactors.STEP.value ** mz.Cycler( - (0, 0, 0, 3, 6, 8, 4, 3), + (0, 1, 12, 11, 10, 9), ) #melody_cycler = mz.Cycler((0.1, 0.5, 0.7)) trigger = mz.PeriodicTriggerSource(mz.Constant(bpm)) @@ -330,6 +361,7 @@ def setup(self): env = RandomDropper(env, p=0.1, drop_to = 0.2) melody = mz.TriggerModulator(env, trigger) melody = melody * sine + melody = Overdrive(melody * 2)/1.3 freq = NicStein(mz.Constant(0.), eps = 0.3)*80 + 250 @@ -350,22 +382,22 @@ def setup(self): weird_thing = mz.ButterworthFilter(werid_thing, f_low=mz.Parameter(1000, knob = "r_mixer_hi"), mode="lp") - self.out = weird_thing + melody + Nick(bpm) - return melody = mz.TriggerModulator(melody_cycler, trigger) voice = mz.SkewedTriangleSource( alpha=mz.Parameter(0.9, lo=0., hi=1., knob="r_mixer_hi"), frequency=freq, ) - melody_trigger = trigger_hi + melody_trigger = trigger src = mz.ADSREnvelopeGenerator(attack=mz.Cycler((3, 8)),#mz.Constant(3), hold=mz.Cycler((1., 4., 8., 20)), release=mz.Constant(5), total_length=mz.Cycler((20000, 26000, 20000))) voice_env = mz.TriggerModulator(src, triggers=melody_trigger) voice_with_env = voice * voice_env - voice_with_env = mz.SimpleDelay(voice_with_env, de, mix=0.2) + #voice_with_env = mz.SimpleDelay(voice_with_env, de, mix=0.2) voice = voice_with_env + self.out = weird_thing + melody + Nick() + return @@ -401,8 +433,9 @@ def setup(self): # f_high=mz.Parameter(5000, lo=1., hi=10000, key="f", clip=True), mode="hp") melody_cycler = 120. * mz.FreqFactors.STEP.value ** mz.Cycler( + (1, 2, 3, 8) #(12,)*8 + (8,)*8 + (5,)*8 + (0, )*4 + (1, 0, 1, 5)) - (12,)*8 + (8,)*8 + (5,)*8 + (0, )*4 + (5,)*8 + (8,)*8 + (12,)*8 + #(12,)*8 + (8,)*8 + (5,)*8 + (0, )*4 + (5,)*8 + (8,)*8 + (12,)*8 ) melody = mz.TriggerModulator(melody_cycler, trigger) freq = mz.Hold(melody) @@ -422,7 +455,7 @@ def setup(self): voice = mz.ButterworthFilter( voice, - f_low=mz.Parameter(5000, lo=1., hi=24000, knob="r_mixer_lo", clip=True), mode="lp") + f_low=mz.Parameter(5000, lo=1., hi=23000, knob="r_mixer_lo", clip=True), mode="lp") # Bassline @@ -480,7 +513,7 @@ def setup(self): kick = mz.SimpleDelay(kick, de, mix=0.2) melody_trigger = mz.PeriodicTriggerSource(bpm*2) - melody_cycler = 150. * mz.FreqFactors.STEP.value ** mz.Cycler((0, 0, 0, 1)) + melody_cycler = 150. * mz.FreqFactors.STEP.value ** RandomCycler((0,0,0,1)) #mz.Cycler((0, 0, 0, 1)) melody = mz.TriggerModulator(melody_cycler, trigger) freq = mz.Hold(melody) @@ -501,7 +534,7 @@ def setup(self): f_low=mz.Constant(6500), mode="lp") - mask = (1 + mz.SineSource(frequency=mz.Constant(0.1)))/2 + mask = (1 + mz.SineSource(frequency=mz.Constant(0.03)))/2 voice = voice_a * mask + voice_b * (1-mask) src = mz.ADSREnvelopeGenerator(attack=mz.Cycler((3, 8)),#mz.Constant(3), @@ -541,6 +574,145 @@ def setup(self): self.voice = mz.ButterworthFilter(self.voide) +class SimpleKick(mz.Module): + bpm: mz.Module + + def setup(self): + base_freq = 320. + length = mz.Constant(2000.) + freq = base_freq * PiecewiseLinearEnvelope( + xs=[0., 0.2, 0.4, 0.7], + ys=[1., 0.7, 0.4, 0.1], + length=length) + src = SignalWithEnvelope( + src=mz.TimeIndependentSineSource(frequency=freq), + env=mz.ADSREnvelopeGenerator(total_length=length*2)) + trigger = mz.PeriodicTriggerSource(self.bpm) + kick = mz.TriggerModulator(src, triggers=trigger) + self.out = kick + + +class Orca2(mz.Module): + # https://www.youtube.com/watch?v=gSFrBFBd7vY + def setup(self): + + bpm = mz.Constant(80) + triggers_3 = mz.PeriodicTriggerSource(bpm, note_value=1/3) + triggers_5 = mz.PeriodicTriggerSource(bpm, note_value=1/5) + triggers_7 = mz.PeriodicTriggerSource(bpm, note_value=1/7) + + base_freq = mz.Constant(160) + + kick = SimpleKick(bpm) + voice_3 = ClickVoice(triggers_3, base_freq) + voice_5 = ClickVoice(triggers_5, base_freq * mz.FreqFactors.STEP.value ** 5) + + voice_7 = ClickVoice(triggers_7, base_freq * mz.FreqFactors.STEP.value ** 7) + + self.out = kick + voice_3 + voice_5 + voice_7 + + +class ClickVoice(mz.BaseModule): + triggers: mz.BaseModule + frequency: mz.Module + + def setup(self): + voice = mz.SkewedTriangleSource( + alpha=mz.Constant(0.5), + frequency=self.frequency, + ) + click_env = mz.ADSREnvelopeGenerator(attack=mz.Constant(0.5), + hold=mz.Constant(1), + release=mz.Constant(3), + total_length=mz.Constant(3000)) # todo depends on bpm maybe + voice_env = mz.TriggerModulator(click_env, triggers=self.triggers) + self.out = voice * voice_env + +class AtmosphericVoice(mz.BaseModule): + pass + +class Orca(mz.Module): + + def setup(self): + bpm = mz.Constant(300) + C = 40 + + trigger_3 = mz.PeriodicTriggerSource(bpm * 105/3 / C) + trigger_5 = mz.PeriodicTriggerSource(bpm * 105/5 / C) + trigger_7 = mz.PeriodicTriggerSource(bpm * 105/7 / C) + + base_freq = 320. + length = mz.Constant(2000.) + freq = base_freq * PiecewiseLinearEnvelope( + xs=[0., 0.2, 0.4, 0.7], + ys=[1., 0.7, 0.4, 0.1], + length=length) + src = SignalWithEnvelope( + src=mz.TimeIndependentSineSource(frequency=freq), + env=mz.ADSREnvelopeGenerator(total_length=length*2)) + #bpm = mz.Constant(130) + trigger = mz.PeriodicTriggerSource(bpm) + kick = mz.TriggerModulator(src, triggers=trigger) +# de = mz.DelayElement(time=2, feedback=0.9, hi_cut=0.7) +# kick = mz.SimpleDelay(kick, de, mix=0.4) +# de = mz.DelayElement(time=5, feedback=0.6, hi_cut=0.7) +# kick = mz.SimpleDelay(kick, de, mix=0.2) + + trigger = mz.PeriodicTriggerSource(bpm) + melody_cycler = 80. * mz.FreqFactors.STEP.value ** mz.Cycler( + (1,)) + + mask = mz.Cycler((1, 0.01, 0.01)) + mask = mz.Hold(mz.TriggerModulator(mask, trigger)) + + melody = mz.TriggerModulator(melody_cycler, trigger) + freq = mz.Hold(melody) + voice = mz.SkewedTriangleSource( + alpha=mz.Parameter(0.9, lo=0., hi=1., knob="r_mixer_hi"), + frequency=freq, + ) * mask + + + melody_cycler = 120. * mz.FreqFactors.STEP.value ** mz.Cycler( + (1,)) + + mask = mz.Cycler((1, 0.01, 0.01, 0.01, 0.01)) + mask = mz.Hold(mz.TriggerModulator(mask, trigger)) + + melody = mz.TriggerModulator(melody_cycler, trigger) + freq = mz.Hold(melody) + voice2 = mz.SkewedTriangleSource( + alpha=mz.Parameter(0.9, lo=0., hi=1., knob="r_mixer_hi"), + frequency=freq, + ) * mask + + melody_cycler = 160. * mz.FreqFactors.STEP.value ** mz.Cycler( + (1,)) + + mask = mz.Cycler((1, 0.01, 0.01, 0.01, 0.01, 0.01, 0.1)) + mask = mz.Hold(mz.TriggerModulator(mask, trigger)) + + melody = mz.TriggerModulator(melody_cycler, trigger) + freq = mz.Hold(melody) + voice3 = mz.SkewedTriangleSource( + alpha=mz.Parameter(0.9, lo=0., hi=1., knob="r_mixer_hi"), + frequency=freq, + ) * mask + + + self.out = voice + kick + voice2 + voice3 + + return +# src = mz.ADSREnvelopeGenerator(attack=mz.Cycler((3, 8)),#mz.Constant(3), +# hold=mz.Cycler((1., 4., 8., 20)), +# release=mz.Constant(5), +# total_length=mz.Cycler((20000, 26000, 20000))) +# voice_env = mz.TriggerModulator(src, triggers=trigger) +# voice_with_env = voice * voice_env + + self.out = mz.SineSource() + + def _run(): t = Testing() diff --git a/playground_v2/mz/base.py b/playground_v2/mz/base.py index f3d469c..241e64a 100644 --- a/playground_v2/mz/base.py +++ b/playground_v2/mz/base.py @@ -1,4 +1,5 @@ """Abstrace base classes.""" + import dataclasses import warnings import enum @@ -46,7 +47,9 @@ def get(self): class ClockSignal(NamedTuple): - # Shape (num_samples,) + """Represents time steps and sample indices.""" + + # Time steps of samples (as real time). Shape (num_samples,) ts: np.ndarray # Shape (num_samples,) sample_indices: np.ndarray @@ -146,8 +149,9 @@ def get_clock_signal_num_samples(self, num_samples: int): return self.get_clock_signal(np.arange(num_samples, dtype=SAMPLE_INDICES_DTYPE)) -class NoCacheKeyError(Exception): - """Used to indicate that a module has no cache key.""" +# TODO: Revisit +# class NoCacheKeyError(Exception): +# """Used to indicate that a module has no cache key.""" # This is extracted from flax, and makes sure BaseModule subclasses get @@ -211,7 +215,7 @@ def setup(self): Notes: - Subclasses should use dataclass attributes to specify instance variables. - Like dataclasses, be careful with list efaults (use dataclasses.field). + Like dataclasses, be careful with list defaults (use dataclasses.field). - Subclasses can overwrite `setup` to provide custom initialization code Functinality: @@ -245,21 +249,20 @@ def __post_init__(self): self.monitor_sender = None - # Need to copy as it's updated inplace! + # Need to copy as it's updated inplace. keys_of_vars_before = set(vars(self)) self.setup() - # TODO: test setup vars_after = vars(self) names_of_variables_added_in_setup = { k for k, v in vars_after.items() if k not in keys_of_vars_before} potential_state_vars = names_of_variables_added_in_setup | non_module_fields - # Note that we directly unpack state variables here via `get()` + # Note that we directly unpack state variables here via `get()`. state_vars = {k: vars_after[k].get() for k in potential_state_vars if isinstance(vars_after[k], Stateful)} # Now set variables to the unpacked version, since we now know which ones they are. - # This allows users to then treat them as normal python types + # This allows users to then treat them as normal python types. for k, v in state_vars.items(): setattr(self, k, v) self._state_vars_names = tuple(state_vars.keys()) @@ -453,11 +456,21 @@ def out_single_value(self, clock_signal: ClockSignal) -> float: return np.mean(o) def out_given_inputs(self, clock_signal: ClockSignal, **inputs): + """Simple subclassing entry points. + + This can be used by subclasses that just need to + combine the numpy arrays that result from module calls. + `**inputs` will contain one element per `mz.Module`. + + TODO: Add example + """ raise NotImplementedError("Must be implemented by subclasses!") def __call__(self, clock_signal: ClockSignal): return self.out(clock_signal) + # Math support ------------------------------------------------------------ + def __mul__(self, other): """Implement module * scalar and module * module.""" return _MathModule(operator.mul, self, other) diff --git a/playground_v2/mz/filters.py b/playground_v2/mz/filters.py index b0a4394..4059ec4 100644 --- a/playground_v2/mz/filters.py +++ b/playground_v2/mz/filters.py @@ -9,6 +9,11 @@ import numpy as np from numpy.polynomial import Polynomial +from mz import pybind_backend + + +Butterworth = pybind_backend.Butterworth + @functools.lru_cache(maxsize=128) def basic_reverb_ir(delay: int, echo: int, p: float): diff --git a/playground_v2/mz/helpers.py b/playground_v2/mz/helpers.py index 8f12ff1..2f5d017 100644 --- a/playground_v2/mz/helpers.py +++ b/playground_v2/mz/helpers.py @@ -1,4 +1,6 @@ import collections +import time +import contextlib from typing import TypeVar import numpy as np @@ -102,3 +104,21 @@ def decorator(cls): def iter_marked_classes(): yield from (marked_class for marked_class in _MARKED_CLASSES) + + +class Timer: + + def __init__(self, emit_every:int = 100): + self._times = [] + self._emit_every = emit_every + + @contextlib.contextmanager + def __call__(self, name: str): + start = time.time() + yield + self._times.append(time.time() - start) + if len(self._times) % self._emit_every == 0: + avg = np.mean(self._times) + print(f">> {name}: {avg:.3e}s avg") + self._times = [] + diff --git a/playground_v2/mz/pybind_backend/__init__.py b/playground_v2/mz/pybind_backend/__init__.py new file mode 100644 index 0000000..ec46af5 --- /dev/null +++ b/playground_v2/mz/pybind_backend/__init__.py @@ -0,0 +1,97 @@ +# NEEDS + +# TODO: Only use the parts of torch that we need. +# pip install torch +# pip install ninja + +import functools +import os +import numpy as np +from torch.utils.cpp_extension import load + +from mz.experimental import subplots +from mz import base +from mz import helpers + + +# We load the CPP only on demand, as it involves a JIT ninja compilation step. +@functools.lru_cache() +def _cpp_backend(): + print("Loading CPP backend...") + pybind_backend_dir = os.path.dirname(os.path.realpath(__file__)) + src_dir = os.path.join(pybind_backend_dir, 'src') + return load( + name="cpp_backend", + sources=[os.path.join(src_dir, "backend.cpp")], + verbose=True) + + +# ------------------------------------------------------------------------------ +# Filters +# ------------------------------------------------------------------------------ + + +class Butterworth(base.Module): + """C++ backed Butterworth filter.""" + + src: base.Module + cutoffs: base.Module = base.Constant(200.) + + # This is only used for bp, in which case it's the high cutoff, + # and `cutoffs` is used as the low cutoff. + cutoffs_bp_hi: base.Module = base.Constant(1000.) + + # TODO: This is still WIP, we currently only support `lp` (low pass). + # TODO: Use enum. + mode: str = "lp" # One of lp, hp, bp. + order: int = 4 + + def setup(self): + + if self.mode == "bp": + raise NotImplementedError("Need High Pass for bandpass!") + elif self.mode == "hp": + raise NotImplementedError("High Pass not implemented yet!") + + # TODO: HP shoul use same number of coeffs, + # but backend not implemented yet. + elif self.mode in ("lp", "hp"): + if self.order % 2 != 0: + raise ValueError + n = self.order // 2 + self._w0 = np.zeros((n,), np.float32) + self._w1 = np.zeros((n,), np.float32) + self._w2 = np.zeros((n,), np.float32) + else: + raise ValueError(f"Invalid mode == {self.mode}!") + + # Will be created whenever num_samples changes. + self._result_buf = None + + # TODO: Only for book keeping, disable + self._timer = helpers.Timer() + + def out(self, clock_signal: base.ClockSignal): + src = self.src(clock_signal) + cutoffs = self.cutoffs(clock_signal) + # Validation. + if cutoffs.shape != src.shape: + raise ValueError("Cutoffs and source must have the same shape!") + if len(src.shape) != 1: + raise ValueError(f"Invalid src shape = {src.shape}!") + # Create output buffer if needed. + if self._result_buf is None or self._result_buf.shape != src.shape: + print("Creating result_buf...") + self._result_buf = np.empty(src.shape, np.float32) + # Run filter. + if self.mode == "bp": + # Fetch highs. + cutoffs_hi = self.cutoffs_bp_hi(clock_signal) + raise NotImplementedError # Swee above. + else: + with self._timer("Time in C++"): + _cpp_backend().butterworth_lowpass_or_highpass( + src, cutoffs, self._result_buf, + self._w0, self._w1, self._w2, + clock_signal.sample_rate, self.mode == "lp") + return self._result_buf diff --git a/playground_v2/mz/pybind_backend/pybind_backend.py b/playground_v2/mz/pybind_backend/pybind_backend.py deleted file mode 100644 index bf09215..0000000 --- a/playground_v2/mz/pybind_backend/pybind_backend.py +++ /dev/null @@ -1,20 +0,0 @@ -# NEEDS -# pip install torch -# pip install ninja - -import os -from torch.utils.cpp_extension import load - -pybind_backend_dir = os.path.dirname(os.path.realpath(__file__)) -src_dir = os.path.join(pybind_backend_dir, 'src') -src = load( - name="src", - sources=[os.path.join(src_dir, "test.cpp")], - verbose=True) - - -import numpy as np -a = np.ones((128,), np.float32) - -b = src.add(a,a) -print(b) diff --git a/playground_v2/mz/pybind_backend/src/backend.cpp b/playground_v2/mz/pybind_backend/src/backend.cpp new file mode 100644 index 0000000..ed88cc9 --- /dev/null +++ b/playground_v2/mz/pybind_backend/src/backend.cpp @@ -0,0 +1,170 @@ + +#include +#include +#include + +#include +#include + + +constexpr float pi = 3.14159265358979323846; + + +/******************************************************************************/ +/* Filters */ +/******************************************************************************/ + + +void butterworth_lowpass_or_highpass( + // Input signal, shape (num_samples,) + py::array_t signal_arr, + // Cutoff freqs, shape (num_samples,) + py::array_t cutoffs_arr, + // Output array, will write here. + py::array_t result_arr, + // Window coefficients, shape (n,) + py::array_t w0_arr, + py::array_t w1_arr, + py::array_t w2_arr, + float sampling_rate, + // This toggles low pass or high pass behavior. + bool low_pass +) { + /** Butterworth Lowpass OR Highpass filter. + + The cutoff frequency may be different for each sample. + */ + + py::buffer_info signal_buf = signal_arr.request(); + py::buffer_info cutoffs_buf = cutoffs_arr.request(); + py::buffer_info result_buf = result_arr.request(); + + int num_samples = signal_buf.shape[0]; + + // Get pointers of signals. + float *signal = (float *) signal_buf.ptr, + *result = (float *) result_buf.ptr, + *cutoffs = (float *) cutoffs_buf.ptr; + + // Get pointers for w0, w1, w2. + py::buffer_info w0_buf = w0_arr.request(), + w1_buf = w1_arr.request(), + w2_buf = w2_arr.request(); + float *w0 = (float *) w0_buf.ptr, + *w1 = (float *) w1_buf.ptr, + *w2 = (float *) w2_buf.ptr; + + // This is ==order/2. + int n = w0_buf.shape[0]; + int actual_n = n*2; + + // Current output. + float x = 0; + float w1sign = low_pass ? 1 : -1; + + float x0 = 0, // Holds x[i] + x1 = 0, // Holds x[i-1] + x2 = 0; // Holds x[i-2] + + float xk0 = 0, // Holds x_k[i] + xk1 = 0, // Holds x_k[i-1] + xk2 = 0; // Holds x_k[i-w] + + float *y0 = w0, + *y1 = w1, + *y2 = w2; + + for (size_t sample_i = 0; sample_i < num_samples; sample_i++) { + x = signal[sample_i]; + const float cutoff_i = cutoffs[sample_i]; + + /* Discretize butterworth via second-order sections. + + Based on [1]. + + We limit ourselfs to n even. We can factor the transfer function H_n,lp(s) of a + n-th order butterworth lowpass into a product of second order sections (see [1]): + + H_n,lp(s) = \prod_{k=0}^{n/2-1} H_2k(s), + + As shown in [1], H_2k(s) becomes the following in the z-domain: + + 1 + 2 z^-1 + z^-2 Y(z) + H_2k(z) = ---------------------------- = ---- + b0_k + b1_k z^-1 + b2_k z^-2 X(z) + + where b0_k = gamma^2 - alpha_k gamma + 1, + b1_k = 2 - 2\gamma^2 + b2_k = gamma^2 + alpha_k gamma + 1, + + gamma, alpha_k as in [1]. + + Converting H_2k(s) into discrete space, we land at: + + b0_k y_k[i] + b1_k y_k[i-1] + b2_k y_k[i-2] + = x_k[i] + 2 x_k[i-1] + x_k[i-2] (1.) + b0_k y_k[i] = - b1_k y_k[i-1] - b2_k y_k[i-2] + x_k[i] + 2 x_k[i-1] + x_k[i-2] + y_k[i] = 1/b0_k * (- b1_k y_k[i-1] - b2_k y_k[i-2] + x_k[i] + 2 x_k[i-1] + x_k[i-2]) + + Here, `i` indexes the time dimension (i.e., it goes from t to t+num_samples). + Note the sub `k` at both y and x, this is due to using SOS. Our input is given as + x_in[n] and our real output we want is y_out[n]. We can calulate y_out by iterating Eq 1. + At each step, we set x_k[i] = y_{k-1}[i], hwere x_0[i] = x_in[i]. + The final y_out becomes y_n. + + For the _high pass_, we have that + + H_n,hp(s) = s^n * H_n,lp(s) + + and thus + + H_2k,hp(s) = H_2k(s) * s^2 + + TODO: Derive code! + + References + ---------- + + [1]: https://tttapa.github.io/Pages/Mathematics/Systems-and-Control-Theory/Digital-filters/Discretization/Discretization-of-a-fourth-order-Butterworth-filter.html + + */ + + x0 = signal[sample_i]; + const float gammaF = 1/(std::tan(pi*cutoff_i/sampling_rate)); + // We set xk[t] = x[t] for i = {t, t-1, t-2}. + xk0 = x0; + xk1 = x1; + xk2 = x2; + for (size_t k=0; k <= (actual_n/2-1); k++){ + const float alpha_k = 2 * std::cos(2*pi * (2.*k + actual_n + 1)/(4.*actual_n)); + const float b0k = gammaF * gammaF - alpha_k * gammaF + 1.; + const float b1k = 2. - 2. * gammaF*gammaF; + const float b2k = gammaF * gammaF + alpha_k * gammaF + 1.; + y0[k] = (-b1k*y1[k] - b2k*y2[k] + xk0 + 2. * xk1 + xk2) / b0k; + + xk0 = y0[k]; + xk1 = y1[k]; + xk2 = y2[k]; + + // Shift ys + y2[k] = y1[k]; + y1[k] = y0[k]; + } + + x2 = x1; + x1 = x0; + + result[sample_i] = y0[actual_n/2-1]; + } +} + +/******************************************************************************/ +/* Exporting to Python */ +/******************************************************************************/ + + +PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { + m.doc() = "Backend"; + m.def("butterworth_lowpass_or_highpass", &butterworth_lowpass_or_highpass, + "Butterworth Low Pass OR High Pass"); +} diff --git a/playground_v2/mz/pybind_backend/src/test.cpp b/playground_v2/mz/pybind_backend/src/test.cpp deleted file mode 100644 index 77f5a3d..0000000 --- a/playground_v2/mz/pybind_backend/src/test.cpp +++ /dev/null @@ -1,38 +0,0 @@ - -#include -#include -#include - - -py::array_t add(py::array_t input1, py::array_t input2) { - py::buffer_info buf1 = input1.request(); - py::buffer_info buf2 = input2.request(); - - if (buf1.size != buf2.size) { - throw std::runtime_error("Input shapes must match"); - } - - /* allocate the buffer */ - py::array_t result = py::array_t(buf1.size); - - py::buffer_info buf3 = result.request(); - - float *ptr1 = (float *) buf1.ptr, - *ptr2 = (float *) buf2.ptr, - *ptr3 = (float *) buf3.ptr; - int X = buf1.shape[0]; - - printf("HI %d", X); - - for (size_t idx = 0; idx < X; idx++) { - ptr3[idx] = ptr1[idx] + ptr2[idx]; - } - - return result; -} - - -PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) { - m.doc() = "pybind11 example plugin"; // optional module docstring - m.def("add", &add, "A function that adds two numbers"); -} \ No newline at end of file