src/pitch/pitch.h: add a note about yinfast
authorPaul Brossier <piem@piem.org>
Tue, 9 May 2017 14:06:23 +0000 (16:06 +0200)
committerPaul Brossier <piem@piem.org>
Tue, 9 May 2017 14:06:23 +0000 (16:06 +0200)
python/demos/demo_fastyin_compare.py [deleted file]
python/demos/demo_yin_compare.py [new file with mode: 0755]
src/pitch/pitch.h

diff --git a/python/demos/demo_fastyin_compare.py b/python/demos/demo_fastyin_compare.py
deleted file mode 100755 (executable)
index 92dc015..0000000
+++ /dev/null
@@ -1,175 +0,0 @@
-#! /usr/bin/env python
-# -*- coding: utf8 -*-
-
-""" Pure python implementation of the sum of squared difference
-
-    sqd_yin: original sum of squared difference [0]
-        d_t(tau) = x ⊗ kernel
-    sqd_yinfast: sum of squared diff using complex domain [0]
-    sqd_yinfftslow: tappered squared diff [1]
-    sqd_yinfft: modified squared diff using complex domain [1]
-
-[0]:http://audition.ens.fr/adc/pdf/2002_JASA_YIN.pdf
-[1]:https://aubio.org/phd/
-"""
-
-import sys
-import numpy as np
-import matplotlib.pyplot as plt
-
-def sqd_yin(samples):
-    """ compute original sum of squared difference
-
-    Brute-force computation (cost o(N**2), slow)."""
-    B = len(samples)
-    W = B//2
-    yin = np.zeros(W)
-    for j in range(W):
-      for tau in range(1, W):
-        yin[tau] += (samples[j] - samples[j+tau])**2
-    return yin
-
-def sqd_yinfast(samples):
-    """ compute approximate sum of squared difference
-
-    Using complex convolution (fast, cost o(n*log(n)) )"""
-    # yin_t(tau) = (r_t(0) + r_(t+tau)(0)) - 2r_t(tau)
-    B = len(samples)
-    W = B//2
-    yin = np.zeros(W)
-    sqdiff = np.zeros(W)
-    kernel = np.zeros(B)
-    # compute r_(t+tau)(0)
-    squares = samples**2
-    for tau in range(W):
-        sqdiff[tau] = squares[tau:tau+W].sum()
-    # add r_t(0)
-    sqdiff += sqdiff[0]
-    # compute r_t(tau) using kernel convolution in complex domain
-    samples_fft = np.fft.fft(samples)
-    kernel[1:W+1] = samples[W-1::-1] # first half, reversed
-    kernel_fft = np.fft.fft(kernel)
-    r_t_tau = np.fft.ifft(samples_fft * kernel_fft).real[W:]
-    # compute yin_t(tau)
-    yin = sqdiff - 2 * r_t_tau
-    return yin
-
-def sqd_yintapered(samples):
-    """ compute tappered sum of squared difference
-
-    Brute-force computation (cost o(N**2), slow)."""
-    B = len(samples)
-    W = B//2
-    yin = np.zeros(W)
-    for tau in range(1, W):
-        for j in range(W - tau):
-            yin[tau] += (samples[j] - samples[j+tau])**2
-    return yin
-
-def sqd_yinfft(samples):
-    """ compute yinfft modified sum of squared differences
-
-    Very fast, improved performance in transients.
-
-    FIXME: biased."""
-    B = len(samples)
-    W = B//2
-    yin = np.zeros(W)
-    def hanningz(W):
-        return .5 * (1. - np.cos(2. * np.pi * np.arange(W) / W))
-    #win = np.ones(B)
-    win = hanningz(B)
-    sqrmag = np.zeros(B)
-    fftout = np.fft.fft(win*samples)
-    sqrmag[0] = fftout[0].real**2
-    for l in range(1, W):
-        sqrmag[l] = fftout[l].real**2 + fftout[l].imag**2
-        sqrmag[B-l] = sqrmag[l]
-    sqrmag[W] = fftout[W].real**2
-    fftout = np.fft.fft(sqrmag)
-    sqrsum = 2.*sqrmag[:W + 1].sum()
-    yin[0] = 0
-    yin[1:] = sqrsum - fftout.real[1:W]
-    return yin / B
-
-def cumdiff(yin):
-    """ compute the cumulative mean normalized difference """
-    W = len(yin)
-    yin[0] = 1.
-    cumsum = 0.
-    for tau in range(1, W):
-        cumsum += yin[tau]
-        if cumsum != 0:
-            yin[tau] *= tau/cumsum
-        else:
-            yin[tau] = 1
-    return yin
-
-def compute_all(x):
-    import time
-    now = time.time()
-
-    yin     = sqd_yin(x)
-    t1 = time.time()
-    print ("yin took %.2fms" % ((t1-now) * 1000.))
-
-    yinfast = sqd_yinfast(x)
-    t2 = time.time()
-    print ("yinfast took: %.2fms" % ((t2-t1) * 1000.))
-
-    yintapered = sqd_yintapered(x)
-    t3 = time.time()
-    print ("yintapered took: %.2fms" % ((t3-t2) * 1000.))
-
-    yinfft  = sqd_yinfft(x)
-    t4 = time.time()
-    print ("yinfft took: %.2fms" % ((t4-t3) * 1000.))
-
-    return yin, yinfast, yintapered, yinfft
-
-def plot_all(yin, yinfast, yintapered, yinfft):
-    fig, axes = plt.subplots(nrows=2, ncols=2, sharex=True, sharey='col')
-
-    axes[0, 0].plot(yin, label='yin')
-    axes[0, 0].plot(yintapered, label='yintapered')
-    axes[0, 0].set_ylim(bottom=0)
-    axes[0, 0].legend()
-    axes[1, 0].plot(yinfast, '-', label='yinfast')
-    axes[1, 0].plot(yinfft, label='yinfft')
-    axes[1, 0].legend()
-
-    axes[0, 1].plot(cumdiff(yin), label='yin')
-    axes[0, 1].plot(cumdiff(yintapered), label='yin tapered')
-    axes[0, 1].set_ylim(bottom=0)
-    axes[0, 1].legend()
-    axes[1, 1].plot(cumdiff(yinfast), '-', label='yinfast')
-    axes[1, 1].plot(cumdiff(yinfft), label='yinfft')
-    axes[1, 1].legend()
-
-    fig.tight_layout()
-
-testfreqs = [441., 800., 10000., 40.]
-
-if len(sys.argv) > 1:
-    testfreqs = map(float,sys.argv[1:])
-
-for f in testfreqs:
-    print ("Comparing yin implementations for sine wave at %.fHz" % f)
-    samplerate = 44100.
-    win_s = 4096
-
-    x = np.cos(2.*np.pi * np.arange(win_s) * f / samplerate)
-
-    n_times = 1#00
-    for n in range(n_times):
-        yin, yinfast, yinfftslow, yinfft = compute_all(x)
-    if 0: # plot difference
-        plt.plot(yin-yinfast)
-        plt.tight_layout()
-        plt.show()
-    if 1:
-        plt.plot(yinfftslow-yinfft)
-        plt.tight_layout()
-        plt.show()
-    plot_all(yin, yinfast, yinfftslow, yinfft)
-plt.show()
diff --git a/python/demos/demo_yin_compare.py b/python/demos/demo_yin_compare.py
new file mode 100755 (executable)
index 0000000..92dc015
--- /dev/null
@@ -0,0 +1,175 @@
+#! /usr/bin/env python
+# -*- coding: utf8 -*-
+
+""" Pure python implementation of the sum of squared difference
+
+    sqd_yin: original sum of squared difference [0]
+        d_t(tau) = x ⊗ kernel
+    sqd_yinfast: sum of squared diff using complex domain [0]
+    sqd_yinfftslow: tappered squared diff [1]
+    sqd_yinfft: modified squared diff using complex domain [1]
+
+[0]:http://audition.ens.fr/adc/pdf/2002_JASA_YIN.pdf
+[1]:https://aubio.org/phd/
+"""
+
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+def sqd_yin(samples):
+    """ compute original sum of squared difference
+
+    Brute-force computation (cost o(N**2), slow)."""
+    B = len(samples)
+    W = B//2
+    yin = np.zeros(W)
+    for j in range(W):
+      for tau in range(1, W):
+        yin[tau] += (samples[j] - samples[j+tau])**2
+    return yin
+
+def sqd_yinfast(samples):
+    """ compute approximate sum of squared difference
+
+    Using complex convolution (fast, cost o(n*log(n)) )"""
+    # yin_t(tau) = (r_t(0) + r_(t+tau)(0)) - 2r_t(tau)
+    B = len(samples)
+    W = B//2
+    yin = np.zeros(W)
+    sqdiff = np.zeros(W)
+    kernel = np.zeros(B)
+    # compute r_(t+tau)(0)
+    squares = samples**2
+    for tau in range(W):
+        sqdiff[tau] = squares[tau:tau+W].sum()
+    # add r_t(0)
+    sqdiff += sqdiff[0]
+    # compute r_t(tau) using kernel convolution in complex domain
+    samples_fft = np.fft.fft(samples)
+    kernel[1:W+1] = samples[W-1::-1] # first half, reversed
+    kernel_fft = np.fft.fft(kernel)
+    r_t_tau = np.fft.ifft(samples_fft * kernel_fft).real[W:]
+    # compute yin_t(tau)
+    yin = sqdiff - 2 * r_t_tau
+    return yin
+
+def sqd_yintapered(samples):
+    """ compute tappered sum of squared difference
+
+    Brute-force computation (cost o(N**2), slow)."""
+    B = len(samples)
+    W = B//2
+    yin = np.zeros(W)
+    for tau in range(1, W):
+        for j in range(W - tau):
+            yin[tau] += (samples[j] - samples[j+tau])**2
+    return yin
+
+def sqd_yinfft(samples):
+    """ compute yinfft modified sum of squared differences
+
+    Very fast, improved performance in transients.
+
+    FIXME: biased."""
+    B = len(samples)
+    W = B//2
+    yin = np.zeros(W)
+    def hanningz(W):
+        return .5 * (1. - np.cos(2. * np.pi * np.arange(W) / W))
+    #win = np.ones(B)
+    win = hanningz(B)
+    sqrmag = np.zeros(B)
+    fftout = np.fft.fft(win*samples)
+    sqrmag[0] = fftout[0].real**2
+    for l in range(1, W):
+        sqrmag[l] = fftout[l].real**2 + fftout[l].imag**2
+        sqrmag[B-l] = sqrmag[l]
+    sqrmag[W] = fftout[W].real**2
+    fftout = np.fft.fft(sqrmag)
+    sqrsum = 2.*sqrmag[:W + 1].sum()
+    yin[0] = 0
+    yin[1:] = sqrsum - fftout.real[1:W]
+    return yin / B
+
+def cumdiff(yin):
+    """ compute the cumulative mean normalized difference """
+    W = len(yin)
+    yin[0] = 1.
+    cumsum = 0.
+    for tau in range(1, W):
+        cumsum += yin[tau]
+        if cumsum != 0:
+            yin[tau] *= tau/cumsum
+        else:
+            yin[tau] = 1
+    return yin
+
+def compute_all(x):
+    import time
+    now = time.time()
+
+    yin     = sqd_yin(x)
+    t1 = time.time()
+    print ("yin took %.2fms" % ((t1-now) * 1000.))
+
+    yinfast = sqd_yinfast(x)
+    t2 = time.time()
+    print ("yinfast took: %.2fms" % ((t2-t1) * 1000.))
+
+    yintapered = sqd_yintapered(x)
+    t3 = time.time()
+    print ("yintapered took: %.2fms" % ((t3-t2) * 1000.))
+
+    yinfft  = sqd_yinfft(x)
+    t4 = time.time()
+    print ("yinfft took: %.2fms" % ((t4-t3) * 1000.))
+
+    return yin, yinfast, yintapered, yinfft
+
+def plot_all(yin, yinfast, yintapered, yinfft):
+    fig, axes = plt.subplots(nrows=2, ncols=2, sharex=True, sharey='col')
+
+    axes[0, 0].plot(yin, label='yin')
+    axes[0, 0].plot(yintapered, label='yintapered')
+    axes[0, 0].set_ylim(bottom=0)
+    axes[0, 0].legend()
+    axes[1, 0].plot(yinfast, '-', label='yinfast')
+    axes[1, 0].plot(yinfft, label='yinfft')
+    axes[1, 0].legend()
+
+    axes[0, 1].plot(cumdiff(yin), label='yin')
+    axes[0, 1].plot(cumdiff(yintapered), label='yin tapered')
+    axes[0, 1].set_ylim(bottom=0)
+    axes[0, 1].legend()
+    axes[1, 1].plot(cumdiff(yinfast), '-', label='yinfast')
+    axes[1, 1].plot(cumdiff(yinfft), label='yinfft')
+    axes[1, 1].legend()
+
+    fig.tight_layout()
+
+testfreqs = [441., 800., 10000., 40.]
+
+if len(sys.argv) > 1:
+    testfreqs = map(float,sys.argv[1:])
+
+for f in testfreqs:
+    print ("Comparing yin implementations for sine wave at %.fHz" % f)
+    samplerate = 44100.
+    win_s = 4096
+
+    x = np.cos(2.*np.pi * np.arange(win_s) * f / samplerate)
+
+    n_times = 1#00
+    for n in range(n_times):
+        yin, yinfast, yinfftslow, yinfft = compute_all(x)
+    if 0: # plot difference
+        plt.plot(yin-yinfast)
+        plt.tight_layout()
+        plt.show()
+    if 1:
+        plt.plot(yinfftslow-yinfft)
+        plt.tight_layout()
+        plt.show()
+    plot_all(yin, yinfast, yinfftslow, yinfft)
+plt.show()
index 253ee64..7ee14d6 100644 (file)
@@ -81,6 +81,11 @@ extern "C" {
 
   see http://recherche.ircam.fr/equipes/pcm/pub/people/cheveign.html
 
+  \b \p yinfast: Yinfast algorithm
+
+  This algorithm is equivalent to the YIN algorithm, but computed in the
+  spectral domain for efficiency. See also `python/demos/demo_yin_compare.py`.
+
   \b \p yinfft : Yinfft algorithm
 
   This algorithm was derived from the YIN algorithm. In this implementation, a