pyfar.dsp.InterpolateSpectrum overshoots

[f-brinkmann]

General

• pyfar version: 0.6.0

Description

pyfar.dsp.InterpolateSpectrum overshoots if not using linear or nearest neigbor interpolation. This is because we internally use scipy.interpolate.interp1d which is a non-monotonic interpolation for higher plynomial orders. (see plot and code below)

Further Issues

• We are using interp1d, which is not maintained anymore.
• Alternative soloutions have non-unique interfaces
  • PChip and CubicSpline take multideimensional input. interp and UnivariateSpline dont.
  • interp does not have an argument for extrapolation.
  • UnivariateSpline can not handle complex input

Suggested solution

• Maintain backwards compatability until pyfar 0.9.0
• New interpolation methods are monotonic, non-monotonic and smoothing splines.
• Linear and nearest neigbor interpolation will be kept from interp1d but will be deprecated if interp1d gets deprecated by scipy
• Extrapolation and smoothing parameters will be moved to the call to increase flexibiliy

Example

import pyfar as pf
import numpy as np
from scipy.interpolate import PchipInterpolator
from scipy.signal import firls
pf.plot.use()

N = 2049                         # number of samples in interpolates Signal
fs = 44100                       # sampling rate
f = pf.dsp.fft.rfftfreq(N, fs)   # frequencies

# pyfar interpolation causes overshoots for everything that is not nearest
# neighbor or linear interpolation (uses scipy.interpolate.interp1d)
data = pf.FrequencyData([1, 10, 10, .1, .1, 1], [0, 1e3, 2e3, 4e3, 8e3, fs/2])

interpolator = pf.dsp.InterpolateSpectrum(
    data, 'magnitude', ('nearest', 'quadratic', 'nearest'),
    fscale='linear')

signal = interpolator(2048, 44100)

# interpolation with scipy.interpolate.PchipInterpolator
interpolator_pchip = PchipInterpolator(
    data.frequencies, np.abs(data.freq).flatten(), extrapolate=False)

signal_pchip = pf.Signal(
    interpolator_pchip(f), fs, N, 'freq')

# least-squares filter design
signal_ls = pf.Signal(
    firls(N, data.frequencies, np.abs(data.freq).flatten(), fs=fs), fs)

# plot
dB = True
ax = pf.plot.freq(signal, dB=dB, label='pyfar')
pf.plot.freq(signal_pchip, dB=dB, label='monotonic cubic spline')
pf.plot.freq(signal_ls, dB=dB, label='least squares FIR')
pf.plot.freq(data, dB=dB, ax=ax, c='g', ls='', marker='.', label='input')
ax.legend()