[FEATURE] fractional octave smoothing in plot_freq, and plot_phase

pytta/__init__.py

@@ -113,7 +113,7 @@
     list_devices, print_devices, get_device_from_user,\
     fft_convolve, find_delay, corr_coef, resample,\
     peak_time, save, load, plot_time, plot_time_dB, plot_freq,\
-    plot_bars, plot_spectrogram, plot_waterfall, SPL
+    plot_bars, plot_spectrogram, plot_waterfall, SPL, smooth_frequency, plot_phase
 
 # Must go on v0.1.0 (DEPRECATED)
 from .functions import fft_degree
@@ -146,12 +146,14 @@
            'resample',
            'corr_coef',
            'peak_time',
+           'smooth_frequency',
            'save',
            'load',
            'weighting',
            'plot_time',
            'plot_time_dB',
            'plot_freq',
+           'plot_phase',
            'plot_bars',
            'plot_spectrogram',
            'plot_waterfall',

pytta/_plot.py

@@ -339,7 +339,8 @@ def _curve_data_extractor_time_dB(sigObjs):
     return curveData
 
 
-def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep):
+def freq(sigObjs, xLabel, yLabel, yLim, xLim, title, decimalSep,
+         phase=False, unwrap=True):
     """
     Plots a signal decibel magnitude in frequency domain.
 
@@ -349,9 +350,6 @@ def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep):
        * sigObjs (), (list):
             a list with SignalObjs.
 
-        * smooth (False), (bool):
-            option for curve smoothing. Uses scipy.signal.savgol_filter.
-            Preliminar implementation. Needs review.
 
         * xLabel ('Time [s]'), (str):
             x axis label.
@@ -376,6 +374,10 @@ def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep):
             may be dot or comma.
 
             >>> decimalSep = ',' # in Brazil
+        * phase (False), (bool):
+            plot the phase in degrees instead of the magnitude in dB.
+        * unwrap (True), (bool):
+            unwrap the phase over frequency.
 
     Return:
     --------
@@ -394,7 +396,7 @@ def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep):
     if xLabel is None:
         xLabel = 'Frequency [Hz]'
     if yLabel is None:
-        yLabel = 'Magnitude [dB]'
+        yLabel = 'Phase [deg]' if phase else 'Magnitude [dB]'
     if title is None:
         title = ''
 
@@ -409,11 +411,13 @@ def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep):
 
     curveData = _curve_data_extractor_freq(sigObjs)
     for data in curveData:
-        with np.errstate(divide='ignore'):
-            dBSignal = 20*np.log10(np.abs(data['y'])/data['dBRef'])
+        if phase:
+            angle = np.angle(data['y'])
+            dBSignal = np.degrees(np.unwrap(angle) if unwrap else angle)
+        else:
+            with np.errstate(divide='ignore'):
+                dBSignal = 20*np.log10(np.abs(data['y'])/data['dBRef'])
 
-        if smooth:
-            dBSignal = ss.savgol_filter(dBSignal, 31, 3)
 
         ax.semilogx(data['x'], dBSignal, label=data['label'])
 

pytta/classes/signal.py

@@ -677,9 +677,10 @@ def plot_freq(self, smooth:bool=False, xLabel:str=None, yLabel:str=None,
         Parameters (default), (type):
         -----------------------------
 
-            * smooth (False), (bool):
-                option for curve smoothing. Uses scipy.signal.savgol_filter.
-                Preliminar implementation. Needs review.
+            * smooth (False), (bool or number):
+                fractional octave smoothing of the magnitude by
+                pytta.smooth_frequency: True for 1/3 octave, a number n for
+                1/n octave;
 
             * xLabel ('Time [s]'), (str):
                 x axis label.
@@ -731,7 +732,8 @@ def plot_freq(self, smooth:bool=False, xLabel:str=None, yLabel:str=None,
                 if self.freqTitle is not None:
                     title = self.freqTitle
 
-        fig = plot.freq((self,), smooth, xLabel, yLabel, yLim, xLim, title,
+        from pytta.functions import _smoothed  # pytta.functions imports this module
+        fig = plot.freq(_smoothed((self,), smooth), xLabel, yLabel, yLim, xLim, title,
                         decimalSep)
         return fig
 
@@ -1147,7 +1149,7 @@ def _ifft(self):
             adjustedFreqSignal = np.zeros(self._freqSignal.shape, dtype=np.complex_)
             # * N/2 for AC
             adjustedFreqSignal[1:, :] = \
-                self._freqSignal * len(self._freqSignal)
+                self._freqSignal[1:, :] * len(self._freqSignal)
             # * N for DC
             adjustedFreqSignal[0, :] = \
                 self._freqSignal[0, :] * self.numSamples

pytta/functions.py

@@ -17,9 +17,11 @@
     >>> pytta.corr_coef(signalObj1, signalObj2)
     >>> pytta.resample(signalObj, newSamplingRate)
     >>> pytta.peak_time(signalObj1, signalObj2, ..., signalObjN)
+    >>> pytta.smooth_frequency(signalObj, nthOct)
     >>> pytta.plot_time(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.plot_time_dB(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.plot_freq(signalObj1, signalObj2, ..., signalObjN)
+    >>> pytta.plot_phase(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.plot_bars(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.save(fileName, obj1, ..., objN)
     >>> pytta.load(fileName)
@@ -37,6 +39,7 @@
 import sounddevice as sd
 import scipy.signal as ss
 import scipy.fftpack as sfft
+from scipy.interpolate import CubicSpline
 import zipfile as zf
 import h5py
 from typing import Union, List
@@ -274,6 +277,68 @@ def peak_time(signal):
     else:
         return peaks_time[0]
 
+def smooth_frequency(signalObj, nthOct=3):
+    """
+    Smooth the magnitude spectrum of a SignalObj over fractions of an octave
+    and return a new SignalObj. The phase is kept.
+
+    The magnitude of freqSignal is interpolated linearly onto a logarithmic
+    frequency grid of 1 cent that spans the 20 octaves below Nyquist, box
+    averaged over 1200 / nthOct cents (made odd; the box shrinks
+    symmetrically at both ends of the grid), and interpolated back onto the
+    frequency bins with a cubic spline. The DC bin is copied unchanged.
+
+    Parameters (default), (type):
+    -----------------------------
+        * signalObj (), (SignalObj):
+            signal whose spectrum is smoothed. It is not modified;
+        * nthOct (3), (int or float):
+            smoothing width of 1 / nthOct octave: 1 for one octave, 3 for a
+            third, 12 for a twelfth.
+
+    The smoothing acts on freqSignal as stored, so the DC and Nyquist bins,
+    which PyTTa scales differently from the others, take part with that
+    scale. This is the 'LogFreqOctave1' smoothing of ita_smooth in the
+    ITA-Toolbox (RWTH Aachen).
+
+    >>> smoothed = pytta.smooth_frequency(ir, nthOct=3)
+    """
+    if not isinstance(signalObj, SignalObj):
+        raise TypeError('signalObj must be a SignalObj.')
+    if isinstance(nthOct, bool) or not isinstance(nthOct, (int, float)) \
+            or nthOct <= 0:
+        raise ValueError('nthOct must be a positive number.')
+    samplingRate = signalObj.samplingRate
+    freq = signalObj.freqVector
+    spectrum = np.array(signalObj.freqSignal, dtype=complex, ndmin=2)
+    # logarithmic grid: 1 cent steps over the 20 octaves below Nyquist
+    gridFreq = (samplingRate / 2) * 2 ** (np.arange(24001) / 1200 - 20)
+    span = int(np.floor(1200 / nthOct + 0.5))
+    width = min(span - 1 + span % 2, gridFreq.size)
+    halfWidth = (width - 1) // 2
+    index = np.arange(gridFreq.size)
+    reach = np.minimum(halfWidth, np.minimum(index, index[::-1]))
+    smoothed = np.empty_like(spectrum)
+    for ch in range(spectrum.shape[1]):
+        magnitude = np.interp(gridFreq, freq, np.abs(spectrum[:, ch]))
+        summed = np.concatenate(([0.], np.cumsum(magnitude)))
+        average = (summed[index + reach + 1] - summed[index - reach]) \
+            / (2 * reach + 1)
+        spline = CubicSpline(gridFreq[1:], average[1:], bc_type='not-a-knot')
+        smoothed[0, ch] = spectrum[0, ch]
+        smoothed[1:, ch] = spline(freq[1:]) \
+            * np.exp(1j * np.angle(spectrum[1:, ch]))
+    newSignal = SignalObj(smoothed, 'freq',
+                          samplingRate=samplingRate,
+                          numSamples=signalObj.numSamples,
+                          signalType=signalObj.signalType,
+                          freqMin=signalObj.freqMin,
+                          freqMax=signalObj.freqMax,
+                          comment=signalObj.comment)
+    newSignal.channels = cp.deepcopy(signalObj.channels)
+    return newSignal
+
+
 def fft_degree(*args,**kwargs):
     """
     DEPRECATED
@@ -423,6 +488,9 @@ def plot_freq(*sigObjs, smooth:bool=False, xLabel:str=None, yLabel:str=None,
         * sigObjs (), (SignalObj):
             non-keyworded input arguments with N SignalObjs.
 
+        * smooth (False), (bool or number):
+            fractional octave smoothing of the magnitude by smooth_frequency:
+            True for 1/3 octave, a number n for 1/n octave;
         * xLabel ('Time [s]'), (str):
             x axis label.
 
@@ -455,12 +523,44 @@ def plot_freq(*sigObjs, smooth:bool=False, xLabel:str=None, yLabel:str=None,
     realSigObjs = \
         _remove_non_(SignalObj, sigObjs, msgPrefix='plot_freq:')
     if len(realSigObjs) > 0:
-        fig = plot.freq(realSigObjs, smooth, xLabel, yLabel, yLim, xLim, title,
-                        decimalSep)
+        fig = plot.freq(_smoothed(realSigObjs, smooth), xLabel, yLabel, yLim,
+                        xLim, title, decimalSep)
         return fig
     else:
         return
 
+def _smoothed(sigObjs, smooth):
+    """The signals, smoothed by smooth_frequency when smooth asks for it:
+    True is 1/3 octave, a number n is 1/n octave."""
+    if smooth is False or smooth is None:
+        return list(sigObjs)
+    nthOct = 3 if smooth is True else smooth
+    return [smooth_frequency(sig, nthOct) for sig in sigObjs]
+
+
+def plot_phase(*sigObjs, unwrap: bool = True, xLabel: str = None,
+               yLabel: str = None, yLim: list = None, xLim: list = None,
+               title: str = None, decimalSep: str = ','):
+    """
+    Plot the phase spectrum of the provided SignalObjs, in degrees, on a
+    logarithmic frequency axis.
+
+    Parameters (default), (type):
+    -----------------------------
+        * sigObjs (), (SignalObj): non-keyworded SignalObjs;
+        * unwrap (True), (bool): unwrap the phase over frequency;
+        * xLabel, yLabel, yLim, xLim, title, decimalSep: as in plot_freq.
+
+    Return:
+    --------
+        matplotlib.figure.Figure object.
+    """
+    realSigObjs = _remove_non_(SignalObj, sigObjs, msgPrefix='plot_phase:')
+    if len(realSigObjs) > 0:
+        return plot.freq(realSigObjs, xLabel, yLabel, yLim, xLim, title,
+                         decimalSep, phase=True, unwrap=unwrap)
+
+
 def plot_bars(*analyses, xLabel:str=None, yLabel:str=None,
               yLim:list=None, xLim:list=None, title:str=None, decimalSep:str=',',
               barWidth:float=0.75, errorStyle:str=None,

tests/golden/make_smooth_frequency_ita.m

@@ -0,0 +1,40 @@
+function make_smooth_frequency_ita(itaRoot)
+% MAKE_SMOOTH_FREQUENCY_ITA  Golden data for tests/test_smooth_frequency.py.
+%
+%   make_smooth_frequency_ita(itaRoot) runs ita_smooth_frequency from the
+%   ITA-Toolbox (RWTH Aachen, https://git.rwth-aachen.de/ita/toolbox) on a
+%   fixed random signal and writes smooth_frequency_ita.mat next to this
+%   file: the time signal, the ITA input spectrum and one output spectrum
+%   per bandwidth (1, 1/3 and 1/12 octave).
+%
+%   The DC and Nyquist bins of the signal are set to zero. PyTTa and the
+%   ITA-Toolbox scale these two bins differently from the others, so with
+%   them at zero the two spectra differ by one constant factor and the
+%   comparison pins the smoothing alone.
+%
+%   Usage:  matlab -batch "make_smooth_frequency_ita('/path/to/ITA-Toolbox')"
+
+addpath(genpath(itaRoot));
+
+fs = 48000;
+n = 2048;
+rng(0, 'twister');
+x = randn(n, 2);
+alt = (-1) .^ (0:n-1)';
+x = x - mean(x, 1);                      % DC bin to zero
+x = x - alt .* mean(x .* alt, 1);        % Nyquist bin to zero
+a = itaAudio(x, fs, 'time');
+
+nthOct = [1 3 12];
+cases = struct('nthOct', {}, 'freqData', {});
+for k = 1:numel(nthOct)
+    r = ita_smooth_frequency(a, 'bandwidth', 1 / nthOct(k));
+    cases(k).nthOct   = nthOct(k);
+    cases(k).freqData = r.freqData;
+end
+freqDataIn = a.freqData;
+
+out = fullfile(fileparts(mfilename('fullpath')), 'smooth_frequency_ita.mat');
+save(out, 'x', 'fs', 'freqDataIn', 'cases', '-v7');
+fprintf('wrote %s (%d cases)\n', out, numel(cases));
+end

tests/test_plot_freq.py

@@ -0,0 +1,44 @@
+"""
+Tests for the curves drawn by pytta.plot_freq with smoothing and by
+pytta.plot_phase. Each test pins one rule on the data of the plotted line.
+"""
+import os
+import unittest
+
+os.environ.setdefault('MPLBACKEND', 'Agg')
+
+import numpy as np
+
+import pytta
+
+FS = 48000
+
+
+class TestPlotFreq(unittest.TestCase):
+
+    def test_smoothing_is_fractional_octave(self):
+        # smooth=n draws the magnitude of smooth_frequency(signal, n), and
+        # smooth=True is a third of an octave.
+        rng = np.random.default_rng(0)
+        sig = pytta.SignalObj(rng.standard_normal(4096), 'time', FS)
+        expected = 20 * np.log10(np.abs(
+            pytta.smooth_frequency(sig, 3).freqSignal[:, 0]))
+        for smooth in (3, True):
+            fig = pytta.plot_freq(sig, smooth=smooth, decimalSep='.')
+            drawn = fig.axes[0].lines[0].get_ydata()
+            self.assertTrue(np.allclose(drawn[1:], expected[1:], atol=1e-4))
+            fig.clf()
+
+    def test_phase_of_a_delay(self):
+        # A delay of d samples has the phase -360 f d / fs degrees.
+        d = 12
+        sig = pytta.SignalObj(np.eye(1, 4096, d)[0], 'time', FS)
+        fig = pytta.plot_phase(sig, decimalSep='.')
+        line = fig.axes[0].lines[0]
+        expected = -360 * line.get_xdata() * d / FS
+        self.assertTrue(np.allclose(line.get_ydata(), expected, atol=1e-2))
+        fig.clf()
+
+
+if __name__ == '__main__':
+    unittest.main()