diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..0a81d84 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, minimum_phase # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +146,7 @@ 'resample', 'corr_coef', 'peak_time', + 'minimum_phase', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..75a625c 100644 --- a/pytta/functions.py +++ b/pytta/functions.py @@ -17,6 +17,7 @@ >>> pytta.corr_coef(signalObj1, signalObj2) >>> pytta.resample(signalObj, newSamplingRate) >>> pytta.peak_time(signalObj1, signalObj2, ..., signalObjN) + >>> pytta.minimum_phase(signalObj) >>> pytta.plot_time(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_time_dB(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_freq(signalObj1, signalObj2, ..., signalObjN) @@ -274,6 +275,47 @@ def peak_time(signal): else: return peaks_time[0] +def minimum_phase(signalObj): + """ + Minimum phase version of a SignalObj: the same magnitude spectrum, with + the phase that puts the energy of each channel as early as possible. + Returns a new SignalObj of the same length. + + The phase comes from the real cepstrum (Oppenheim and Schafer, Discrete + Time Signal Processing, ch. 13): the logarithm of the magnitude is taken + to the time domain, folded onto positive times and brought back. The + cepstrum is computed on the signal length, so the magnitude is kept + exactly; its time aliasing leaves a small error in the phase, which + zero padding the signal before the call reduces. + + >>> mp = pytta.minimum_phase(ir) + """ + if not isinstance(signalObj, SignalObj): + raise TypeError('signalObj must be a SignalObj.') + N = signalObj.numSamples + spectrum = np.fft.rfft(np.asarray(signalObj.timeSignal, dtype=float), + axis=0) + magnitude = np.abs(spectrum) + # ponytail: spectral zeros floored at -200 dB re the peak, so log is finite + magnitude = np.maximum(magnitude, 1e-10 * magnitude.max(axis=0)) + cepstrum = np.fft.irfft(np.log(magnitude), n=N, axis=0) + fold = np.zeros(N) + fold[0] = 1 + fold[1:(N + 1) // 2] = 2 + if N % 2 == 0: + fold[N // 2] = 1 + minPhase = np.fft.irfft(np.exp(np.fft.rfft(cepstrum * fold[:, None], + axis=0)), n=N, axis=0) + newSignal = SignalObj(minPhase, 'time', + samplingRate=signalObj.samplingRate, + 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 diff --git a/tests/test_minimum_phase.py b/tests/test_minimum_phase.py new file mode 100644 index 0000000..61551ee --- /dev/null +++ b/tests/test_minimum_phase.py @@ -0,0 +1,85 @@ +""" +Tests for pytta.minimum_phase. + +The minimum phase version of a signal has the same magnitude spectrum and +the phase that puts its energy as early as possible. Each test pins one +rule. The first uses a short filter whose minimum phase version is known +exactly: zeros inside the unit circle give a minimum phase filter, and its +time reversal has the same magnitude and maximum phase. +""" +import os +import unittest + +os.environ.setdefault('MPLBACKEND', 'Agg') + +import numpy as np + +import pytta + +FS = 48000 +N = 1024 +TOL32 = 1e-6 # rounding of the float32 SignalObj storage + + +class TestMinimumPhase(unittest.TestCase): + + def test_known_filter_comes_back_in_minimum_phase(self): + b = np.real(np.poly([0.9, -0.7, 0.5 * np.exp(0.4j), + 0.5 * np.exp(-0.4j)])) + data = np.zeros((N, 2)) + data[:b.size, 0] = b # minimum phase already + data[:b.size, 1] = b[::-1] # same magnitude, maximum phase + sig = pytta.SignalObj(data, 'time', FS) + out = pytta.minimum_phase(sig).timeSignal + for ch in range(2): + err = np.max(np.abs(out[:, ch] - data[:, 0])) + print(f'\n channel {ch}: max |error| = {err:.1e}') + self.assertLess(err, TOL32) + + def test_magnitude_is_kept(self): + # A delayed, decaying noise burst: far from minimum phase. + rng = np.random.default_rng(0) + t = np.arange(N) + x = rng.standard_normal(N) * np.exp(-(t - 100) / 80) * (t >= 100) + sig = pytta.SignalObj(x, 'time', FS) + out = pytta.minimum_phase(sig) + X = np.abs(np.fft.rfft(sig.timeSignal[:, 0])) + Y = np.abs(np.fft.rfft(out.timeSignal[:, 0])) + loud = X > 1e-3 * X.max() + err = np.max(np.abs(Y[loud] / X[loud] - 1)) + print(f'\n max relative magnitude change = {err:.1e}') + self.assertLess(err, 1e-5) + + def test_energy_comes_as_early_as_possible(self): + # The cumulative energy of the minimum phase version is at least + # that of the input at every sample. The cepstrum is computed on the + # signal length, so its time aliasing allows a small deficit: + # 2e-6 of the total energy on this signal, 3e-4 on the burst of + # test_magnitude_is_kept. + rng = np.random.default_rng(1) + t = np.arange(N) + x = rng.standard_normal(N) * np.exp(-(t - 300) / 60) * (t >= 300) + sig = pytta.SignalObj(x, 'time', FS) + y = pytta.minimum_phase(sig).timeSignal[:, 0].astype(float) + x = sig.timeSignal[:, 0].astype(float) + deficit = np.max(np.cumsum(x ** 2) - np.cumsum(y ** 2)) / np.sum(x ** 2) + print(f'\n largest energy deficit = {deficit:.1e} of the total') + self.assertLess(deficit, 1e-3) + self.assertLess(np.argmax(np.abs(y)), 50) + + def test_returns_new_object_and_keeps_metadata(self): + rng = np.random.default_rng(2) + sig = pytta.SignalObj(rng.standard_normal((N, 2)), 'time', FS, + signalType='energy') + before = sig.timeSignal.copy() + out = pytta.minimum_phase(sig) + self.assertTrue(np.array_equal(sig.timeSignal, before)) + self.assertIsNot(out, sig) + self.assertEqual((out.samplingRate, out.numSamples, out.numChannels, + out.signalType), (FS, N, 2, 'energy')) + with self.assertRaises(TypeError): + pytta.minimum_phase(before) + + +if __name__ == '__main__': + unittest.main()