[aguirreSL]

Adds to pytta.functions three spectrum utilities that the ITA-Toolbox has:

1. psd(signal, signal2=None, block_size=None, fft_size=None): the power spectral density of each channel, or the cross power spectral density with signal2, as ita_psd. It uses Welch averaging with a symmetric Hann window of the block size (the whole signal by default) and half-block overlap, one-sided, in units²/Hz. The cross spectrum follows the convention of MATLAB cpsd, which ita_psd calls, X conj(Y). That is the conjugate of scipy.signal.csd.
2. interpolate_spectrum(signal, n_samples, method='spline', abs_phase=False): the DFT interpolated onto the bins of another length, as ita_interpolate_spk. The interpolation is by spline, as MATLAB interp1, or linear, on real and imaginary parts or on magnitude and unwrapped phase. It is meant for impulse responses.
3. extrapolate_low(signal, freq_range, slope=0): the spectrum below freq_range replaced by its mean magnitude there, changed by slope dB per decade. The phase is linear, from the mean group delay of the range, and joins the phase at its lower end. This is the lower part of ita_spk_extrapolate, whose upper part is not implemented in the ITA-Toolbox either.

Four other spectrum functions of the ITA-Toolbox are not included here. ita_split_frequencies relies on the channel naming of the ITA-Toolbox. ita_extrapolate needs a rational fit of poles and zeros. ita_smooth_notches and ita_lifter need smooth_frequency (#50), the cepstrum (#80), start_ir (#96) and time_window (#49), so they can follow once those are merged. ita_decimate_spk was also left out, since what its phase reconstruction should preserve is not clear from the code.

The tests in tests/test_spectrum.py print the measured differences:

1. The power and cross power spectral densities of a noise and of the noise through a short FIR agree with MATLAB cpsd, called as ita_psd calls it, within 3.4e-8 relative (golden data from tests/golden/make_psd_matlab.m). The integral of the density is the variance within 0.4 %.
2. The 256-point DFT of a one-pole filter a^n, interpolated onto the bins of 1024 points, follows its true response 1 / (1 - a e^(-jω)) within 1.4e-6. Interpolating onto the same length returns the signal.
3. The extrapolation keeps a pure delay (3.6e-16), and a slope of 6 dB per decade leaves the magnitude 6.0000 dB lower a decade below the range.

Checked with Python 3.12, numpy 1.26 and scipy 1.17. The 7 failures of development (3 in test_h5io.py, 4 in test_pyttasave.py) are the same with and without this change.

Co-Authored-By: Claude Fable 5.1 noreply@anthropic.com