ref: 48fa81d3de6e167612b3d2bde13af325052d0ba5
parent: c268bf2c950b5e3243fa03f15c2b4f1ef095175a
author: Paul Brossier <piem@piem.org>
date: Tue May 9 12:06:23 EDT 2017
src/pitch/pitch.h: add a note about yinfast
--- a/python/demos/demo_fastyin_compare.py
+++ /dev/null
@@ -1,175 +1,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()
--- /dev/null
+++ b/python/demos/demo_yin_compare.py
@@ -1,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()
--- a/src/pitch/pitch.h
+++ b/src/pitch/pitch.h
@@ -81,6 +81,11 @@
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