diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..5dbff89 100644 --- a/pytta/__init__.py +++ b/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', diff --git a/pytta/_plot.py b/pytta/_plot.py index c0a0970..1036651 100644 --- a/pytta/_plot.py +++ b/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']) diff --git a/pytta/classes/signal.py b/pytta/classes/signal.py index 79151c3..b7a46ac 100644 --- a/pytta/classes/signal.py +++ b/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 diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..8465f45 100644 --- a/pytta/functions.py +++ b/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, diff --git a/tests/golden/make_smooth_frequency_ita.m b/tests/golden/make_smooth_frequency_ita.m new file mode 100644 index 0000000..8bf226c --- /dev/null +++ b/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 diff --git a/tests/golden/smooth_frequency_ita.mat b/tests/golden/smooth_frequency_ita.mat new file mode 100644 index 0000000..5b8a80d Binary files /dev/null and b/tests/golden/smooth_frequency_ita.mat differ diff --git a/tests/test_plot_freq.py b/tests/test_plot_freq.py new file mode 100644 index 0000000..c07e4b7 --- /dev/null +++ b/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() diff --git a/tests/test_smooth_frequency.py b/tests/test_smooth_frequency.py new file mode 100644 index 0000000..60842f8 --- /dev/null +++ b/tests/test_smooth_frequency.py @@ -0,0 +1,144 @@ +""" +Tests for pytta.smooth_frequency. + +Each test pins one rule. The spectra are built directly in the frequency +domain, so the input magnitude is known bin by bin. The last test compares +the whole thing with the ITA-Toolbox (RWTH Aachen) on data generated by +golden/make_smooth_frequency_ita.m. + +The smoothing is a moving average over a logarithmic frequency grid of +1 cent, 1200 / nthOct cents wide, applied to the magnitude of freqSignal; +the phase is kept. The grid runs 20 octaves below Nyquist, so the first +bins take part through linear interpolation from DC. +""" +import os +import unittest + +import numpy as np +import scipy.io as sio + +import pytta + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden', + 'smooth_frequency_ita.mat') +FS = 48000 +N = 2048 +K = N // 2 + 1 +TOL32 = 1e-6 # rounding of the complex64 freqSignal storage + + +# SignalObj built from a half spectrum of K rows. Built inline in each test: +# SignalObj inspects the caller's frame (classes/_instanceinfo.py) and fails +# when the caller's first local variable is an array. +FREQ = dict(samplingRate=FS, numSamples=N, signalType='power') + + +def cents(k, k0): + """Distance in cents from bin k0 to bin k.""" + return 1200 * np.log2(k / k0) + + +class TestSmoothFrequency(unittest.TestCase): + + def test_constant_magnitude_keeps_spectrum_and_phase(self): + # A flat magnitude is its own moving average on any grid, and the + # phase is not touched: the output equals the input, bin by bin. + rng = np.random.default_rng(1) + phase = rng.uniform(-np.pi, np.pi, K) + phase[0] = phase[-1] = 0 + spec = 0.3 * np.exp(1j * phase) + out = pytta.smooth_frequency(pytta.SignalObj(spec, 'freq', **FREQ), nthOct=3) + worst = np.max(np.abs(out.freqSignal[:, 0] - spec)) / 0.3 + print(f'\n flat spectrum, worst rel change = {worst:.2e}') + self.assertLess(worst, TOL32) + + def test_single_bin_spreads_into_a_flat_band_of_the_given_width(self): + # One bin at k0 = 100 is, on the cent grid, a triangle of base + # 1200 log2(101/99) ~ 34.6 cents and height 1. A box average of + # w cents turns it into a plateau of height area / w wherever the + # box holds the whole triangle: inside +-(h - 18) cents of k0, with + # h = (w - 1) / 2, and zero beyond +-(h + 18) cents. w is + # round(1200 / nthOct) made odd: 1199, 399 and 99 cents. + k0 = 100 + spec = np.zeros(K, dtype=complex) + spec[k0] = 1j # phase pi/2: the magnitude is read back below + area = 600 * np.log2((k0 + 1) / (k0 - 1)) + k = np.arange(1, K) + for nthOct, w in [(1, 1199), (3, 399), (12, 99)]: + out = pytta.smooth_frequency(pytta.SignalObj(spec, 'freq', **FREQ), + nthOct=nthOct) + y = np.abs(out.freqSignal[1:, 0]) + h = (w - 1) / 2 + inside = np.abs(cents(k, k0)) <= h - 40 + outside = np.abs(cents(k, k0)) >= h + 40 + plateau = y[inside] + flat = np.ptp(plateau) / plateau.mean() + err = plateau.mean() / (area / w) - 1 + leak = np.max(np.abs(y[outside])) / plateau.mean() + print(f'\n 1/{nthOct} oct: {inside.sum()} bins on plateau, ' + f'flatness {flat:.1e}, height err {err:+.3f}, leak {leak:.1e}') + self.assertLess(flat, 10 * TOL32) + self.assertLess(abs(err), 0.03) + self.assertLess(leak, TOL32) + self.assertEqual(out.freqSignal[0, 0], 0) + + def test_dc_and_nyquist_bins_are_kept(self): + # DC is copied through. At Nyquist the box shrinks to one grid + # point, so the magnitude there is also unchanged. + rng = np.random.default_rng(2) + spec = rng.standard_normal(K) + 1j * rng.standard_normal(K) + spec[0] = 2.0 + spec[-1] = -1.5 + out = pytta.smooth_frequency(pytta.SignalObj(spec, 'freq', **FREQ), nthOct=1) + self.assertAlmostEqual(out.freqSignal[0, 0], 2.0, delta=TOL32) + self.assertAlmostEqual(out.freqSignal[-1, 0], -1.5, delta=TOL32) + + def test_returns_new_object_and_smooths_channels_independently(self): + rng = np.random.default_rng(3) + sig = pytta.SignalObj(rng.standard_normal((N, 2)), 'time', FS) + before = sig.freqSignal.copy() + out = pytta.smooth_frequency(sig, nthOct=3) + self.assertTrue(np.array_equal(sig.freqSignal, before)) + self.assertIsNot(out, sig) + self.assertEqual(out.samplingRate, FS) + self.assertEqual(out.numSamples, N) + self.assertEqual(out.numChannels, 2) + single = pytta.smooth_frequency( + pytta.SignalObj(sig.timeSignal[:, 1], 'time', FS), nthOct=3) + worst = (np.max(np.abs(out.freqSignal[:, 1] - single.freqSignal[:, 0])) + / np.max(np.abs(single.freqSignal))) + self.assertLess(worst, TOL32) + + def test_rejects_invalid_arguments(self): + sig = pytta.SignalObj(np.ones(N), 'time', FS) + with self.assertRaises(ValueError): + pytta.smooth_frequency(sig, nthOct=0) + with self.assertRaises(ValueError): + pytta.smooth_frequency(sig, nthOct=-3) + with self.assertRaises(TypeError): + pytta.smooth_frequency(np.ones(N), nthOct=3) + + def test_matches_ita_toolbox(self): + # Sample by sample against ita_smooth_frequency on the stored + # signal, at 1, 1/3 and 1/12 octave. PyTTa and ITA scale the AC + # bins by different constants, measured here from the two input + # spectra (the signal has zero DC and Nyquist, see the generator). + self.assertTrue(os.path.exists(GOLDEN), + 'run golden/make_smooth_frequency_ita.m') + g = sio.loadmat(GOLDEN, squeeze_me=True, struct_as_record=False) + sig = pytta.SignalObj(g['x'], 'time', int(g['fs'])) + ratio = np.abs(g['freqDataIn'][1:-1]) / np.abs(sig.freqSignal[1:-1]) + self.assertLess(np.ptp(ratio) / ratio.mean(), TOL32, + 'input spectra are not one scale apart') + scale = ratio.mean() + for case in np.atleast_1d(g['cases']): + out = pytta.smooth_frequency(sig, nthOct=int(case.nthOct)) + peak = np.max(np.abs(case.freqData)) + diff = np.max(np.abs(out.freqSignal * scale - case.freqData)) / peak + print(f'\n 1/{int(case.nthOct):<2d} oct: rel peak diff vs ITA = ' + f'{diff:.2e} (tol {TOL32})') + self.assertLess(diff, TOL32) + + +if __name__ == '__main__': + unittest.main()