Skip to content

Commit

Permalink
Merge pull request #453 from pyfar/hotfix/interpolate_spectrum
Browse files Browse the repository at this point in the history 
Hotfix: Spectrum interpolation on logarithmic frequency axis avoiding zero frequency.
  • Loading branch information
@mberz
mberz committed Mar 30, 2023
2 parents bb8f9a6 + c1661de commit 85739c9
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 84 deletions.
111 changes: 54 additions & 57 deletions pyfar/dsp/interpolation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -640,7 +640,7 @@ class InterpolateSpectrum():
Parameters
----------
data : FrequencyData
Input data to be interpolated. `data.fft_norm` must be `'none'`.
Input data to be interpolated.
method : string
Specifies the input data for the interpolation
Expand Down Expand Up @@ -678,10 +678,7 @@ class InterpolateSpectrum():
``'linear'``
Interpolate on a linear frequency axis.
``'log'``
Interpolate on a logarithmic frequency axis. Note that 0 Hz can
not be interpolated on a logarithmic scale because the logarithm
of 0 does not exist. Frequencies of 0 Hz are thus replaced by the
next highest frequency before interpolation.
Interpolate on a logarithmic frequency axis.
The default is ``'linear'``.
clip : bool, tuple
Expand Down Expand Up @@ -717,39 +714,43 @@ class InterpolateSpectrum():
>>> import pyfar as pf
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>>
>>> pf.plot.use()
>>> _, ax = plt.subplots(2, 3)
>>>
>>> # generate data
>>> data = pf.FrequencyData([1, 0], [5e3, 20e3])
>>> interpolator = pf.dsp.InterpolateSpectrum(
... data, 'magnitude', ('nearest', 'linear', 'nearest'))
>>> # interpolate 64 samples at a sampling rate of 44100
>>> signal = interpolator(64, 44100)
>>> # add linear phase
>>> signal = pf.dsp.linear_phase(signal, 32)
>>> # plot input and output data
>>> with pf.plot.context():
>>> _, ax = plt.subplots(2, 2)
>>>
>>> # interpolate and plot
>>> for ff, fscale in enumerate(["linear", "log"]):
>>> interpolator = pf.dsp.InterpolateSpectrum(
... data, 'magnitude', ('nearest', 'linear', 'nearest'),
... fscale)
>>>
>>> # interpolate to 64 samples linear phase impulse response
>>> signal = interpolator(64, 44100)
>>> signal = pf.dsp.linear_phase(signal, 32)
>>>
>>> # time signal (linear and logarithmic amplitude)
>>> pf.plot.time(signal, ax=ax[0, 0])
>>> pf.plot.time(signal, ax=ax[1, 0], dB=True)
>>> pf.plot.time(signal, ax=ax[ff, 0], unit='ms', dB=True)
>>> # frequency plot (linear x-axis)
>>> pf.plot.freq(signal, dB=False, freq_scale="linear",
... ax=ax[0, 1])
>>> pf.plot.freq(
... signal, dB=False, freq_scale="linear", ax=ax[ff, 1])
>>> pf.plot.freq(data, dB=False, freq_scale="linear",
... ax=ax[0, 1], c='r', ls='', marker='.')
>>> ax[0, 1].set_xlim(0, signal.sampling_rate/2)
... ax=ax[ff, 1], c='r', ls='', marker='.')
>>> ax[ff, 1].set_xlim(0, signal.sampling_rate/2)
>>> ax[ff, 1].set_title(
... f"Interpolated on {fscale} frequency scale")
>>> # frequency plot (log x-axis)
>>> pf.plot.freq(signal, dB=False, ax=ax[1, 1], label='input')
>>> pf.plot.freq(data, dB=False, ax=ax[1, 1],
>>> pf.plot.freq(signal, dB=False, ax=ax[ff, 2], label='input')
>>> pf.plot.freq(data, dB=False, ax=ax[ff, 2],
... c='r', ls='', marker='.', label='output')
>>> min_freq = np.min([signal.sampling_rate / signal.n_samples,
... data.frequencies[0]])
>>> ax[1, 1].set_xlim(min_freq, signal.sampling_rate/2)
>>> ax[1, 1].legend(loc='best')
>>> ax[ff, 2].set_xlim(2e3, signal.sampling_rate/2)
>>> ax[ff, 2].legend(loc='best')
"""

def __init__(self, data, method, kind, fscale='linear',
clip=False, group_delay=None, unit='samples'):
def __init__(self, data, method, kind, fscale='linear', clip=False):

# check input ---------------------------------------------------------
# ... data
Expand Down Expand Up @@ -789,6 +790,7 @@ def __init__(self, data, method, kind, fscale='linear',
self._method = method
self._clip = clip
self._fscale = fscale
self._kind = kind

# flatten input data to work with scipy interpolators
self._cshape = data.cshape
Expand All @@ -805,33 +807,42 @@ def __init__(self, data, method, kind, fscale='linear',
self._data = [np.abs(data.freq)]

# frequencies for interpolation (store for testing)
self._f_in = self._get_frequencies(data.frequencies.copy())
self._f_in = data.frequencies.copy()

def __call__(self, n_samples, sampling_rate, show=False):
"""
Interpolate a Signal with n_samples length.
(see class docstring) for more information.
"""

# length of half sided spectrum and highest frequency
n_fft = n_samples//2 + 1
f_max = sampling_rate / n_samples * (n_fft - 1)
# get the frequency values
if self._fscale == "linear":
# linearly spaced frequencies
self._f_query = pf.dsp.fft.rfftfreq(n_samples, sampling_rate)
self._f_base = self._f_in
else:
# logarithmically scaled frequencies between 0 and log10(n_fft)
self._f_query = np.log10(np.arange(1, n_fft+1))
self._f_base = np.log10(self._f_in / f_max * (n_fft - 1) + 1)

# frequency range
self._freq_range = [self._f_in[0], self._f_in[-1]]
self._freq_range = [self._f_base[0], self._f_base[-1]]

# get the interpolators
self._interpolators = []
for d in self._data:
interpolators = []
for idx, k in enumerate(kind):
for idx, k in enumerate(self._kind):
if idx == 1:
interpolators.append(interp1d(self._f_in, d, k))
interpolators.append(interp1d(self._f_base, d, k))
else:
interpolators.append(interp1d(
self._f_in, d, k, fill_value="extrapolate"))
self._f_base, d, k, fill_value="extrapolate"))
self._interpolators.append(interpolators)

def __call__(self, n_samples, sampling_rate, show=False):
"""
Interpolate a Signal with n_samples length.
(see class docstring) for more information.
"""

# get the query frequencies (store for testing)
self._f_query = self._get_frequencies(
pf.dsp.fft.rfftfreq(n_samples, sampling_rate))

# get interpolation ranges
id_below = self._f_query < self._freq_range[0]
id_within = np.logical_and(self._f_query >= self._freq_range[0],
Expand Down Expand Up @@ -890,17 +901,3 @@ def __call__(self, n_samples, sampling_rate, show=False):
ax[1, 1].legend(loc='best')

return signal

def _get_frequencies(self, frequencies):
"""
Return frequencies for creating or quering interpolation objects.
In case logfrequencies are requested, 0 Hz entries are replaced by
the next highest frequency, because the logarithm of 0 does not exist.
"""
if self._fscale == "log":
if frequencies[0] == 0:
frequencies[0] = frequencies[1]
frequencies = np.log(frequencies)

return frequencies
50 changes: 23 additions & 27 deletions tests/test_dsp_interpolation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,37 +325,33 @@ def test_interpolate_spectrum_clip():
np.all(np.abs(signal_clip.freq) <= 2)


def test_interpolate_spectrum_fscale():
@pytest.mark.parametrize(
'fscale,n_samples,sampling_rate,f_in,f_base,f_query',
[('linear', 10, 40, [0, 10, 20], [0, 10, 20],
pf.dsp.fft.rfftfreq(10, 40)),
('log', 10, 40, [0, 10, 20], [0, .544, .778],
[0, .301, .477, .602, .699, .778])])
def test_interpolate_spectrum_fscale(
fscale, n_samples, sampling_rate, f_in, f_base, f_query):
"""
Test frequency vectors for linear and logarithmic frequency interpolation.
"""

# test parametres and data
f_in_lin = [0, 10, 20]
f_in_log = np.log([10, 10, 20])
n_samples = 10
sampling_rate = 40
f_query_lin = pf.dsp.fft.rfftfreq(n_samples, sampling_rate)
f_query_log = f_query_lin.copy()
f_query_log[0] = f_query_log[1]
f_query_log = np.log(f_query_log)
data = pf.FrequencyData([1, 1, 1], f_in_lin)

# generate interpolator with linear frequency
interpolator_lin = InterpolateSpectrum(
data, "magnitude", ("linear", "linear", "linear"), fscale="linear")
_ = interpolator_lin(n_samples, sampling_rate)
# generate interpolator with logarithmic frequency
interpolator_log = InterpolateSpectrum(
data, "magnitude", ("linear", "linear", "linear"), fscale="log")
_ = interpolator_log(n_samples, sampling_rate)

# test frequency vectors
npt.assert_allclose(interpolator_lin._f_in, f_in_lin)
npt.assert_allclose(interpolator_lin._f_query, f_query_lin)

npt.assert_allclose(interpolator_log._f_in, f_in_log)
npt.assert_allclose(interpolator_log._f_query, f_query_log)
# test data
data = pf.FrequencyData([1, 1, 1], f_in)

# generate interpolator
interpolator = InterpolateSpectrum(
data, "magnitude", ("linear", "linear", "linear"), fscale=fscale)
signal = interpolator(n_samples, sampling_rate)

# test time and frequency vectors
npt.assert_allclose(signal.time, pf.signals.impulse(n_samples).time)
npt.assert_almost_equal(interpolator._f_in, f_in)
npt.assert_almost_equal(interpolator._f_base, f_base, 3)
npt.assert_almost_equal(interpolator._f_query, f_query, 3)
npt.assert_almost_equal(interpolator._f_query[[0, -1]],
interpolator._f_base[[0, -1]])


def test_interpolate_spectrum_show():
Expand Down

0 comments on commit 85739c9

Please sign in to comment.