Skip to content

Resampling

Sinc interpolation

The sample_quality settings from 3 to 10 are interpolations based on the windowed sinc method.

https://ccrma.stanford.edu/~jos/lumped/Windowed_Sinc_Interpolation.html

In the RGC sfz software, we have a number of quality settings available from a drop-down menu.
These sinc settings are labelled: 08 12 16 24 36 48 60 72
One can easily guess that these designate a number of points. The numbers are multiples of 4, which makes them practical for SIMD processing.

For implementation purposes, we want to precompute the sinc function in a large table. A window must be applied to the windowed sinc, in order to squash the edges. A Kaiser window is most appropriate, which permits to keep aliasing ripples under a determined amplitude threshold. It has a parameter Beta, permitting to establish a compromise: the higher Beta, the lower is the alias magnitude, but the selectivity of the "brick wall" filter worsens near the cutoff point.

At first it seems an alright choice of Beta may be from about 6 to 10, as table size increases from 8 to 72.

response.py: tool for plotting frequency responses of sinc interpolators dynamically, when size and Beta are varied, based on the deip.pdf paper.

Click here to view
#!/usr/bin/env python

import numpy as np
from scipy.signal import *
from scipy.fft import *
import matplotlib.pyplot as plt
import matplotlib.ticker as plticker
import matplotlib.figure as plfigure
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg, NavigationToolbar2Tk
import sys
#from cmd import Cmd
import tkinter as tk
from tkinter import ttk
import locale

# Set number of kernel points to compute
N = 1024

# Set type of interpolation
Interpolator = {'type': 'winsinc', 'points': 8, 'beta': 3.0}
#Interpolator = {'type': 'hermite3'}
#Interpolator = {'type': 'bspline3'}

# Nyquist frequency (only for plot, 1.0 for normalized in pi*rad/s)
NyquistF = 1.0
#NyquistF = 0.5 * 44100.0


def linear(x):
    y = 1.0-np.abs(x)
    y = np.maximum(0.0, y)
    return y

def bspline3(x):
    x = np.abs(x)
    x2 = x * x
    x3 = x2 * x
    y = 0
    p1 = (2./3.) - x2 + (1./2.) * x3
    p2 = (4./3.) - (2.) * x + x2 - (1./6.) * x3
    y = np.where(x < 2., p2, y)
    y = np.where(x < 1., p1, y)
    return y

def hermite3(x):
    x = np.abs(x)
    x2 = x * x
    x3 = x2 * x
    y = 0
    q = (5./2.) * x2
    p1 = (1.) - q + (3./2.) * x3
    p2 = (2.) - (4.) * x + q - (1./2.) * x3
    y = np.where(x < 2., p2, y)
    y = np.where(x < 1., p1, y)
    return y

def winsinc(x, beta):
    y = np.where(x == 0, 1.0, np.sin(np.pi * x) / (np.pi * x))
    y *= kaiser(len(x), beta)
    return y

###
def kernel(itp):
    type = itp['type']
    if type == 'hermite3':
        Ex = 4
        X = np.linspace(-Ex/2.0, Ex/2.0, N)
        Y = hermite3(X)
    elif type == 'bspline3':
        Ex = 4
        X = np.linspace(-Ex/2.0, Ex/2.0, N)
        Y = bspline3(X)
    elif type == 'linear':
        Ex = 2
        X = np.linspace(-Ex/2.0, Ex/2.0, N)
        Y = linear(X)
    elif type == 'winsinc':
        Ex = int(itp['points'])
        X = np.linspace(-Ex/2.0, Ex/2.0, N)
        Y = winsinc(X, float(itp['beta']))
    else:
        raise ValueError('Unknown type of interpolation')
    return Ex, X, Y

###
def plot(itp):
    print('Interpolator: %s' % (Interpolator))
    try:
        Ex, X, Y = kernel(Interpolator)
    except ValueError as e:
        print('Error:', str(e))
        return
    W, H = freqz(Y, worN=64*N, fs=1.0)

    fig = plt.gcf()
    fig.clf()
    fig.set_figwidth(15)
    ax1, ax2 = fig.subplots(2)

    ax1.grid(alpha=0.25)
    ax1.plot(X, Y)
    ax1.set_xlim(-Ex/2.0, Ex/2.0)
    ax1.xaxis.set_major_locator(plticker.MultipleLocator(base=1.0))

    ax2.grid(alpha=0.25)
    ax2.plot(NyquistF*W*N*2.0/Ex, 20*np.log10(np.abs(H)/np.sum(Y)))
    ax2.set_xlim(0, 5*NyquistF)
    ax2.set_ylim(-120, 12)
    ax2.xaxis.set_major_locator(plticker.MultipleLocator(base=NyquistF))
    ax2.yaxis.set_major_locator(plticker.MultipleLocator(base=12.0))
    if NyquistF == 1.0:
        ax2.set_xlabel('Frequency (π rad/s)')
    else:
        ax2.set_xlabel('Frequency (Hz)')
    ax2.set_ylabel('Magnitude (dB)')

    plt.draw()

def main(args):
    plot(Interpolator)

    class Application(tk.Frame):
        def __init__(self, master=None):
            super().__init__(master)
            self.master = master
            self.pack()
            self.create_widgets()

        def create_widgets(self):
            CHOICES = ['linear', 'bspline3', 'hermite3', 'winsinc']

            self.choice = ttk.Combobox(self, values=CHOICES)
            self.choice.current(CHOICES.index(Interpolator['type']))
            self.choice.bind("<<ComboboxSelected>>", self.on_change_type)
            self.choice.pack(side="top")

            self.beta = tk.Scale(self.master, from_=1.0, to=16.0, digits=3, resolution=0.01, orient=tk.HORIZONTAL, command=self.on_change_beta)
            self.beta.set(Interpolator['beta'])
            self.beta.pack(side="top")

            self.numpoints = tk.Scale(self.master, from_=8, to=72, resolution=4, orient=tk.HORIZONTAL, command=self.on_change_numpoints)
            self.numpoints.set(Interpolator['points'])
            self.numpoints.pack(side="top")

            self.var_normalize = tk.IntVar()
            self.normalize = tk.Checkbutton(self.master, text='Normalize', command=self.on_change_normalized, variable=self.var_normalize)
            if NyquistF == 1.0:
                self.normalize.select()
            else:
                self.normalize.deselect()
            self.normalize.pack(side="top")

            # self.quit = tk.Button(self, text="QUIT", fg="red", command=self.master.destroy)
            # self.quit.pack(side="bottom")

            self.canvas = FigureCanvasTkAgg(plt.gcf(), master=self.master)
            self.canvas.get_tk_widget().pack(side="bottom")

        def on_change_type(self, event):
            Interpolator['type'] = self.choice.get()
            plot(Interpolator)

        def on_change_beta(self, value):
            Interpolator['beta'] = float(value)
            plot(Interpolator)

        def on_change_numpoints(self, value):
            Interpolator['points'] = int(value)
            plot(Interpolator)

        def on_change_normalized(self):
            value = self.var_normalize.get()
            global NyquistF
            NyquistF = (value == 0) and (0.5 * 44100) or 1.0
            plot(Interpolator)

    locale.setlocale(locale.LC_NUMERIC, 'C')
    root = tk.Tk()
    app = Application(master=root)
    app.mainloop()

###
if __name__ == '__main__':
    main(sys.argv)