From 6e94eff7d1e6413cf403d9ccb6378eb2d452c047 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sat, 31 Oct 2026 20:02:21 +0100 Subject: [PATCH] [FEATURE] add minimum_phase for SignalObj Return the minimum phase version of a signal: the same magnitude spectrum, with the phase that puts the energy of each channel as early as possible. The phase comes from the real cepstrum folded onto positive times. The cepstrum is computed on the signal length, so the magnitude is kept exactly and the time aliasing of the cepstrum leaves a small phase error, which zero padding the input reduces. Tests pin a filter whose minimum phase version is known exactly, from its minimum and its maximum phase forms, the magnitude of a delayed noise burst, the early energy of the result and the metadata. Co-Authored-By: Claude Fable 5.1 --- pytta/__init__.py | 3 +- pytta/functions.py | 42 ++++++++++++++++++ tests/test_minimum_phase.py | 85 +++++++++++++++++++++++++++++++++++++ 3 files changed, 129 insertions(+), 1 deletion(-) create mode 100644 tests/test_minimum_phase.py 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()