fir filtering demo

tests/circuitpython-manual/synthio/note/fir.py

@@ -9,6 +9,7 @@ import audiocore
 import synthio
 from ulab import numpy as np
 import adafruit_wave as wave
+import mkfilter
 
 random.seed(9)
 
@@ -16,89 +17,37 @@ envelope = synthio.Envelope(
     attack_time=0.1, decay_time=0.05, release_time=0.2, attack_level=0.8, sustain_level=0.8
 )
 
-h = np.array(
-    [
-        -0.001229734800309099,
-        -0.008235561806605458,
-        -0.015082497016061390,
-        -0.020940136918319988,
-        -0.024981800822463429,
-        -0.026464233332370746,
-        -0.024803890156806906,
-        -0.019642276775473012,
-        -0.010893620860173042,
-        0.001230341899766145,
-        0.016221637398855598,
-        0.033304135659230648,
-        0.051486665261155681,
-        0.069636961761409016,
-        0.086570197432542767,
-        0.101144354207918147,
-        0.112353938422488253,
-        0.119413577288191297,
-        0.121823886314051028,
-        0.119413577288191297,
-        0.112353938422488253,
-        0.101144354207918147,
-        0.086570197432542767,
-        0.069636961761409016,
-        0.051486665261155681,
-        0.033304135659230648,
-        0.016221637398855598,
-        0.001230341899766145,
-        -0.010893620860173042,
-        -0.019642276775473012,
-        -0.024803890156806906,
-        -0.026464233332370746,
-        -0.024981800822463429,
-        -0.020940136918319988,
-        -0.015082497016061390,
-        -0.008235561806605458,
-        -0.001229734800309099,
-    ]
-)
+SAMPLE_SIZE = 1024
+bend_out = np.linspace(0, 32767, num=SAMPLE_SIZE, endpoint=True, dtype=np.int16)
 
-filter_coeffs = np.array(h[::-1] * 32768, dtype=np.int16)
-
-synth = synthio.Synthesizer(sample_rate=48000, filter=filter_coeffs)
+filter_rectangular = mkfilter.LPF(48000, 800, 13)
+filter_rectangular_big = mkfilter.LPF(48000, 800, 59)
+filter_blackman = mkfilter.LPF(48000, 800, 59, win=mkfilter.blackman)
+print(filter_blackman)
 
 
 def synthesize(synth):
     n = synthio.Note(
         frequency=120,
         envelope=envelope,
-        filter=False,
+        filter=True,
+        bend=synthio.LFO(bend_out, once=True, rate=1 / 2, scale=5),
     )
 
+    synth.press(n)
     print(synth, n)
-    synth.press((n,))
-    for _ in range(20):
-        n.frequency *= 1.0595
-        yield 36
-    synth.release_all()
-    yield 36
-
-    n.filter = True
-    n.frequency = 120
-    synth.press((n,))
-    for _ in range(20):
-        n.frequency *= 1.0595
-        yield 36
+    yield 2 * 48000 // 256
     synth.release_all()
     yield 36
 
 
-def chain(*args):
-    for a in args:
-        yield from a
-
-
-# sox -r 48000 -e signed -b 16 -c 1 tune.raw tune.wav
 with wave.open("fir.wav", "w") as f:
     f.setnchannels(1)
     f.setsampwidth(2)
     f.setframerate(48000)
-    for n in chain(synthesize(synth)):
-        for i in range(n):
-            result, data = audiocore.get_buffer(synth)
-            f.writeframes(data)
+    for filter_coeffs in [None, filter_rectangular, filter_rectangular_big, filter_blackman]:
+        synth = synthio.Synthesizer(sample_rate=48000, filter=filter_coeffs)
+        for n in synthesize(synth):
+            for i in range(n):
+                result, data = audiocore.get_buffer(synth)
+                f.writeframes(data)

tests/circuitpython-manual/synthio/note/mkfilter.py

@@ -0,0 +1,105 @@
+try:
+    from ulab import numpy as np
+except ImportError:
+    import numpy as np
+
+
+def lpf(fS, f, N, win=lambda N: 1):
+    if not (N & 1):
+        raise ValueError("filter length must be odd")
+    h = np.sinc(2 * f / fS * (np.arange(N) - (N - 1) / 2))
+    h = h * win(N)
+    return h * (1 / np.sum(h))
+
+
+def hpf(fS, f, N, win=lambda N: 1):
+    if not (N & 1):
+        raise ValueError("filter length must be odd")
+    h = -lpf(fS, f, N)
+    h = h * win(N)
+    h[(N - 1) // 2] += 1
+    return h
+
+
+def brf(fS, fL, NL, fH, NH, win=lambda N: 1):
+    hlpf = lpf(fS, fL, NL, win)
+    hhpf = hpf(fS, fH, NH, win)
+
+    if NH > NL:
+        h = hhpf
+        h[(NH - NL) // 2 : (NH - NL) // 2 + NL] += hlpf
+    else:
+        h = hlpf
+        h[(NL - NH) // 2 : (NL - NH) // 2 + NH] += hhpf
+
+    return h
+
+
+def bpf(fS, fL, NL, fH, NH, win=lambda N: 1):
+    hlpf = lpf(fS, fL, NL, win)
+    hhpf = hpf(fS, fH, NH, win)
+    return np.convolve(hlpf, hhpf)
+
+
+def blackman(M):
+    n = np.arange(1 - M, M, 2)
+    return 0.42 + 0.5 * np.cos(np.pi * n / (M - 1)) + 0.08 * np.cos(2.0 * np.pi * n / (M - 1))
+
+
+def tosynthio(coeffs):
+    result = np.array(coeffs * 32767, dtype=np.int16)
+    return trim_zeros(result)
+
+
+def trim_zeros(arr):
+    i = 0
+    j = len(arr) - 1
+    while i < len(arr) and arr[i] == 0:
+        i += 1
+    while j > i and arr[j] == 0:
+        j -= 1
+    return arr[i : j + 1]
+
+
+# fiiir.com uses factor 4.6 for blackman window, 0.91 for rectangular
+def ntaps(fS, fB, factor=4.6):
+    b = fB / fS
+    return round(factor / b) | 1
+
+
+def LPF(*args, **kw):
+    return tosynthio(lpf(*args, **kw))
+
+
+def HPF(*args, **kw):
+    return tosynthio(hpf(*args, **kw))
+
+
+def BRF(*args, **kw):
+    return tosynthio(brf(*args, **kw))
+
+
+def BPF(*args, **kw):
+    return tosynthio(bpf(*args, **kw))
+
+
+if __name__ == "__main__":
+    print("lpf(24000, 2040, 13) # 1920Hz transition window")
+    print(list(lpf(24000, 2040, 13)))
+
+    print("hpf(24000, 9600, 13) # 960Hz transition window")
+    print(list(hpf(24000, 9600, 23)))
+
+    print("bpf(24000, 1200, 11, 3960, 15) #  2400Hz, 1600Hz transition windows")
+    print(list(bpf(24000, 1200, 11, 3960, 15)))
+
+    print("brf(24000, 960, 19, 2400, 13) # 1200, 1800Hz transition windows")
+    brf_tst = brf(24000, 960, 19, 2400, 13)
+    print(brf_tst)
+
+    print("brf(24000, 960, 13, 2400, 19) # 1200, 1800Hz transition windows")
+    brf_tst = brf(24000, 960, 13, 2400, 19)
+    print(brf_tst)
+
+    print("lpf(1, 0.1, 59, blackman) # 1920Hz transition window, blackman")
+    print(lpf(1, 0.1, 59, blackman))