diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..fdd8c9f 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, enhanced_spectrogram # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +146,7 @@ 'resample', 'corr_coef', 'peak_time', + 'enhanced_spectrogram', 'save', 'load', 'weighting', diff --git a/pytta/_plot.py b/pytta/_plot.py index c0a0970..1c9dc6a 100644 --- a/pytta/_plot.py +++ b/pytta/_plot.py @@ -26,6 +26,7 @@ import scipy.signal as ss import copy as cp from plotly import io +from pytta.utils.spectrogram import enhanced_map io.renderers.default = 'browser' @@ -804,7 +805,7 @@ def _curve_data_extractor_bars(analyses): def spectrogram(sigObjs, winType, winSize, overlap, xLabel, yLabel, xLim, yLim, - title, decimalSep): + title, decimalSep, enhanced=None): """ Plots a signal spectrogram in frequency domain. @@ -846,6 +847,10 @@ def spectrogram(sigObjs, winType, winSize, may be dot or comma. >>> decimalSep = ',' # in Brazil + * enhanced (None), (str): + 'readable' or 'sharp' plots the enhanced spectrogram + (pytta.enhanced_spectrogram) on a logarithmic frequency axis; + winType, winSize and overlap are then not used. Return: -------- @@ -867,13 +872,23 @@ def spectrogram(sigObjs, winType, winSize, projection='rectilinear') ax.set_snap(False) - _spectrogram, _specTime, _specFreq = \ - _calc_spectrogram(data['timeSignal'], data['timeVector'], - data['samplingRate'], overlap, winType, - winSize, data['dBRef']) - - pcmesh = ax.pcolormesh(_specTime, _specFreq, _spectrogram, - cmap=plt.jet(), vmin=-120) + if enhanced: + _specTime, _specFreq, meanSquare = \ + enhanced_map(data['timeSignal'], data['samplingRate'], + enhanced) + _spectrogram = 10 * np.log10(meanSquare / data['dBRef'] ** 2 + + 1e-12) + pcmesh = ax.pcolormesh(_specTime, _specFreq, _spectrogram, + cmap=plt.jet(), shading='nearest', + vmin=_spectrogram.max() - 80) + ax.set_yscale('log') + else: + _spectrogram, _specTime, _specFreq = \ + _calc_spectrogram(data['timeSignal'], data['timeVector'], + data['samplingRate'], overlap, winType, + winSize, data['dBRef']) + pcmesh = ax.pcolormesh(_specTime, _specFreq, _spectrogram, + cmap=plt.jet(), vmin=-120) if xLim is None: xLim = (data['timeVector'][0], data['timeVector'][-1]) @@ -883,8 +898,15 @@ def spectrogram(sigObjs, winType, winSize, yLim = (data['minFreq'], data['maxFreq']) ax.set_ylim(yLim) - ax.yaxis.set_major_locator(ticker.MaxNLocator(min_n_ticks=8)) - ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useOffset=True)) + if enhanced: + octaves = [f for f in (31.5, 63, 125, 250, 500, 1000, 2000, 4000, + 8000, 16000) if yLim[0] <= f <= yLim[1]] + ax.set_yticks(octaves) + ax.set_yticklabels([str(f) for f in octaves]) + ax.yaxis.set_minor_formatter(ticker.NullFormatter()) + else: + ax.yaxis.set_major_locator(ticker.MaxNLocator(min_n_ticks=8)) + ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useOffset=True)) ax.xaxis.set_major_locator(ticker.MaxNLocator(min_n_ticks=10)) ax.xaxis.set_major_formatter(ticker.ScalarFormatter(useOffset=True)) for item in (ax.get_xticklabels() + ax.get_yticklabels()): diff --git a/pytta/classes/signal.py b/pytta/classes/signal.py index 79151c3..d358fbd 100644 --- a/pytta/classes/signal.py +++ b/pytta/classes/signal.py @@ -738,7 +738,8 @@ def plot_freq(self, smooth:bool=False, xLabel:str=None, yLabel:str=None, def plot_spectrogram(self, winType:str='hann', winSize:int=1024, overlap:float=0.5, xLabel:str=None, yLabel:str=None, yLim:list=None, xLim:list=None, title:str=None, - decimalSep:str=','): + decimalSep:str=',', + enhanced:str=None): """Plots a signal spectrogram. xLabel, yLabel, and title are saved for the next plots when provided. @@ -778,6 +779,10 @@ def plot_spectrogram(self, winType:str='hann', winSize:int=1024, may be dot or comma. >>> decimalSep = ',' # in Brazil + * enhanced (None), (str): + 'readable' or 'sharp' plots the enhanced spectrogram + (pytta.enhanced_spectrogram) on a logarithmic frequency + axis; winType, winSize and overlap are then not used. Return: -------- @@ -807,7 +812,7 @@ def plot_spectrogram(self, winType:str='hann', winSize:int=1024, figs = plot.spectrogram((self,), winType, winSize, overlap, xLabel, yLabel, xLim, yLim, - title, decimalSep) + title, decimalSep, enhanced) return figs def calib_voltage(self, chIndex, refSignalObj, refVrms=1, refFreq=1000): diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..e214f03 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.enhanced_spectrogram(signalObj) >>> pytta.plot_time(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_time_dB(signalObj1, signalObj2, ..., signalObjN) >>> pytta.plot_freq(signalObj1, signalObj2, ..., signalObjN) @@ -49,6 +50,7 @@ import copy as cp from warnings import warn from pytta import _plot as plot +from pytta.utils.spectrogram import enhanced_map # For backwards compatibility purposes. Planned to get out of here from pytta.utils.maths import fft_degree as new_fft_degree @@ -274,6 +276,62 @@ def peak_time(signal): else: return peaks_time[0] +def enhanced_spectrogram(signalObj, smoothing='readable', n_frames=2000, + f_min=20, bins_per_octave=192): + """ + Enhanced spectrogram: a spectrogram with no window to choose. + + Nine Blackman-Harris windows from 8 to 512 ms (geometric spacing) are + reassigned (Auger and Flandrin, 1995), the energy of each frame is held + over the frame step of its window, each map is smoothed, and the maps are + combined by geometric mean, so that only the energy that all the windows + place at the same point of the time-frequency plane remains (the idea of + Cheung and Lim, 1991, applied to reassigned maps). + + Parameters (default), (type): + ----------------------------- + * signalObj (), (SignalObj): + signal to analyse. It is not modified; + * smoothing ('readable'), (str): + 'readable': Gaussian smoothing of 4 ms in time and 1.45 % of the + frequency, for continuous lines; 'sharp': 1 ms and 1 Hz, for the + thinnest lines; + * n_frames (2000), (int): + largest number of time columns. The time step is 1 ms or more; + * f_min (20), (float): + lowest frequency of the output, in Hz; + * bins_per_octave (192), (int): + frequency resolution of the logarithmic output grid. + + Return: + ------- + * t: time of each column, in s, shape (M,); + * f: frequency of each row, in Hz, shape (F,); + * L: level of each cell in dB re the channel's dBRef (dB SPL for a + channel in Pa), shape (F, M, numChannels). A cell holds the mean + square value of its share of the signal, so the cells of a column + add up to the level of that instant; a steady tone spread by the + smoothing over several cells reads lower in each of them (about + 10 dB with 'readable'). Silence sits at -120 dB. + + Author: Sergio Aguirre, 2026. Released in PyTTa under the MIT license. + The same method is in the SQAT toolbox (GPL-3.0) as + gui/SQAT_GUI_enhanced_stft.m (https://github.com/ggrecow/SQAT). + + >>> t, f, L = pytta.enhanced_spectrogram(signalObj, smoothing='sharp') + """ + if not isinstance(signalObj, SignalObj): + raise TypeError('signalObj must be a SignalObj.') + L = [] + for chIndex, chNum in enumerate(signalObj.channels.mapping): + t, f, meanSquare = enhanced_map(signalObj.timeSignal[:, chIndex], + signalObj.samplingRate, smoothing, + n_frames, f_min, bins_per_octave) + dBRef = signalObj.channels[chNum].dBRef + L.append(10 * np.log10(meanSquare / dBRef ** 2 + 1e-12)) + return t, f, np.stack(L, axis=2) + + def fft_degree(*args,**kwargs): """ DEPRECATED @@ -533,7 +591,8 @@ def plot_bars(*analyses, xLabel:str=None, yLabel:str=None, def plot_spectrogram(*sigObjs, winType:str='hann', winSize:int=1024, overlap:float=0.5, xLabel:str=None, yLabel:str=None, yLim:list=None, xLim:list=None, title:str=None, - decimalSep:str=','): + decimalSep:str=',', + enhanced:str=None): """ Plots provided SignalObjs spectrogram. @@ -575,6 +634,10 @@ def plot_spectrogram(*sigObjs, winType:str='hann', winSize:int=1024, may be dot or comma. >>> decimalSep = ',' # in Brazil + * enhanced (None), (str): + 'readable' or 'sharp' plots the enhanced spectrogram + (pytta.enhanced_spectrogram) on a logarithmic frequency axis; + winType, winSize and overlap are then not used. Return: -------- @@ -586,7 +649,7 @@ def plot_spectrogram(*sigObjs, winType:str='hann', winSize:int=1024, if len(realSigObjs) > 0: figs = plot.spectrogram(realSigObjs, winType, winSize, overlap, xLabel, yLabel, xLim, yLim, - title, decimalSep) + title, decimalSep, enhanced) return figs else: return diff --git a/pytta/utils/__init__.py b/pytta/utils/__init__.py index ff5c7c8..f7a4a87 100755 --- a/pytta/utils/__init__.py +++ b/pytta/utils/__init__.py @@ -15,6 +15,7 @@ * colore * freq * maths + * spectrogram Created on Tue May 5 00:31:42 2020 diff --git a/pytta/utils/spectrogram.py b/pytta/utils/spectrogram.py new file mode 100644 index 0000000..7c6f7a7 --- /dev/null +++ b/pytta/utils/spectrogram.py @@ -0,0 +1,144 @@ +""" +Enhanced spectrogram on NumPy arrays. + +Author: Sergio Aguirre, 2026. Released in PyTTa under the MIT license. +The same method is in the SQAT toolbox (GPL-3.0) as +gui/SQAT_GUI_enhanced_stft.m (https://github.com/ggrecow/SQAT). +The SignalObj interface is pytta.enhanced_spectrogram. +""" +import numpy as np +from scipy.ndimage import convolve1d + +__all__ = ['enhanced_map'] + + +def enhanced_map(x, fs, smoothing='readable', n_frames=2000, f_min=20, + bins_per_octave=192): + """ + Enhanced spectrogram of one channel, as a map of mean square values. + + Nine Blackman-Harris windows from 8 to 512 ms (geometric spacing) are + reassigned (Auger and Flandrin, 1995), the energy of each frame is held + over the frame step of its window, each map is smoothed, and the maps are + combined by geometric mean, so that only the energy that all the windows + place at the same point of the time-frequency plane remains (the idea of + Cheung and Lim, 1991, applied to reassigned maps). + + Returns t (M,) in s, f (F,) in Hz on a logarithmic grid, and the map + (F, M) of mean square values in the unit of x squared: the cells of a + column add up to the mean square value of that instant. See + pytta.enhanced_spectrogram for the parameters. + """ + if f_min >= fs / 2: + raise ValueError(f'f_min ({f_min} Hz) must be below fs/2 ' + f'({fs / 2} Hz).') + if smoothing not in ('readable', 'sharp'): + raise ValueError("smoothing must be 'readable' or 'sharp'.") + x = np.asarray(x, dtype=float).ravel() + numSamples = x.size + durations = 8e-3 * 64 ** (np.arange(9) / 8) + windows = 2 * _round_half_away(durations * fs / 2).astype(int) + hop = max(1, int(_round_half_away(0.001 * fs)), + int(np.ceil(numSamples / n_frames))) + steps = np.maximum(1, _round_half_away(windows / 8).astype(int)) + numColumns = len(range(0, numSamples, hop)) + numRows = int(np.ceil(np.log2((fs / 2) / f_min) * bins_per_octave)) + 1 + f = f_min * 2 ** (np.arange(numRows) / bins_per_octave) + cellWidth = f * (2 ** (1 / bins_per_octave) - 1) + if smoothing == 'readable': + sigmaTime = 0.004 * fs / hop + sigmaFreq = 2 * (2 ** (1 / 96) - 1) * f / cellWidth + else: + sigmaTime = 0.001 * fs / hop + sigmaFreq = 1 / cellWidth + timeKernel = _gaussian_kernel(max(sigmaTime, 0.3)) + sigmaFreq = np.maximum(sigmaFreq, 0.3) + logSum = np.zeros((numColumns, numRows)) + for N, step in zip(windows, steps): + energy = np.concatenate(([0.], np.cumsum(np.concatenate( + (np.zeros(N), x ** 2, np.zeros(N)))))) + frameMax = np.max(energy[N:] - energy[:-N]) + R = _reassigned_map(x, fs, N, step, hop, numColumns, numRows, + f_min, bins_per_octave, frameMax) / hop + numHold = 2 * int(np.floor(step / hop / 2)) + 1 + if numHold > 1: + R = convolve1d(R, np.ones(numHold) / numHold, axis=0, + mode='constant') + R = convolve1d(R, timeKernel, axis=0, mode='constant') + R = _smooth_frequency_rows(R, sigmaFreq) + with np.errstate(divide='ignore'): + logSum += np.log(R + 1e-6 * R.max()) + t = np.arange(numColumns) * hop / fs + return t, f, np.exp(logSum / windows.size).T + + +def _round_half_away(x): + """Round halves away from zero, as MATLAB does.""" + x = np.asarray(x, dtype=float) + return np.sign(x) * np.floor(np.abs(x) + 0.5) + + +def _gaussian_kernel(sigma): + half = int(np.ceil(3 * sigma)) + kernel = np.exp(-np.arange(-half, half + 1) ** 2 / (2 * sigma ** 2)) + return kernel / kernel.sum() + + +def _smooth_frequency_rows(R, sigma): + """Gaussian smoothing along the frequency axis (columns of R) with a + width that changes with frequency: one kernel per geometric level.""" + level = _round_half_away(np.log(sigma) / np.log(1.25)).astype(int) + out = R.copy() + numFreqs = R.shape[1] + for lev in np.unique(level): + cols = np.flatnonzero(level == lev) + kernel = _gaussian_kernel(1.25 ** lev) + half = kernel.size // 2 + c1 = max(0, cols[0] - half) + c2 = min(numFreqs, cols[-1] + half + 1) + smoothed = convolve1d(R[:, c1:c2], kernel, axis=1, mode='constant') + out[:, cols] = smoothed[:, cols - c1] + return out + + +def _reassigned_map(x, fs, N, step, hop, numColumns, numRows, f_min, + bins_per_octave, frameMax): + """Reassigned spectrogram of one unit-energy Blackman-Harris window of N + samples and frame step `step`, accumulated on the output grid.""" + half = N // 2 + padded = np.concatenate((np.zeros(half), x, np.zeros(half + N))) + n = np.arange(-half, half) + phase = 2 * np.pi * (n + half) / N + a = (0.35875, 0.48829, 0.14128, 0.01168) + w = a[0] - a[1] * np.cos(phase) + a[2] * np.cos(2 * phase) \ + - a[3] * np.cos(3 * phase) + wd = a[1] * (2 * np.pi / N) * np.sin(phase) \ + - a[2] * (4 * np.pi / N) * np.sin(2 * phase) \ + + a[3] * (6 * np.pi / N) * np.sin(3 * phase) + norm = 1 / np.linalg.norm(w) + w, wd = w * norm, wd * norm + wt = n * w + freqBins = np.arange(half + 1) * fs / N + starts = np.arange(0, x.size, step) + R = np.zeros(numColumns * numRows) + chunk = max(1, int(2e6 // N)) + for b in range(0, starts.size, chunk): + s = starts[b:b + chunk] + frames = padded[s[:, None] + np.arange(N)[None, :]] + Xh = np.fft.rfft(frames * w, axis=1) + Xt = np.fft.rfft(frames * wt, axis=1) + Xd = np.fft.rfft(frames * wd, axis=1) + power = np.abs(Xh) ** 2 + r = np.conj(Xh) / np.maximum(power, np.finfo(float).eps) + tHat = s[:, None] + np.real(Xt * r) + fHat = freqBins[None, :] - fs / (2 * np.pi) * np.imag(Xd * r) + col = _round_half_away(tHat / hop) + with np.errstate(divide='ignore'): + row = _round_half_away(bins_per_octave * np.log2( + np.maximum(fHat, np.finfo(float).eps) / f_min)) + ok = (power > 1e-10 * frameMax) & (col >= 0) & (col < numColumns) \ + & (row >= 0) & (row < numRows) & (fHat > 0) + cells = (col[ok] + row[ok] * numColumns).astype(int) + R += np.bincount(cells, weights=power[ok] * (2 * step / N), + minlength=R.size) + return R.reshape(numRows, numColumns).T diff --git a/tests/golden/enhanced_spectrogram_sqat.mat b/tests/golden/enhanced_spectrogram_sqat.mat new file mode 100644 index 0000000..deaa08b Binary files /dev/null and b/tests/golden/enhanced_spectrogram_sqat.mat differ diff --git a/tests/golden/make_enhanced_spectrogram_sqat.m b/tests/golden/make_enhanced_spectrogram_sqat.m new file mode 100644 index 0000000..c794135 --- /dev/null +++ b/tests/golden/make_enhanced_spectrogram_sqat.m @@ -0,0 +1,45 @@ +function make_enhanced_spectrogram_sqat(sqatRoot) +% MAKE_ENHANCED_SPECTROGRAM_SQAT Golden data for tests/test_enhanced_spectrogram.py. +% +% make_enhanced_spectrogram_sqat(sqatRoot) runs SQAT_GUI_enhanced_stft +% from SQAT (https://github.com/ggrecow/SQAT, gui/) on a fixed test signal +% with both smoothings and the default grid, and writes +% enhanced_spectrogram_sqat.mat next to this file. +% +% The signal (0.6 s at 16 kHz, in Pa) holds a 1 kHz tone at 80 dB SPL, a +% 250 Hz tone modulated at 4 Hz, a click at 0.3 s and white noise at +% 40 dB SPL. It is rounded to single precision, as PyTTa stores it, so +% both sides read the same numbers. The maps are stored on every 16th +% frequency and every 10th column, enough to pin the method without +% storing the whole grid; t and f are stored whole. +% +% Usage: matlab -batch "make_enhanced_spectrogram_sqat('/path/to/SQAT')" + +addpath(fullfile(sqatRoot, 'gui')); + +fs = 16000; +n = round(0.6 * fs); +tt = (0:n-1)' / fs; +rng(0, 'twister'); +x = 0.2 * sqrt(2) * sin(2*pi*1000*tt) ... + + 0.05 * (1 + 0.8*sin(2*pi*4*tt)) .* sin(2*pi*250*tt) ... + + 2e-3 * randn(n, 1); +x(round(0.3 * fs) + 1) = x(round(0.3 * fs) + 1) + 1; +x = double(single(x)); + +modes = {'readable', 'sharp'}; +L = cell(1, 2); +for k = 1:2 + [t, f, Lk] = SQAT_GUI_enhanced_stft(x, fs, modes{k}); + rows = 1:16:numel(f); + cols = 1:10:numel(t); + L{k} = Lk(rows, cols); +end +L_readable = L{1}; +L_sharp = L{2}; + +out = fullfile(fileparts(mfilename('fullpath')), 'enhanced_spectrogram_sqat.mat'); +x = single(x); +save(out, 'x', 'fs', 't', 'f', 'rows', 'cols', 'L_readable', 'L_sharp', '-v7'); +fprintf('wrote %s (%d x %d per map)\n', out, numel(rows), numel(cols)); +end diff --git a/tests/test_enhanced_spectrogram.py b/tests/test_enhanced_spectrogram.py new file mode 100644 index 0000000..2060d90 --- /dev/null +++ b/tests/test_enhanced_spectrogram.py @@ -0,0 +1,133 @@ +""" +Tests for pytta.enhanced_spectrogram. + +The enhanced spectrogram is the geometric mean of nine reassigned +spectrograms with Blackman-Harris windows of 8 to 512 ms, smoothed +('readable' or 'sharp'), on a logarithmic frequency grid, in dB re the +channel's dBRef (20 uPa for a channel in Pa). Each test pins one rule. The +first compares the whole method with SQAT_GUI_enhanced_stft of SQAT on data +generated by golden/make_enhanced_spectrogram_sqat.m. +""" +import os +import unittest + +os.environ.setdefault("MPLBACKEND", "Agg") + +import numpy as np +import scipy.io as sio + +import pytta + +GOLDEN = os.path.join(os.path.dirname(__file__), 'golden', + 'enhanced_spectrogram_sqat.mat') +FS = 16000 + + +def in_pascal(signalObj): + for ch in signalObj.channels.mapping: + signalObj.channels[ch].unit = 'Pa' + return signalObj + + +def tone(freq, rms, seconds): + tt = np.arange(int(seconds * FS)) / FS + return rms * np.sqrt(2) * np.sin(2 * np.pi * freq * tt) + + +class TestEnhancedSpectrogram(unittest.TestCase): + + def test_matches_sqat(self): + # Against SQAT on the stored signal, both smoothings, default grid: + # t and f whole, the maps on every 16th row and 10th column. Both + # sides compute in double from the same float32 samples, so the + # tolerance is far below any visible change and far above the + # differences of FFT summation order (measured: ~4e-10 dB). + self.assertTrue(os.path.exists(GOLDEN), + 'run golden/make_enhanced_spectrogram_sqat.m') + g = sio.loadmat(GOLDEN, squeeze_me=True) + sig = in_pascal(pytta.SignalObj(g['x'], 'time', int(g['fs']))) + rows = np.atleast_1d(g['rows']) - 1 + cols = np.atleast_1d(g['cols']) - 1 + for smoothing in ('readable', 'sharp'): + t, f, L = pytta.enhanced_spectrogram(sig, smoothing=smoothing) + self.assertLess(np.max(np.abs(t - g['t'])), 1e-12) + self.assertLess(np.max(np.abs(f / g['f'] - 1)), 1e-12) + ref = g['L_' + smoothing] + got = L[np.ix_(rows, cols)][:, :, 0] + diff = np.max(np.abs(got - ref)) + print(f'\n {smoothing:8s} max |dL| vs SQAT = {diff:.2e} dB ' + f'on {ref.size} cells (tol 1e-6)') + self.assertLess(diff, 1e-6) + + def test_steady_tone_column_adds_up_to_its_level(self): + # The cells of a column hold mean square pressure, so a column in the + # middle of a 94 dB SPL tone (1 Pa rms) adds up to 94 dB. Predicted + # before running: every window reassigns a steady tone to its true + # frequency, so the nine maps agree and the geometric mean keeps the + # energy; within 0.2 dB with either smoothing (measured: 93.94 dB + # readable, 94.06 dB sharp). + sig = in_pascal(pytta.SignalObj(tone(1000, 1.0, 1.0), 'time', FS)) + for smoothing in ('readable', 'sharp'): + t, f, L = pytta.enhanced_spectrogram(sig, smoothing=smoothing) + column = L[:, len(t) // 2, 0] + level = 10 * np.log10(np.sum(10 ** (column / 10))) + print(f'\n {smoothing:8s} column level = {level:.2f} dB SPL') + self.assertLess(abs(level - 94), 0.2) + + def test_tone_and_click_land_where_they_are(self): + # A 1 kHz tone peaks on the frequency row nearest 1 kHz, and a click + # at 0.5 s peaks within one column (1 ms) of 0.5 s on a row at 4 kHz. + sig = in_pascal(pytta.SignalObj(tone(1000, 0.1, 1.0), 'time', FS)) + t, f, L = pytta.enhanced_spectrogram(sig, smoothing='sharp') + peak = f[np.argmax(L[:, len(t) // 2, 0])] + self.assertEqual(peak, f[np.argmin(np.abs(np.log2(f / 1000)))]) + click = np.zeros(FS) + click[FS // 2] = 1 + sig = in_pascal(pytta.SignalObj(click, 'time', FS)) + t, f, L = pytta.enhanced_spectrogram(sig, smoothing='sharp') + row = np.argmin(np.abs(np.log2(f / 4000))) + self.assertLessEqual(abs(t[np.argmax(L[row, :, 0])] - 0.5), 1e-3) + + def test_silence_sits_at_the_floor(self): + # log of an empty map: every cell at -120 dB, no NaN or -inf. + sig = pytta.SignalObj(np.zeros(FS // 4), 'time', FS) + t, f, L = pytta.enhanced_spectrogram(sig) + self.assertTrue(np.all(L == -120)) + + def test_channels_are_independent_and_input_is_kept(self): + rng = np.random.default_rng(0) + data = np.column_stack([tone(500, 0.1, 0.5), + 0.05 * rng.standard_normal(FS // 2)]) + sig = in_pascal(pytta.SignalObj(data, 'time', FS)) + before = sig.timeSignal.copy() + t, f, L = pytta.enhanced_spectrogram(sig) + self.assertTrue(np.array_equal(sig.timeSignal, before)) + self.assertEqual(L.shape, (f.size, t.size, 2)) + single = in_pascal(pytta.SignalObj(sig.timeSignal[:, 1], 'time', FS)) + _, _, L1 = pytta.enhanced_spectrogram(single) + self.assertTrue(np.array_equal(L[:, :, 1], L1[:, :, 0])) + + def test_plot_shows_the_enhanced_map(self): + # plot_spectrogram(enhanced=...) draws exactly the map of + # enhanced_spectrogram, on a logarithmic frequency axis. + sig = in_pascal(pytta.SignalObj(tone(1000, 0.1, 0.5), 'time', FS)) + figs = pytta.plot_spectrogram(sig, enhanced='sharp', decimalSep='.') + ax = figs[0].axes[0] + _, _, L = pytta.enhanced_spectrogram(sig, smoothing='sharp') + drawn = np.asarray(ax.collections[0].get_array()).reshape(L.shape[:2]) + self.assertTrue(np.array_equal(drawn, L[:, :, 0])) + self.assertEqual(ax.get_yscale(), 'log') + figs[0].clf() + + def test_rejects_invalid_arguments(self): + sig = pytta.SignalObj(np.zeros(FS // 4), 'time', FS) + with self.assertRaises(ValueError): + pytta.enhanced_spectrogram(sig, f_min=FS / 2) + with self.assertRaises(ValueError): + pytta.enhanced_spectrogram(sig, smoothing='blurry') + with self.assertRaises(TypeError): + pytta.enhanced_spectrogram(np.zeros(FS), smoothing='sharp') + + +if __name__ == '__main__': + unittest.main()