Add more chroma features #2

examples/chromatest.py

@@ -43,7 +43,7 @@ def main():
 
     # 2) analyze previously built test signal
 
-    chroma = Chroma(sr, cp)
+    chroma = Chroma(sr, cp, feature='cent')
 
     chromagram = chroma.chroma(x)[-1]  # pick last chromagram entry
 

src/python/qdft/chroma.py

@@ -9,13 +9,14 @@
 
 import numpy
 
+from .fafe import QFAFE
 from .qdft import QDFT
 from .scale import Scale
 
 
 class Chroma:
 
-    def __init__(self, samplerate, concertpitch=440, bandwidth=('A0', 'C#8'), decibel=True):
+    def __init__(self, samplerate, concertpitch=440, bandwidth=('A0', 'C#8'), decibel=True, feature=None):
 
         scale = Scale(concertpitch)
 
@@ -39,6 +40,7 @@ def __init__(self, samplerate, concertpitch=440, bandwidth=('A0', 'C#8'), decibe
         self.concertpitch = concertpitch
         self.bandwidth = bandwidth
         self.decibel = decibel
+        self.feature = feature
         self.semitones = semitones
         self.frequencies = frequencies
         self.notes = notes
@@ -49,77 +51,35 @@ def __init__(self, samplerate, concertpitch=440, bandwidth=('A0', 'C#8'), decibe
 
     def chroma(self, samples):
 
-        stash = { 'cents': None }
-
-        def analysis(dfts, mode=None):
-            stash['cents'] = self.analyze(dfts, mode)
-
-        # TODO: analyze raw dfts
-        # dfts = self.qdft.qdft(samples, analysis)
-
-        # TODO: analyze windowed dfts
         dfts = self.qdft.qdft(samples)
-        stash['cents'] = self.analyze(dfts, 'p')
 
-        magnis = numpy.abs(dfts)
-        cents = stash['cents']
+        magnitudes = numpy.abs(dfts)
+        features = None
 
         if self.decibel:
 
             with numpy.errstate(all='ignore'):
-                magnis = 20 * numpy.log10(magnis)
-
-        chromagram = magnis + 1j * cents
-
-        chromagram = chromagram[..., ::2]
-        assert chromagram.shape[-1] == self.frequencies.shape[-1]
-
-        return chromagram
+                magnitudes = 20 * numpy.log10(magnitudes)
 
-    def analyze(self, dfts, mode=None):
+        if str(self.feature).lower() in 'phase':
 
-        l = numpy.roll(dfts, +1, axis=-1)
-        m = dfts
-        r = numpy.roll(dfts, -1, axis=-1)
+            features = numpy.angle(dfts)
 
-        if mode is None:
+        if str(self.feature).lower() in 'hz':
 
-            with numpy.errstate(all='ignore'):
-                drifts = -numpy.real((r - l) / (2 * m - r - l))
-
-        elif str(mode).lower() == 'p':
+            fafe = QFAFE(self.qdft)
+            features = fafe.hz(dfts)
 
-            p = 1.36
+        if str(self.feature).lower() in 'cent':
 
-            l = numpy.abs(l)
-            m = numpy.abs(m)
-            r = numpy.abs(r)
+            fafe = QFAFE(self.qdft)
+            features = fafe.cent(dfts)
 
-            with numpy.errstate(all='ignore'):
-                drifts = p * (r - l) / (m + r + l)
+        chromagram = (magnitudes + 1j * features) \
+                     if features is not None \
+                     else magnitudes
 
-        elif str(mode).lower() == 'q':
-
-            q = 0.55
-
-            with numpy.errstate(all='ignore'):
-                drifts = -numpy.real(q * (r - l) / (2 * m + r + l))
-
-        else:
-
-            drifts = numpy.zeros(dfts.shape)
-
-        drifts[...,  0] = 0
-        drifts[..., -1] = 0
-
-        oldfreqs = self.qdft.frequencies
-        oldbins = numpy.arange(oldfreqs.size)
-        newbins = oldbins + drifts
-        # TODO: is approximation possible? https://en.wikipedia.org/wiki/Cent_(music)
-        newfreqs = self.qdft.bandwidth[0] * numpy.power(2, newbins / self.qdft.resolution)
-        # TODO: does interp make sense?
-        # newfreqs = numpy.interp(newbins, oldbins, oldfreqs)
-
-        cents = 1200 * numpy.log2(newfreqs / oldfreqs)
+        chromagram = chromagram[..., ::2]
+        assert chromagram.shape[-1] == self.frequencies.shape[-1]
 
-        return cents
+        return chromagram

src/python/qdft/fafe.py

@@ -0,0 +1,95 @@
+"""
+Copyright (c) 2023 Juergen Hock
+
+SPDX-License-Identifier: MIT
+
+Fast, Accurate Frequency Estimator according to [1].
+
+[1] Eric Jacobsen and Peter Kootsookos
+    Fast, Accurate Frequency Estimators
+    IEEE Signal Processing Magazine (2007)
+    https://ieeexplore.ieee.org/document/4205098
+
+Source: https://github.com/jurihock/qdft
+"""
+
+
+import numpy
+
+
+class FAFE:
+
+    def __init__(self, mode=None):
+
+        self.mode = mode
+
+    def __call__(self, dfts):
+
+        dfts = numpy.atleast_2d(dfts)
+        assert dfts.ndim == 2
+
+        l = numpy.roll(dfts, +1, axis=-1)
+        m = dfts
+        r = numpy.roll(dfts, -1, axis=-1)
+
+        if self.mode is None:
+
+            with numpy.errstate(all='ignore'):
+                drifts = -numpy.real((r - l) / (2 * m - r - l))
+
+        elif str(self.mode).lower() == 'p':
+
+            p = 1.36  # TODO: hann?
+
+            l = numpy.abs(l)
+            m = numpy.abs(m)
+            r = numpy.abs(r)
+
+            with numpy.errstate(all='ignore'):
+                drifts = p * (r - l) / (m + r + l)
+
+        elif str(self.mode).lower() == 'q':
+
+            q = 0.55  # TODO: hann?
+
+            with numpy.errstate(all='ignore'):
+                drifts = -numpy.real(q * (r - l) / (2 * m + r + l))
+
+        else:
+
+            drifts = numpy.zeros(dfts.shape)
+
+        drifts[...,  0] = 0
+        drifts[..., -1] = 0
+
+        return drifts
+
+
+class QFAFE:
+
+    def __init__(self, qdft):
+
+        self.qdft = qdft
+        self.fafe = FAFE('p' if qdft.window is not None else None)
+
+    def hz(self, dfts):
+
+        oldfreqs = self.qdft.frequencies
+
+        oldbins = numpy.arange(oldfreqs.size)
+        newbins = oldbins + self.fafe(dfts)
+
+        # TODO: is approximation possible? https://en.wikipedia.org/wiki/Cent_(music)
+        newfreqs = self.qdft.bandwidth[0] * numpy.power(2, newbins / self.qdft.resolution)
+
+        # TODO: does interp make sense?
+        # newfreqs = numpy.interp(newbins, oldbins, oldfreqs)
+
+        return newfreqs
+
+    def cent(self, dfts):
+
+        newfreqs = self.hz(dfts)
+        oldfreqs = self.qdft.frequencies
+
+        return 1200 * numpy.log2(newfreqs / oldfreqs)

src/python/qdft/qdft.py

@@ -67,7 +67,7 @@ def __init__(self, samplerate, bandwidth, resolution=24, latency=0, window=(+0.5
         self.outputs = outputs
         self.kernels = kernels
 
-    def qdft(self, samples, callback=None):
+    def qdft(self, samples):
 
         samples = numpy.atleast_1d(samples).astype(float)
         assert samples.ndim == 1
@@ -103,14 +103,6 @@ def qdft(self, samples, callback=None):
 
             outputs[k] = dfts[k][-1]
 
-        if callback is not None:
-
-            assert callable(callback)
-
-            dfts[0].setflags(write=False)
-            try: callback(dfts[0])
-            finally: dfts[0].setflags(write=True)
-
         if window is not None:
 
             a, b = window[0], window[1] / 2