From cc291091d17d4544d4ecdcea68619a79b739c585 Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 00:56:00 +0100 Subject: [PATCH] [FEATURE] add psd, interpolate_spectrum and extrapolate_low Adds to pytta.functions three spectrum utilities that the ITA-Toolbox has: psd (ita_psd, the power or cross power spectral density by Welch averaging with the window, overlap and convention of MATLAB cpsd), interpolate_spectrum (ita_interpolate_spk, the DFT interpolated onto the bins of another length, by spline or linearly, on real and imaginary parts or on magnitude and phase) and extrapolate_low (the lower extrapolation of ita_spk_extrapolate, mean magnitude with an optional slope and the linear phase of the mean group delay). Tests check the densities against MATLAB cpsd on golden data from tests/golden/make_psd_matlab.m (3.4e-8 relative) and Parseval, that the interpolated spectrum of a one-pole filter follows its true response (1.4e-6), and that the extrapolation keeps a pure delay and gives the slope asked for. Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/functions.rst | 4 +- pytta/__init__.py | 6 +- pytta/functions.py | 107 +++++++++++++++++++++++++++++++++ tests/golden/make_psd_matlab.m | 12 ++++ tests/golden/psd_matlab.mat | Bin 0 -> 6258 bytes tests/test_spectrum.py | 80 ++++++++++++++++++++++++ 6 files changed, 207 insertions(+), 2 deletions(-) create mode 100644 tests/golden/make_psd_matlab.m create mode 100644 tests/golden/psd_matlab.mat create mode 100644 tests/test_spectrum.py diff --git a/docs/pyttadoc/functions.rst b/docs/pyttadoc/functions.rst index add57c0..e6cf5d1 100644 --- a/docs/pyttadoc/functions.rst +++ b/docs/pyttadoc/functions.rst @@ -26,4 +26,6 @@ Functions .. autofunction:: pytta.plot_bars .. autofunction:: pytta.plot_spectrogram .. autofunction:: pytta.plot_waterfall - +.. autofunction:: pytta.psd +.. autofunction:: pytta.interpolate_spectrum +.. autofunction:: pytta.extrapolate_low diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..3114eed 100644 --- a/pytta/__init__.py +++ b/pytta/__init__.py @@ -113,7 +113,8 @@ 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,\ + psd, interpolate_spectrum, extrapolate_low # Must go on v0.1.0 (DEPRECATED) from .functions import fft_degree @@ -146,6 +147,9 @@ 'resample', 'corr_coef', 'peak_time', + 'psd', + 'interpolate_spectrum', + 'extrapolate_low', 'save', 'load', 'weighting', diff --git a/pytta/functions.py b/pytta/functions.py index 9b39239..1c656e4 100644 --- a/pytta/functions.py +++ b/pytta/functions.py @@ -23,6 +23,9 @@ >>> pytta.plot_bars(signalObj1, signalObj2, ..., signalObjN) >>> pytta.save(fileName, obj1, ..., objN) >>> pytta.load(fileName) + >>> pytta.psd(signalObj) + >>> pytta.interpolate_spectrum(signalObj, nSamples) + >>> pytta.extrapolate_low(signalObj, freqRange) For further information, check the function specific documentation. @@ -274,6 +277,110 @@ def peak_time(signal): else: return peaks_time[0] +def psd(signal, signal2=None, block_size=None, fft_size=None): + """ + Power spectral density of each channel, or the cross power spectral + density with the same channel of signal2, as ita_psd: Welch averaging + with a symmetric Hann window of block_size samples (the whole signal by + default) and half-block overlap, one-sided, in units^2/Hz, with the + convention of MATLAB cpsd, X conj(Y) (the conjugate of + scipy.signal.csd). + + Returns + ------- + (numpy.ndarray, numpy.ndarray) + Frequencies (Hz) and the density, (bins, channels). + + >>> freq, P = pytta.psd(noise, block_size=1024) + """ + x = np.asarray(signal.timeSignal, dtype=float) + y = x if signal2 is None else np.asarray(signal2.timeSignal, dtype=float) + if y.shape != x.shape: + raise ValueError('both signals need the same length and channels.') + block = x.shape[0] if block_size is None else int(round(block_size)) + block += block % 2 + if block > x.shape[0]: + x = np.vstack([x, np.zeros((block - x.shape[0], x.shape[1]))]) + y = np.vstack([y, np.zeros((block - y.shape[0], y.shape[1]))]) + nfft = block if fft_size is None else int(fft_size) + window = np.hanning(block + 2)[1:-1] # MATLAB hanning(block) + f, P = ss.csd(x, y, fs=signal.samplingRate, window=window, noverlap=block // 2, + nfft=nfft, axis=0, detrend=False, scaling='density') + return f, np.conj(P) if signal2 is not None else P.real # MATLAB cpsd: X conj(Y) + + +def interpolate_spectrum(signal, n_samples, method='spline', abs_phase=False): + """ + signal with its spectrum interpolated onto the bins of n_samples, as + ita_interpolate_spk: the DFT is interpolated between DC and the highest + frequency (spline, as MATLAB interp1, or linear), on its real and + imaginary parts or (abs_phase=True) on magnitude and unwrapped phase. + Meant for energy signals such as impulse responses. + + >>> longer = pytta.interpolate_spectrum(ir, 2 ** 16) + """ + from scipy.interpolate import CubicSpline + x = np.asarray(signal.timeSignal, dtype=float) + X = np.fft.rfft(x, axis=0) + old = np.linspace(0, 1, X.shape[0]) + new = np.linspace(0, 1, n_samples // 2 + 1) + if method == 'spline': + interp = lambda v: CubicSpline(old, v, axis=0)(new) + elif method == 'linear': + interp = lambda v: np.column_stack([np.interp(new, old, v[:, c]) for c in range(v.shape[1])]) + else: + raise ValueError("method must be 'spline' or 'linear'.") + if abs_phase: + Y = interp(np.abs(X)) * np.exp(1j * interp(np.unwrap(np.angle(X), axis=0))) + else: + Y = interp(X.real) + 1j * interp(X.imag) + newSignal = SignalObj(np.fft.irfft(Y, n_samples, axis=0), 'time', + samplingRate=signal.samplingRate, signalType=signal.signalType) + newSignal.channels = cp.deepcopy(signal.channels) + return newSignal + + +def extrapolate_low(signal, freq_range, slope=0.): + """ + Spectrum below freq_range replaced by an extrapolation, as the lower + part of ita_spk_extrapolate: the mean magnitude of freq_range, changed + by slope dB per decade towards low frequencies, and a linear phase from + the mean group delay of freq_range, joined to the phase at its lower + end. + + Parameters + ---------- + freq_range : (float, float) + Frequencies (Hz) whose magnitude and group delay are averaged. + slope : float + dB per decade below freq_range. Default 0. + + >>> fixed = pytta.extrapolate_low(ir, (80, 160)) + """ + x = np.asarray(signal.timeSignal, dtype=float) + n = x.shape[0] + X = np.fft.rfft(x, axis=0) + f = np.fft.rfftfreq(n, 1 / signal.samplingRate) + lo = int(np.flatnonzero(f > min(freq_range))[0]) + hi = int(np.flatnonzero(f > max(freq_range))[0]) + df = f[1] + band = X[lo:hi + 1] + gdelay = -np.diff(np.unwrap(np.angle(X), axis=0), axis=0) / (2 * np.pi * df) + tau = np.mean(gdelay[lo:hi + 1], axis=0) + fl = f[:lo + 1, None] + phase = -2 * np.pi * tau * fl + phase = phase - (phase[-1] - np.angle(band[0])) + with np.errstate(divide='ignore', invalid='ignore'): + shape = 10 ** (slope * (np.log10(fl) - np.log10(fl[-1])) / 20) + shape[0] = 0 if slope > 0 else shape[min(1, lo)] # DC: limit, or held for slope <= 0 + X[:lo + 1] = np.mean(np.abs(band), axis=0) * shape * np.exp(1j * phase) + X[0] = X[0].real + newSignal = SignalObj(np.fft.irfft(X, n, axis=0), 'time', + samplingRate=signal.samplingRate, signalType=signal.signalType) + newSignal.channels = cp.deepcopy(signal.channels) + return newSignal + + def fft_degree(*args,**kwargs): """ DEPRECATED diff --git a/tests/golden/make_psd_matlab.m b/tests/golden/make_psd_matlab.m new file mode 100644 index 0000000..f2db99d --- /dev/null +++ b/tests/golden/make_psd_matlab.m @@ -0,0 +1,12 @@ +% Golden data for tests/test_spectrum.py: MATLAB cpsd as ita_psd of the +% ITA-Toolbox calls it (symmetric Hann window of the block size, half-block +% overlap, one-sided density). The noise is rand with the Mersenne twister, +% which numpy.random.RandomState(1).random_sample rebuilds. +fs = 8000; n = 4096; block = 512; +rng(1, 'twister'); +x = rand(n, 1) - 0.5; +y = filter([1 0.5 0.25], 1, x); +w = hanning(block); +Pxx = cpsd(x, x, w, block/2, block, fs); +Pxy = cpsd(x, y, w, block/2, block, fs); +save('psd_matlab.mat', 'fs', 'block', 'Pxx', 'Pxy'); diff --git a/tests/golden/psd_matlab.mat b/tests/golden/psd_matlab.mat new file mode 100644 index 0000000000000000000000000000000000000000..016d857b62806b49ab09ac717b347894584aa41c GIT binary patch literal 6258 zcma){H$Ce5&X@?d)xFa8<@c2gPT`X<6p`3W*|Vj_l73AHpNjX z4M3fEQxe-mNUcbeif-dIzS=B$n8CAXw|&3% zDjx^F2%c~g1XSRm0=?QVNk2|;2wnim*XIwZZHb+_3yvHQWp*_2SVcW0yF7P`v>)SF zv#4mR^geuh5z{CdK7be^M7cB4zt~DrW36|zHSfMnjC3+8QEth`^cB?8izHP!a}e{$ zkF6ybOCYXIRAt%MRjR*m8vJuR!=|({DXWu^-1!0|3jM2J$uz9xCxzylERUKY+g&<> zEtgRNi3entaJ{-1bS||4Ap2Q^H8}*JNn4L2iwjJN*af^~J7k-B7cMm|28+(J7e5m$ zv@MIhx>{x?sVLm4K_{X`WAaqbwj`)am^ufBeH@|InlQW|QSL>cdw3v9=4qaGjMbce z5Z7)|ch`j@$+T%Kmvm6J#qth5bo+iH*N(?Hix(aD)MS@-iFtj?|Jchh`bnI!u@Wes zim82k^S6-xjb3G&G-LGVSH0jcHA|Xb&DP~tRDfAv-t?THdIIsYK3VCxHfA`rozM8? z)d{B1FWWc3a$VSQXaGI>mg#ZUe3w-_M{#|eXBn3vI4j@JFnqtKrR<-uT9K)^XvHCB zu9j!tBov`y-~L|UVMUF>l5DHqER#PDzkZL|HI8&kOKv$4H;I0cWx{WwEMC*8k8u+A z*ip5qwh?WzWf)v=U zJyp>od~Ex7YCg1sg5UxYi6J=!RoK;Utjsqg^&wUm7 z`8S3OvDMyF5Q*PcdM)nO5@}KH!7wr5n81n#ExB63+V3Z7IHUUbv#9NrL%CCk6YR)j}p1a^8)rCLbX2>|Nrx-~YEW+dht zO-N+$T214aO78U9?9%n?aRA<=<=FYeLRJNfITebaDBcV^qrQkj4<+eWUZQ^gjV9YJ zWD36~Pdi)fnCi{+$LiRxh})U-;!Y%>KEy}qJsy$%-Na62%-5eBaFz&msq22YDL929 zcH_D>+>yS0A>YsWM22o>Gprp=0`*tjvE|V+xeUvl;Y}PVy_)OZ3QsJJqIci(+kkR; z_lZ3hzE+aW_q?wQ!r*s#a1waQZ&Y1AvOCJ2Z3fS?ev%H-H{TcD!X3Y+Hc;3qJa_3J z&1*aTKwAkS!M)nYGKtzp;5&O`5m{tmpDCjeFZ_w$;A2JvPDH@6kIyB%-HT`EbXjxI0!CB?U8{4b^w3s1Mb z-IjE+81orC{;B zbuTv^tt_TzU9Ef>Rmv7+J0<>_VNx1zy?XXkvF_8Yf9EK&v^PjjX{29{kB@WKjw*-m z_-`xhad8`-dJlVk5BZXt`O@px``b1tvPaN?legPlQaBGQ=u zkd?U((fQA@;!kXZX{B|l&R`hV$xgvWp%T$1u;~dUlxf>r#ZTL2LjHa2pNgGr%$b^! zO!YF^cK4+PEqiCw)bD(GyRyh|pqWpS*@%3vlXri8?%-;Ww{*n3f%|M4PtWM1lTfMB zz>E}D^eUy%ocWFjgimqHojxl22$I6?8?R}6Wuw!?N-ioo*;6E{LEoC@(x`2L(?!yWPrLKZ4cf1sny5KPrXK~6@&3V=zr+U78KbJfB`3~H3ZEwE#ND5c%tS5_3wSZ~Io!+{iclgeO zU%ENhw_Ody(gnZJkyow^K<=L@*}vobbp3>Ze)8d2CfrdotwSI|rmWU&t0>~#J2-EV zK^_tPp>!j%RDxj_T6Z_1fg!PCBYBfZJ%BhA zn>_Z9_`YLPwK@8OMbZ`}C)Ji}3u1QbH_4=WF2a^~BD%TEtU3Az#m2Uz_vy(yp>0!> zl3<;f6#NrfxowZ4ahA3D$E_82YyWOxuSb)fDa3c={4e5*F=3&CD^)taoum zB|IfKl&sccYQvF7;9Ib%?>gA)A<*f04eXRZ_3yE!dEO}B>*mE-nJCt#DNT+jb0q4L zjXCkyqeMjiyyf++Wq(i#n6=I_e-n@K64L!rx7=^6`kmh|P%_FHgF%gqIv@)cI)&v&4=F7v)72{EIq zCKZ(o8YExy0HXn%^hvx|U7fATF?huQx}7>@2o0fd=5Nw>YOE?8&}l!Vy+g#yL?YoE8P1!bk4CG=^!$aDH765n5SgZ_c!!jj3G0YFHV7=Vxs&X=cS9Jl=3h$!O;MK}kMk&-|F81U`B2bYm zjonmjwyb@cfX96kNK1En*$w_qX*&-x?w(A{SeNz8Ie*xNsU@MTYD^#y6~@>F7$Y^k z3Y^DSUHM3#2V-o}pP!-|<+GSEcfu3QL6qMzFyA<%s!7+T^}u0;Sl;ye{eN?h-}L-61USX zFwfJQ95TCUCx`>Byw*djg!(Daj`(${cc`67RxCYY{Aq-{tC(YUgM zTsD1O0Uz|LQ*ot962--Z5^e}QZ?Gju`E0A%^oU4lv#qlk@wlZOYsu_fccf%M?2Bwk zNjX(wN1JOA6}~g;8gny$@>?eU`F+4e;kY;Zp+?0@(puf3q)q46>A&H_fq@_?Ra&5+ z@?oz0aem6xIMZ`BnS#uO{$BU;&oAjKoA?*M`iY^FnNNi3(s;-vu!PF&P|Ig?{VF`b3+tMLPSRXxksZYg*?WzB3}b$!F`8&_0(a}s{J*AAd> zMkF_VU7Vif`n=sTAqX&%ZSRf_6tk`YLjEm_8}vwzRuV|Bs#_=f2WnK=;Zy&HRaxg> zHLjnX$UB?UVE}&;+I{lmlcPCH6-=NRNX!n^$jb=L5AGo*o9Yh^^C6uFq5W3+9wLW5jysW_vPEV6&+?st={0MmN5^1z-B9ercb%58CQ|Qv@2v$|x zL`3?PpU}Cn@oVD3`}vyZJ%kiKe%HKClmwP_vUQt8kLoUc+X2o%5-e&^O();(9QDd| z6_MSP8+kvl47Ma92jC{r3WzT5Q@Q}E9H^t74|E!Q);FU3j?)e?I414(tG%P$L1rY% zD=!d=2bDsvqb_vI6j1%F&Y(iVk3RHfV+nnc3m%oo$~!mOyQ#E9=;$lmkahn+{nKgO zUI80XMSag5&KjB+hCz;sba5ksQnd|2bItOPk>LRVrDqIvhCPPY$>L9y0)oRB6k!`k zN20Ah&_wflU1;Q2=lqxHfMsRytGjv-<0N|IC9eyt;D_RmVTMM$jE2qwpKg)|>b%)S zjNk~@u9xz_Z)+PFL+(qOrciH}dY79}r+(V3tRf3_3>7^{x2w*lC||V)r-KbL=gvme zWrlj}qBY@n7h&MEE-Ql`@e5*0S{8D#mXJ*XBQ@sL3VdZYNc5^xMBGLlKBc%>L(OE6 zB3(e^>QQU$KqtqbqDIFo>#g-Q#67qcn-VQvSC4FRA+)qG)GmgXa%p_4`4ZeluO%oZ zzE0{7CO^#7FcM!>i9KiN)?61h0Gp61DXwgE6=JA8Bs0lNE$GO$rOBNQ?K`85@n3(QR*22D{3^qiLyhwTQG{i(a806O zi+}Q?DIIZ|;1?I}kVAZ9|34z!zdZCvFuZs}^#n{UKL;j4Bo#C9RXdL*y*LK$hAyix z$I1#a9<3Xp{ZGMHo0ASaqYgrDAkO{1&#h18cd`!6S>Wu~CebI7LH4vb_1N z>|Uw;q$p*;Ui5ijp@FIk+FSS@Tk(Q(p{3Eh(TVp0<)Z%jTGdk~vW`x!FWK00yKfPg znJ(kY8pfpU2tTjonMMom!5}C_X&(NP6u}j`s{O}a9O*A9m1R=`k_FZ@=BF%`(u6PRl+ zZ*bVH^0n3^_!9;H`4elke~6w7-)qKg(>cQJI#tW^QtY2!d&V%1lzoTGBjWykefU^@ z#E_9nV`sJQr;X6ss|No(1-nIo^s3h$gO|Xd*nnJKq-U?_F8d!R6=}KcG9mBFkOl#L zGiP=t*`vZAuP;bgG@p2O)^aGM<`aD%m@m)?I75WDG}+%;eSchm?Y`gls8+Rqc#^po zJRw3*;3k-DgTX{-S{QbVmH3VX*uiII4q0*ASEBwdo+N42LM956*B1R^u1$vkGz3CG zohp$8Grd#pkD2KDtchpn?3eHlk2~1@%{G__qpZz?XY3CeFplfgBZY{2oy{{iEsi=c zxWT*@TR}yVhk4gq?{Q(nbfd*wGeiUX_v!-fDI^mz2?HBu=m%G0ABBaPHl%i{2ohc3 zk?X~`6iY`0lvw^z=@Q)E*`vwey*A0mHk~PyC1by#cu_v7g(S7HP@X1Xq@D z%thPa`6N^u;$Q0oHgusw-C03l83D7oT7i2yW{G_tAJxH`Ei&-BLn&uD+-1>+fL$tT zk7`9%#J^GnLo&Ut6W;o$G*=VvR1e-+;smP7*hB}GfD-Of{qTQjRiMe*j}GFL``SP1 za84s)9f101HxDmFJ{q+HhMt-DZsuq>Ln$#^X{?&|_+9R!tV3og@mb&GDD$>Qt7Fp8WnN4=Q;P$$aj7LSzw|I-?;2YEEsm&J#)S*O0P3#kMBNWxe!v z-#S1ztfLY32!lNcK7599O!w%h~N2H?_kj@Vr`ewU#m0&_IRa$SwD;2tz z>y`19HF<*e)rn8*C!fRG>fOp)vFqhL4l+|G&@mzx{r0Vp(Hu6cp@a PpJaOrB&4?~Z;Jm1