From: Paul Brossier Date: Tue, 9 May 2017 14:06:23 +0000 (+0200) Subject: src/pitch/pitch.h: add a note about yinfast X-Git-Tag: 0.4.6~51^2~6 X-Git-Url: https://git.aubio.org/?a=commitdiff_plain;h=48fa81d3de6e167612b3d2bde13af325052d0ba5;p=aubio.git src/pitch/pitch.h: add a note about yinfast --- diff --git a/python/demos/demo_fastyin_compare.py b/python/demos/demo_fastyin_compare.py deleted file mode 100755 index 92dc015d..00000000 --- a/python/demos/demo_fastyin_compare.py +++ /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 index 00000000..92dc015d --- /dev/null +++ b/python/demos/demo_yin_compare.py @@ -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() diff --git a/src/pitch/pitch.h b/src/pitch/pitch.h index 253ee648..7ee14d63 100644 --- a/src/pitch/pitch.h +++ b/src/pitch/pitch.h @@ -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