audio_samples 2.0.0

//! Spectral analysis and frequency-domain transformations for [`AudioSamples`].
//!
//! This module implements the [`AudioTransforms`] trait, providing FFT-based
//! spectral analysis and spectrogram computation for audio data. All spectral
//! operations delegate to the [`spectrograms`] crate.
//!
//! Frequency-domain analysis is central to audio processing --- spectrograms,
//! MFCCs, chromagrams, and pitch detection all begin with a transform into the
//! spectral domain.
//!
//! All operations are available on any [`AudioSamples<T>`] where `T` is a
//! supported sample type (`u8`, `i16`, `I24`, `i32`, `f32`, `f64`).
//! Most spectral methods require mono input; multi-channel signals must be
//! mixed or channel-selected before use. The one exception is
//! [`fft`](AudioTransforms::fft), which transforms every channel
//! independently.
//!
//! ### Available operations
//!
//! | Method | Description | Multi-channel? |
//! |--------|-------------|:--------------:|
//! | [`fft`](AudioTransforms::fft) | Full-signal FFT | Yes |
//! | [`stft`](AudioTransforms::stft) / [`istft`](AudioTransforms::istft) | Forward / inverse STFT | No |
//! | [`linear_spectrogram`](crate::operations::AudioTransforms::linear_spectrogram)(AudioTransforms::linear_spectrogram) | Linearly-spaced spectrogram | No |
//! | [`log_frequency_spectrogram`](crate::operations::AudioTransforms::log_frequency_spectrogram)(AudioTransforms::log_frequency_spectrogram) | Log-Hz spectrogram | No |
//! | [`mel_spectrogram`](crate::operations::AudioTransforms::mel_spectrogram)(AudioTransforms::mel_spectrogram) | Mel-scaled spectrogram | No |
//! | [`mfcc`](AudioTransforms::mfcc) | Mel-Frequency Cepstral Coefficients | No |
//! | [`chromagram`](AudioTransforms::chromagram) | Chroma (pitch-class) features | No |
//! | [`gammatone_spectrogram`](crate::operations::AudioTransforms::gammatone_spectrogram)(AudioTransforms::gammatone_spectrogram) | Gammatone (auditory) spectrogram | No |
//! | [`constant_q_transform`](AudioTransforms::constant_q_transform) | Constant-Q Transform | No |
//! | [`cqt_spectrogram`](crate::operations::AudioTransforms::cqt_spectrogram)(AudioTransforms::cqt_spectrogram) | CQT spectrogram | No |
//! | [`power_spectral_density`](AudioTransforms::power_spectral_density) | Welch PSD estimate | No |
//! | [`magphase`](AudioTransforms::magphase) | Magnitude / phase decomposition | — |
//!
//! ## Example
//!
//! ```rust
//! use audio_samples::{AudioSamples, AudioTransforms, sample_rate};
//! use ndarray::array;
//! use std::num::NonZeroUsize;
//!
//! let data  = array![1.0f32, 0.0, -1.0, 0.0];
//! let audio = AudioSamples::new_mono(data, sample_rate!(44100)).unwrap();
//!
//! let spectrum = audio.fft(NonZeroUsize::new(4).unwrap()).unwrap();
//! assert_eq!(spectrum.shape()[0], 1); // one row per channel
//! ```
//!
//! ## Error Handling
//!
//! All operations return [`crate::AudioSampleResult`]. Mono-only methods
//! return [`crate::AudioSampleError::Layout`]
//! ([`LayoutError::ChannelCountUnsupported`])
//! when given multi-channel input. Failures originating inside [`spectrograms`]
//! are propagated as [`crate::AudioSampleError::Processing`].
//!
//! ## See Also
//!
//! - [`AudioSamples`]: The core audio data type.
//! - [`AudioTransforms`]: The trait defining all spectral operations.
//!
//! [`AudioTransforms`]: crate::operations::traits::AudioTransforms
use std::num::{NonZeroU32, NonZeroUsize};

use crate::operations::traits::AudioTransforms;
use crate::{
    AudioSampleError, AudioSampleResult, AudioSamples, AudioTypeConversion, ChannelRequirement,
    LayoutError, ParameterError, StandardSample,
};

use ndarray::{Array1, Array2, Axis};
use non_empty_slice::NonEmptySlice;
use num_complex::Complex;
use spectrograms::{
    AmpScaleSpec, ChromaParams, Chromagram, CqtParams, CqtResult, CqtSpectrogram, FftPlanner,
    Gammatone, GammatoneParams, GammatoneSpectrogram, LinearHz, LinearSpectrogram, LogHz,
    LogHzParams, LogHzSpectrogram, LogParams, MelParams, MelSpectrogram, Mfcc, MfccParams,
    Spectrogram, SpectrogramParams, StftParams, StftPlan, StftResult, WindowType, fft_convolve,
    fft_deconvolve,
};

use std::cell::RefCell;

thread_local! {
    /// Thread-local [`FftPlanner`] reused across all FFT-based transforms in
    /// this module. `spectrograms::FftPlanner` wraps a real-FFT planner that
    /// memoizes plans (twiddles, scratch) by size internally, so reconstructing
    /// it per call discarded that cache. Keeping one planner per thread lets
    /// repeated same-size transforms reuse cached plans. The planner is `!Sync`
    /// (it owns a real-FFT planner), so `thread_local` is the correct scope.
    static FFT_PLANNER: RefCell<FftPlanner> = RefCell::new(FftPlanner::new());
}

/// Runs `f` with the thread-local cached [`FftPlanner`] borrowed mutably.
///
/// The borrow is released as soon as `f` returns, so callers must finish all
/// planner use inside the closure.
#[inline]
fn with_fft_planner<R>(f: impl FnOnce(&mut FftPlanner) -> R) -> R {
    FFT_PLANNER.with(|p| f(&mut p.borrow_mut()))
}

/// Power spectral density estimate produced by [`AudioTransforms::power_spectral_density`].
///
/// Pairs a frequency axis with the estimated power-per-Hz at each bin. The two
/// slices are always the same length: `frequencies()[i]` is the centre frequency
/// of bin `i` in Hz, and `density()[i]` is the estimated power spectral density
/// at that frequency.
#[derive(Clone, Debug, PartialEq)]
pub struct Psd {
    frequencies: Vec<f64>,
    density: Vec<f64>,
}

impl Psd {
    /// Constructs a [`Psd`] from a frequency axis and a matching density vector.
    #[inline]
    pub(crate) fn new(frequencies: Vec<f64>, density: Vec<f64>) -> Self {
        Self {
            frequencies,
            density,
        }
    }

    /// The frequency axis, in Hz. Same length as [`density`](Self::density).
    #[inline]
    #[must_use]
    pub fn frequencies(&self) -> &[f64] {
        &self.frequencies
    }

    /// The estimated power spectral density (power per Hz) at each frequency bin.
    /// Same length as [`frequencies`](Self::frequencies).
    #[inline]
    #[must_use]
    pub fn density(&self) -> &[f64] {
        &self.density
    }

    /// Consumes the estimate, returning the `(frequencies, density)` vectors.
    #[inline]
    #[must_use]
    pub fn into_parts(self) -> (Vec<f64>, Vec<f64>) {
        (self.frequencies, self.density)
    }
}

impl<T> AudioTransforms for AudioSamples<'_, T>
where
    T: StandardSample,
{
    /// Computes the Fast Fourier Transform of the audio signal.
    ///
    /// Each channel is transformed independently. The output has one row per
    /// channel containing the complex spectral bins.
    ///
    /// # Arguments
    /// - `n_fft` — FFT length in samples. If longer than the signal the input
    ///   is zero-padded internally.
    ///
    /// # Returns
    /// An `Array2<Complex<f64>>` where each row is the FFT of the
    /// corresponding channel.
    ///
    /// # Errors
    /// Returns an error if the FFT computation fails.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, sample_rate};
    /// use ndarray::array;
    /// use std::num::NonZeroUsize;
    ///
    /// let data  = array![1.0f32, 0.0, -1.0, 0.0];
    /// let audio = AudioSamples::new_mono(data, sample_rate!(44100)).unwrap();
    /// let spectrum = audio.fft(NonZeroUsize::new(4).unwrap()).unwrap();
    /// assert_eq!(spectrum.shape()[0], 1); // one row per channel
    /// ```
    #[inline]
    fn fft(&self, n_fft: NonZeroUsize) -> AudioSampleResult<Array2<Complex<f64>>> {
        let working_samples = self.to_format::<f64>();

        let mut channel_ffts: Vec<Array1<Complex<f64>>> =
            Vec::with_capacity(self.num_channels().get() as usize);

        with_fft_planner(|fft_planner| -> AudioSampleResult<()> {
            for ch in working_samples.channels() {
                let ch_slice = ch
                    .as_slice()
                    .expect("Can always get a mono channel as a slice");
                // safety: ch_slice is non-empty since audio has samples
                let working_ch_slice = unsafe { NonEmptySlice::new_unchecked(ch_slice) };
                let channel_fft: Array1<Complex<f64>> = fft_planner.fft(working_ch_slice, n_fft)?;
                channel_ffts.push(channel_fft);
            }
            Ok(())
        })?;

        ndarray::stack(
            Axis(0),
            &channel_ffts.iter().map(|a| a.view()).collect::<Vec<_>>(),
        )
        .map_err(Into::into)
    }

    /// Computes the Short-Time Fourier Transform (STFT) of a mono signal.
    ///
    /// The signal is divided into overlapping, windowed frames and each frame
    /// is transformed to the frequency domain. The returned [`StftResult`]
    /// carries both the complex matrix and the parameters required by
    /// [`istft`] for reconstruction.
    ///
    /// # Arguments
    /// - `params` — STFT configuration (FFT size, hop size, window function,
    ///   and centering behaviour). Field-level constraints are documented on
    ///   [`StftParams`].
    ///
    /// # Returns
    /// An [`StftResult`] containing the complex STFT matrix and metadata.
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying [`spectrograms`] STFT computation.
    ///
    /// [`istft`]: Self::istft
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let params = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let result = audio.stft(&params).unwrap();
    /// assert!(result.data.nrows() > 0); // frequency bins
    /// ```
    #[inline]
    fn stft(&self, params: &StftParams) -> AudioSampleResult<StftResult> {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "stft",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }

        let samples = self.to_format::<f64>();
        let samples_slice = samples
            .as_slice()
            .expect("Safe since we have ensured mono audio samples");
        // safety: samples_slice is non-empty since audio has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        let spect_params = SpectrogramParams::new(params.to_owned(), self.sample_rate_hz())?;
        let mut stft_plan = StftPlan::new(&spect_params)?;

        stft_plan
            .compute(samples_slice, &spect_params)
            .map_err(Into::into)
    }

    /// Reconstructs a time-domain signal from an [`StftResult`].
    ///
    /// Uses overlap-add synthesis with the window and hop parameters stored
    /// in the [`StftResult`]. The output is a mono signal at the sample rate
    /// that was recorded during the forward transform.
    ///
    /// # Arguments
    /// - `stft` — the [`StftResult`] produced by a prior call to [`stft`].
    ///
    /// # Returns
    /// A mono [`AudioSamples`] at the original sample rate.
    ///
    /// # Errors
    /// Returns an error if the reconstruction fails (e.g. mismatched
    /// parameters inside the [`StftResult`]).
    ///
    /// [`stft`]: Self::stft
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let params = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let stft_result = audio.stft(&params).unwrap();
    /// let reconstructed = AudioSamples::<f64>::istft(stft_result).unwrap();
    /// assert!(reconstructed.samples_per_channel().get() > 0);
    /// ```
    #[inline]
    fn istft(stft: StftResult) -> AudioSampleResult<AudioSamples<'static, T>> {
        let sample_rate = stft.sample_rate;
        // safety: sample_rate is non-zero since it was validated during STFT computation on the ``spectrograms`` side
        let sample_rate = unsafe { NonZeroU32::new_unchecked(sample_rate as u32) };
        let data = &stft.data;
        let n_fft = stft.params.n_fft();
        let hop_size = stft.params.hop_size();
        let window = stft.params.window();
        let centre = stft.params.centre();
        let data = spectrograms::istft(data, n_fft, hop_size, window, centre)?;

        let audio = AudioSamples::from_mono_vec::<f64>(data, sample_rate);
        Ok(audio)
    }

    /// Computes a linearly-spaced spectrogram.
    ///
    /// Prefer the typed convenience methods —
    /// [`linear_magnitude_spectrogram`](crate::operations::AudioTransforms::linear_magnitude_spectrogram), [`linear_power_spectrogram`](crate::operations::AudioTransforms::linear_power_spectrogram), or
    /// [`linear_db_spectrogram`](crate::operations::AudioTransforms::linear_db_spectrogram) — for the most common amplitude scales.
    ///
    /// # Arguments
    /// - `params` — spectrogram parameters (window, hop, FFT size).
    /// - `db` — required when `AmpScale` is `Decibels`; ignored otherwise.
    ///
    /// # Returns
    /// A `Spectrogram<LinearHz, AmpScale>`.
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying spectrogram computation.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{Magnitude, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let spect = audio.linear_spectrogram::<Magnitude>(&params, None).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn linear_spectrogram<AmpScale>(
        &self,
        params: &SpectrogramParams,
        db: Option<&LogParams>,
    ) -> AudioSampleResult<Spectrogram<LinearHz, AmpScale>>
    where
        AmpScale: AmpScaleSpec + 'static,
    {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "linear_frequency_spectrogram",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }
        let samples = self.cast_as_f64();
        let samples_slice = samples.as_slice().expect("Safe since we have ensured mono");
        // safety: Guaranteed non-empty since audio is mono and has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        LinearSpectrogram::<AmpScale>::compute(samples_slice, params, db).map_err(Into::into)
    }

    /// Computes a log-frequency-spaced spectrogram.
    ///
    /// Prefer the typed convenience methods —
    /// [`loghz_power_spectrogram`](crate::operations::AudioTransforms::loghz_power_spectrogram), [`loghz_magnitude_spectrogram`](crate::operations::AudioTransforms::loghz_magnitude_spectrogram), or
    /// [`loghz_db_spectrogram`](crate::operations::AudioTransforms::loghz_db_spectrogram) — for the most common amplitude scales.
    ///
    /// # Arguments
    /// - `params` — spectrogram parameters (window, hop, FFT size).
    /// - `loghz` — log-Hz frequency-axis configuration (min/max
    ///   frequency, number of bins). Field-level constraints are
    ///   documented on [`LogHzParams`].
    /// - `db` — required when `AmpScale` is `Decibels`; ignored otherwise.
    ///
    /// # Returns
    /// A `Spectrogram<LogHz, AmpScale>`.
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying spectrogram computation.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{LogHzParams, Magnitude, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let loghz = LogHzParams::new(nzu!(64), 80.0, 8_000.0).unwrap();
    /// let spect = audio.log_frequency_spectrogram::<Magnitude>(&params, &loghz, None).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn log_frequency_spectrogram<AmpScale>(
        &self,
        params: &SpectrogramParams,
        loghz: &LogHzParams,
        db: Option<&LogParams>,
    ) -> AudioSampleResult<Spectrogram<LogHz, AmpScale>>
    where
        AmpScale: AmpScaleSpec,
    {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "log_frequency_spectrogram",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }
        let samples: AudioSamples<'static, f64> = self.cast_as_f64();
        let samples_slice = samples.as_slice().expect("Safe since we have ensured mono");
        // safety: Guaranteed non-empty since audio is mono and has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        LogHzSpectrogram::<AmpScale>::compute(samples_slice, params, loghz, db).map_err(Into::into)
    }

    /// Computes a mel-scaled spectrogram.
    ///
    /// The mel scale approximates human auditory perception by compressing
    /// high frequencies relative to low ones.  Prefer the typed
    /// convenience methods — [`mel_mag_spectrogram`](crate::operations::AudioTransforms::mel_mag_spectrogram),
    /// [`mel_power_spectrogram`](crate::operations::AudioTransforms::mel_power_spectrogram), or [`mel_db_spectrogram`](crate::operations::AudioTransforms::mel_db_spectrogram) — for the
    /// most common amplitude scales.
    ///
    /// # Arguments
    /// - `params` — spectrogram parameters (window, hop, FFT size).
    /// - `mel` — mel filter-bank configuration (number of bands,
    ///   frequency range). Field-level constraints are documented on
    ///   [`MelParams`].
    /// - `db` — required when `AmpScale` is `Decibels`; ignored otherwise.
    ///
    /// # Returns
    /// A [`MelSpectrogram<AmpScale>`].
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying spectrogram computation.
    ///
    /// ## See Also
    /// - [Mel scale — Wikipedia](https://en.wikipedia.org/wiki/Mel_scale)
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{Magnitude, MelParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let mel = MelParams::new(nzu!(40), 0.0, 8_000.0).unwrap();
    /// let spect = audio.mel_spectrogram::<Magnitude>(&params, &mel, None).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn mel_spectrogram<AmpScale>(
        &self,
        params: &SpectrogramParams,
        mel: &MelParams,
        db: Option<&LogParams>, // only used when AmpScale = Decibels
    ) -> AudioSampleResult<MelSpectrogram<AmpScale>>
    where
        AmpScale: AmpScaleSpec,
    {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "mel_spectrogram",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }
        let samples = self.cast_as_f64();
        let samples_slice = samples.as_slice().expect("Safe since we have ensured mono");
        // safety: Guaranteed non-empty since audio is mono and has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        MelSpectrogram::<AmpScale>::compute(samples_slice, params, mel, db).map_err(Into::into)
    }

    /// Computes Mel-Frequency Cepstral Coefficients (MFCCs).
    ///
    /// MFCCs are a compact spectral representation widely used in speech
    /// recognition and audio classification.  They are derived by
    /// applying a DCT to log-mel filter-bank energies.
    ///
    /// # Arguments
    /// - `stft_params` — STFT configuration used for the underlying
    ///   spectrogram.
    /// - `n_mels` — number of mel filter-bank bands.
    /// - `mfcc_params` — MFCC-specific configuration (number of
    ///   coefficients, etc.). Field-level constraints are documented on
    ///   [`MfccParams`].
    ///
    /// # Returns
    /// An [`Mfcc`] containing the MFCC matrix.
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying computation.
    ///
    /// ## See Also
    /// - [MFCC — Wikipedia](https://en.wikipedia.org/wiki/Mel-frequency_cepstral_coefficients)
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{MfccParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let mfcc_params = MfccParams::speech_standard();
    /// let result = audio.mfcc(&stft, nzu!(40), &mfcc_params).unwrap();
    /// assert!(result.data.nrows() > 0); // MFCC coefficients
    /// ```
    #[inline]
    fn mfcc(
        &self,
        stft_params: &StftParams,
        n_mels: NonZeroUsize,
        mfcc_params: &MfccParams,
    ) -> AudioSampleResult<Mfcc> {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "mfcc",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }

        let samples = self.cast_as_f64();
        let samples_slice = samples.as_slice().expect("Safe since we have ensured mono");
        // safety: Guaranteed non-empty since audio is mono and has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        let sample_rate_f = self.sample_rate_hz();
        spectrograms::mfcc(
            samples_slice,
            stft_params,
            sample_rate_f,
            n_mels,
            mfcc_params,
        )
        .map_err(Into::into)
    }

    /// Computes chromagram (pitch-class energy) features.
    ///
    /// A chromagram projects the spectrum onto the twelve pitch classes
    /// (C, C♯, D, … , B), collapsing octave differences.  The result
    /// is useful for harmonic and key detection.
    ///
    /// # Arguments
    /// - `stft_params` — STFT configuration used for the underlying
    ///   spectrogram.
    /// - `cfg` — chromagram configuration (tuning, normalization,
    ///   etc.). Field-level constraints are documented on
    ///   [`ChromaParams`].
    ///
    /// # Returns
    /// A [`Chromagram`].
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying computation.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{ChromaParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let chroma_params = ChromaParams::music_standard();
    /// let result = audio.chromagram(&stft, &chroma_params).unwrap();
    /// assert_eq!(result.data.nrows(), 12); // twelve pitch classes (C through B)
    /// ```
    #[inline]
    fn chromagram(
        &self,
        stft_params: &StftParams,
        cfg: &ChromaParams,
    ) -> AudioSampleResult<Chromagram> {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "chromagram",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }

        let samples = self.cast_as_f64();
        let samples_slice = samples.as_slice().expect("Safe since we have ensured mono");
        // safety: Guaranteed non-empty since audio is mono and has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        let sample_rate_f = self.sample_rate_hz();
        spectrograms::chromagram(samples_slice, stft_params, sample_rate_f, cfg).map_err(Into::into)
    }

    /// Computes a gammatone-filtered spectrogram.
    ///
    /// Gammatone filters model the bandpass response of the human
    /// cochlea.  The filter centre frequencies are spaced according to
    /// the ERB (Equivalent Rectangular Bandwidth) scale.  Prefer the
    /// typed convenience methods —
    /// [`gammatone_magnitude_spectrogram`](crate::operations::AudioTransforms::gammatone_magnitude_spectrogram),
    /// [`gammatone_power_spectrogram`](crate::operations::AudioTransforms::gammatone_power_spectrogram), or
    /// [`gammatone_db_spectrogram`](crate::operations::AudioTransforms::gammatone_db_spectrogram) — for the most common amplitude
    /// scales.
    ///
    /// # Arguments
    /// - `params` — spectrogram parameters (window, hop, FFT size).
    /// - `gammatone_params` — gammatone filter-bank configuration
    ///   (number of bands, frequency range). Field-level constraints are
    ///   documented on [`GammatoneParams`].
    /// - `db` — required when `AmpScale` is `Decibels`; ignored otherwise.
    ///
    /// # Returns
    /// A `Spectrogram<Gammatone, AmpScale>`.
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying computation.
    ///
    /// ## See Also
    /// - [Gammatone filter — Wikipedia](https://en.wikipedia.org/wiki/Gammatone)
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{GammatoneParams, Magnitude, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let gammatone = GammatoneParams::new(nzu!(32), 80.0, 8_000.0).unwrap();
    /// let spect = audio.gammatone_spectrogram::<Magnitude>(&params, &gammatone, None).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn gammatone_spectrogram<AmpScale>(
        &self,
        params: &SpectrogramParams,
        gammatone_params: &GammatoneParams,
        db: Option<&LogParams>,
    ) -> AudioSampleResult<Spectrogram<Gammatone, AmpScale>>
    where
        AmpScale: AmpScaleSpec,
    {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "gammatone_spectrogram",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }

        let samples = self.cast_as_f64();
        let samples_slice = samples.as_slice().expect("Safe since we have ensured mono");
        // safety: Guaranteed non-empty since audio is mono and has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        GammatoneSpectrogram::<AmpScale>::compute(samples_slice, params, gammatone_params, db)
            .map_err(Into::into)
    }

    /// Computes the Constant-Q Transform (CQT) of a mono signal.
    ///
    /// The CQT uses a bank of bandpass filters whose centre frequencies
    /// are spaced logarithmically with a constant ratio
    /// Q = f / Δf.  This gives it the same frequency resolution as the
    /// musical scale, making it preferred for pitch and harmonic
    /// analysis.
    ///
    /// # Arguments
    /// - `params` — CQT parameters (frequency range, bins per octave,
    ///   etc.). Field-level constraints are documented on [`CqtParams`].
    /// - `hop_size` — hop length in samples between successive frames.
    ///
    /// # Returns
    /// A [`CqtResult`] containing the CQT matrix.
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying CQT computation.
    ///
    /// ## See Also
    /// - [Constant-Q transform — Wikipedia](https://en.wikipedia.org/wiki/Constant-Q_transform)
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::CqtParams;
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// // 12 bins per octave, 7 octaves, starting at 32.7 Hz (C1)
    /// let cqt_params = CqtParams::new(nzu!(12), nzu!(7), 32.7).unwrap();
    /// let result = audio.constant_q_transform(&cqt_params, nzu!(256)).unwrap();
    /// assert!(result.data.nrows() > 0);
    /// ```
    #[inline]
    fn constant_q_transform(
        &self,
        params: &CqtParams,
        hop_size: NonZeroUsize,
    ) -> AudioSampleResult<CqtResult> {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "constant_q_transform",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }

        let samples = self.cast_as_f64();
        let samples_slice = samples.as_slice().expect("Safe since we have ensured mono");
        // safety: Guaranteed non-empty since audio is mono and has samples
        let samples_slice = unsafe { NonEmptySlice::new_unchecked(samples_slice) };
        spectrograms::cqt(
            samples_slice,
            f64::from(self.sample_rate().get()),
            params,
            hop_size,
        )
        .map_err(Into::into)
    }

    /// Computes a CQT-based spectrogram.
    ///
    /// Applies the CQT to the signal and returns the result as a typed
    /// spectrogram.  Prefer the typed convenience methods —
    /// [`cqt_magnitude_spectrogram`](crate::operations::AudioTransforms::cqt_magnitude_spectrogram), [`cqt_power_spectrogram`](crate::operations::AudioTransforms::cqt_power_spectrogram), or
    /// [`cqt_db_spectrogram`](crate::operations::AudioTransforms::cqt_db_spectrogram) — for the most common amplitude scales.
    ///
    /// # Arguments
    /// - `params` — spectrogram parameters (window, hop, FFT size).
    /// - `cqt` — CQT parameters (frequency range, bins per octave,
    ///   etc.). Field-level constraints are documented on [`CqtParams`].
    /// - `db` — required when `AmpScale` is `Decibels`; ignored otherwise.
    ///
    /// # Returns
    /// A [`CqtSpectrogram<AmpScale>`].
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if the signal is multi-channel.
    /// - Errors from the underlying computation.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{CqtParams, Magnitude, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let cqt = CqtParams::new(nzu!(12), nzu!(7), 32.7).unwrap();
    /// let spect = audio.cqt_spectrogram::<Magnitude>(&params, &cqt, None).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn cqt_spectrogram<AmpScale>(
        &self,
        params: &SpectrogramParams,
        cqt: &CqtParams,
        db: Option<&LogParams>,
    ) -> AudioSampleResult<CqtSpectrogram<AmpScale>>
    where
        AmpScale: AmpScaleSpec,
    {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "cqt",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }

        let working_samples = self.as_f64();
        let working_samples_slice = working_samples
            .as_slice()
            .expect("Safe since we have ensured mono");

        // safety: working_samples_slice is non-empty since audio is mono and has samples
        let working_samples_slice = unsafe { NonEmptySlice::new_unchecked(working_samples_slice) };

        CqtSpectrogram::<AmpScale>::compute(working_samples_slice, params, cqt, db)
            .map_err(Into::into)
    }

    /// Estimates the power spectral density using Welch's method.
    ///
    /// The signal is split into overlapping segments; each is windowed
    /// with a Hanning window and FFT'd, and the resulting periodograms
    /// are averaged.  The final values are normalised to power per Hz.
    ///
    /// # Arguments
    /// - `window_size` — length of each segment in samples.  Must not
    ///   exceed the signal length.
    /// - `overlap` — fractional overlap between adjacent segments, in
    ///   the range `[0, 1)`.
    ///
    /// # Returns
    /// A [`Psd`] whose [`frequencies`](Psd::frequencies) and
    /// [`density`](Psd::density) slices are of equal length. `frequencies[i]`
    /// is the centre frequency of bin `i` in Hz; `density[i]` is the estimated
    /// power spectral density at that frequency.
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Parameter`] if the signal is
    ///   multi-channel, if `overlap` is outside `[0, 1)`, or if
    ///   `window_size` exceeds the signal length.
    ///
    /// ## See Also
    /// - [Welch's method — Wikipedia](https://en.wikipedia.org/wiki/Welch%27s_method)
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// let psd = audio.power_spectral_density(nzu!(1024), 0.5).unwrap();
    /// assert_eq!(psd.frequencies().len(), psd.density().len());
    /// assert!(!psd.frequencies().is_empty());
    /// ```
    #[inline]
    fn power_spectral_density(
        &self,
        window_size: NonZeroUsize,
        overlap: f64,
    ) -> AudioSampleResult<Psd> {
        if self.is_multi_channel() {
            return Err(AudioSampleError::Layout(
                LayoutError::channel_count_unsupported(
                    "power_spectral_density",
                    ChannelRequirement::Mono,
                    self.num_channels().get(),
                ),
            ));
        }
        if !(0.0..1.0).contains(&overlap) {
            return Err(AudioSampleError::Parameter(ParameterError::out_of_range(
                "overlap",
                overlap.to_string(),
                "0.0",
                "1.0",
                "overlap must be in [0, 1)",
            )));
        }
        if self.samples_per_channel() < window_size {
            return Err(AudioSampleError::Parameter(ParameterError::invalid_value(
                "window_size",
                "window_size must not exceed the signal length",
            )));
        }

        let samples = self.to_format::<f64>();
        let signal = samples
            .as_slice()
            .expect("Safe since we have ensured mono audio samples");

        let win = window_size.get();
        let hop = ((1.0 - overlap) * win as f64).floor().max(1.0) as usize;
        let sample_rate = self.sample_rate_hz();

        let mut sum: Vec<f64> = Vec::new();
        let mut segment_count = 0u64;

        with_fft_planner(|planner| -> AudioSampleResult<()> {
            let mut start = 0;
            while start + win <= signal.len() {
                let segment = &signal[start..start + win];
                // safety: segment length == win == window_size >= 1
                let segment_slice = unsafe { NonEmptySlice::new_unchecked(segment) };
                let power = planner.power_spectrum(
                    segment_slice,
                    window_size,
                    Some(WindowType::Hanning),
                )?;

                if sum.is_empty() {
                    sum = power.to_vec();
                } else {
                    for (acc, &val) in sum.iter_mut().zip(power.iter()) {
                        *acc += val;
                    }
                }
                segment_count += 1;
                start += hop;
            }
            Ok(())
        })?;

        // Average
        let count = segment_count as f64;
        for val in &mut sum {
            *val /= count;
        }

        // Convert power to power spectral density (power / Hz)
        let freq_resolution = sample_rate / win as f64;
        for val in &mut sum {
            *val /= freq_resolution;
        }

        // Build frequency axis
        let frequencies: Vec<f64> = (0..sum.len())
            .map(|i| i as f64 * sample_rate / win as f64)
            .collect();

        Ok(Psd::new(frequencies, sum))
    }

    /// Shorthand for [`linear_spectrogram`](crate::operations::AudioTransforms::linear_spectrogram) with `Magnitude` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let spect = audio.linear_magnitude_spectrogram(&params).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn linear_magnitude_spectrogram(
        &self,
        params: &SpectrogramParams,
    ) -> AudioSampleResult<spectrograms::LinearMagnitudeSpectrogram> {
        self.linear_spectrogram::<spectrograms::Magnitude>(params, None)
    }

    /// Shorthand for [`linear_spectrogram`](crate::operations::AudioTransforms::linear_spectrogram) with `Power` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let spect = audio.linear_power_spectrogram(&params).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn linear_power_spectrogram(
        &self,
        params: &SpectrogramParams,
    ) -> AudioSampleResult<spectrograms::LinearPowerSpectrogram> {
        self.linear_spectrogram::<spectrograms::Power>(params, None)
    }

    /// Shorthand for [`linear_spectrogram`](crate::operations::AudioTransforms::linear_spectrogram) with `Decibels` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{LogParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let db = LogParams::new(-80.0).unwrap();
    /// let spect = audio.linear_db_spectrogram(&params, &db).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn linear_db_spectrogram(
        &self,
        params: &SpectrogramParams,
        db: &LogParams,
    ) -> AudioSampleResult<spectrograms::LinearDbSpectrogram> {
        self.linear_spectrogram::<spectrograms::Decibels>(params, Some(db))
    }

    /// Shorthand for [`log_frequency_spectrogram`](crate::operations::AudioTransforms::log_frequency_spectrogram) with `Power` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{LogHzParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let loghz = LogHzParams::new(nzu!(64), 80.0, 8_000.0).unwrap();
    /// let spect = audio.loghz_power_spectrogram(&params, &loghz).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn loghz_power_spectrogram(
        &self,
        params: &SpectrogramParams,
        loghz: &LogHzParams,
    ) -> AudioSampleResult<spectrograms::LogHzPowerSpectrogram> {
        self.log_frequency_spectrogram::<spectrograms::Power>(params, loghz, None)
    }

    /// Shorthand for [`log_frequency_spectrogram`](crate::operations::AudioTransforms::log_frequency_spectrogram) with `Magnitude` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{LogHzParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let loghz = LogHzParams::new(nzu!(64), 80.0, 8_000.0).unwrap();
    /// let spect = audio.loghz_magnitude_spectrogram(&params, &loghz).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn loghz_magnitude_spectrogram(
        &self,
        params: &SpectrogramParams,
        loghz: &LogHzParams,
    ) -> AudioSampleResult<spectrograms::LogHzMagnitudeSpectrogram> {
        self.log_frequency_spectrogram::<spectrograms::Magnitude>(params, loghz, None)
    }

    /// Shorthand for [`log_frequency_spectrogram`](crate::operations::AudioTransforms::log_frequency_spectrogram) with `Decibels` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{LogHzParams, LogParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let loghz = LogHzParams::new(nzu!(64), 80.0, 8_000.0).unwrap();
    /// let db = LogParams::new(-80.0).unwrap();
    /// let spect = audio.loghz_db_spectrogram(&params, &loghz, &db).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn loghz_db_spectrogram(
        &self,
        params: &SpectrogramParams,
        loghz: &LogHzParams,
        db: &LogParams,
    ) -> AudioSampleResult<spectrograms::LogHzDbSpectrogram> {
        self.log_frequency_spectrogram::<spectrograms::Decibels>(params, loghz, Some(db))
    }

    /// Shorthand for [`mel_spectrogram`](crate::operations::AudioTransforms::mel_spectrogram) with `Magnitude` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{MelParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let mel = MelParams::new(nzu!(40), 0.0, 8_000.0).unwrap();
    /// let spect = audio.mel_mag_spectrogram(&params, &mel).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn mel_mag_spectrogram(
        &self,
        params: &SpectrogramParams,
        mel: &MelParams,
    ) -> AudioSampleResult<spectrograms::MelMagnitudeSpectrogram> {
        self.mel_spectrogram(params, mel, None)
    }

    /// Shorthand for [`mel_spectrogram`](crate::operations::AudioTransforms::mel_spectrogram) with `Decibels` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{LogParams, MelParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let mel = MelParams::new(nzu!(40), 0.0, 8_000.0).unwrap();
    /// let db = LogParams::new(-80.0).unwrap();
    /// let spect = audio.mel_db_spectrogram(&params, &mel, &db).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn mel_db_spectrogram(
        &self,
        params: &SpectrogramParams,
        mel: &MelParams,
        db: &LogParams,
    ) -> AudioSampleResult<spectrograms::LogMelSpectrogram> {
        self.mel_spectrogram(params, mel, Some(db))
    }

    /// Shorthand for [`mel_spectrogram`](crate::operations::AudioTransforms::mel_spectrogram) with `Power` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{MelParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let mel = MelParams::new(nzu!(40), 0.0, 8_000.0).unwrap();
    /// let spect = audio.mel_power_spectrogram(&params, &mel).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn mel_power_spectrogram(
        &self,
        params: &SpectrogramParams,
        mel: &MelParams,
    ) -> AudioSampleResult<spectrograms::MelPowerSpectrogram> {
        self.mel_spectrogram(params, mel, None)
    }

    /// Shorthand for [`gammatone_spectrogram`](crate::operations::AudioTransforms::gammatone_spectrogram) with `Magnitude` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{GammatoneParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let gammatone = GammatoneParams::new(nzu!(32), 80.0, 8_000.0).unwrap();
    /// let spect = audio.gammatone_magnitude_spectrogram(&params, &gammatone).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn gammatone_magnitude_spectrogram(
        &self,
        params: &SpectrogramParams,
        gammatone_params: &GammatoneParams,
    ) -> AudioSampleResult<spectrograms::GammatoneMagnitudeSpectrogram> {
        self.gammatone_spectrogram::<spectrograms::Magnitude>(params, gammatone_params, None)
    }

    /// Shorthand for [`gammatone_spectrogram`](crate::operations::AudioTransforms::gammatone_spectrogram) with `Power` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{GammatoneParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let gammatone = GammatoneParams::new(nzu!(32), 80.0, 8_000.0).unwrap();
    /// let spect = audio.gammatone_power_spectrogram(&params, &gammatone).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn gammatone_power_spectrogram(
        &self,
        params: &SpectrogramParams,
        gammatone_params: &GammatoneParams,
    ) -> AudioSampleResult<spectrograms::GammatonePowerSpectrogram> {
        self.gammatone_spectrogram::<spectrograms::Power>(params, gammatone_params, None)
    }

    /// Shorthand for [`gammatone_spectrogram`](crate::operations::AudioTransforms::gammatone_spectrogram) with `Decibels` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{GammatoneParams, LogParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(100), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let gammatone = GammatoneParams::new(nzu!(32), 80.0, 8_000.0).unwrap();
    /// let db = LogParams::new(-80.0).unwrap();
    /// let spect = audio.gammatone_db_spectrogram(&params, &gammatone, &db).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn gammatone_db_spectrogram(
        &self,
        params: &SpectrogramParams,
        gammatone_params: &GammatoneParams,
        db: &LogParams,
    ) -> AudioSampleResult<spectrograms::GammatoneDbSpectrogram> {
        self.gammatone_spectrogram::<spectrograms::Decibels>(params, gammatone_params, Some(db))
    }

    /// Shorthand for [`cqt_spectrogram`](crate::operations::AudioTransforms::cqt_spectrogram) with `Magnitude` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{CqtParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let cqt = CqtParams::new(nzu!(12), nzu!(7), 32.7).unwrap();
    /// let spect = audio.cqt_magnitude_spectrogram(&params, &cqt).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn cqt_magnitude_spectrogram(
        &self,
        params: &SpectrogramParams,
        cqt: &CqtParams,
    ) -> AudioSampleResult<spectrograms::CqtMagnitudeSpectrogram> {
        self.cqt_spectrogram::<spectrograms::Magnitude>(params, cqt, None)
    }

    /// Shorthand for [`cqt_spectrogram`](crate::operations::AudioTransforms::cqt_spectrogram) with `Power` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{CqtParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let cqt = CqtParams::new(nzu!(12), nzu!(7), 32.7).unwrap();
    /// let spect = audio.cqt_power_spectrogram(&params, &cqt).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn cqt_power_spectrogram(
        &self,
        params: &SpectrogramParams,
        cqt: &CqtParams,
    ) -> AudioSampleResult<spectrograms::CqtPowerSpectrogram> {
        self.cqt_spectrogram::<spectrograms::Power>(params, cqt, None)
    }

    /// Shorthand for [`cqt_spectrogram`](crate::operations::AudioTransforms::cqt_spectrogram) with `Decibels` amplitude scale.
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use spectrograms::{CqtParams, LogParams, SpectrogramParams, StftParams, WindowType};
    /// use std::time::Duration;
    ///
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(200), sample_rate!(44100), 0.8);
    /// let stft = StftParams::new(nzu!(1024), nzu!(256), WindowType::Hanning, true).unwrap();
    /// let params = SpectrogramParams::new(stft, audio.sample_rate_hz()).unwrap();
    /// let cqt = CqtParams::new(nzu!(12), nzu!(7), 32.7).unwrap();
    /// let db = LogParams::new(-80.0).unwrap();
    /// let spect = audio.cqt_db_spectrogram(&params, &cqt, &db).unwrap();
    /// assert!(spect.data().nrows() > 0);
    /// ```
    #[inline]
    fn cqt_db_spectrogram(
        &self,
        params: &SpectrogramParams,
        cqt: &CqtParams,
        db: &LogParams,
    ) -> AudioSampleResult<spectrograms::CqtDbSpectrogram> {
        self.cqt_spectrogram::<spectrograms::Decibels>(params, cqt, Some(db))
    }

    /// Linear convolution of this signal with `other` via FFT.
    ///
    /// See [`AudioTransforms::convolve`](crate::operations::AudioTransforms::convolve).
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if either signal is multi-channel.
    /// - Errors propagated from the underlying FFT.
    #[inline]
    fn convolve(&self, other: &Self) -> AudioSampleResult<AudioSamples<'static, Self::Sample>> {
        let result_f64 = convolve_or_deconvolve_mono(self, other, |a, b| Ok(fft_convolve(a, b)?))?;
        from_f64_mono::<T>(result_f64, self.sample_rate())
    }

    /// Regularised spectral-division deconvolution.
    ///
    /// See [`AudioTransforms::deconvolve`](crate::operations::AudioTransforms::deconvolve).
    ///
    /// # Errors
    /// - [`crate::AudioSampleError::Layout`] if either signal is multi-channel.
    /// - Errors propagated from the underlying FFT.
    #[inline]
    fn deconvolve(
        &self,
        denominator: &Self,
        regularization: f64,
    ) -> AudioSampleResult<AudioSamples<'static, Self::Sample>> {
        let result_f64 = convolve_or_deconvolve_mono(self, denominator, |a, b| {
            Ok(fft_deconvolve(a, b, regularization)?)
        })?;
        from_f64_mono::<T>(result_f64, self.sample_rate())
    }

    /// Decomposes a complex spectrogram into magnitude and phase.
    ///
    /// Given a complex matrix `D`, returns `(S, P)` such that
    /// `D = S * P` elementwise, where `S` contains magnitudes raised to
    /// `power` and `P` contains unit-magnitude complex phase factors.
    /// Bins where the magnitude is zero are assigned a phase of `1 + 0i`.
    ///
    /// # Arguments
    /// - `complex_spect` — the complex STFT or FFT matrix.
    /// - `power` — exponent applied to the magnitude values.  `None`
    ///   defaults to 1 (raw magnitude).
    ///
    /// # Returns
    /// `(magnitude, phase)` — the magnitude matrix (real-valued) and
    /// the phase matrix (complex unit-magnitude).
    ///
    /// # Examples
    /// ```
    /// use audio_samples::{AudioSamples, AudioTransforms, nzu, sample_rate, sine_wave};
    /// use std::time::Duration;
    ///
    /// // Use a short signal so n_fft (1024) is larger than the signal length (441 samples).
    /// let audio = sine_wave::<f64>(440.0, Duration::from_millis(10), sample_rate!(44100), 0.8);
    /// let spectrum = audio.fft(nzu!(1024)).unwrap();
    /// let (mag, phase) = AudioSamples::<f64>::magphase(&spectrum, None);
    /// assert_eq!(mag.shape(), phase.shape());
    /// ```
    #[inline]
    fn magphase(
        complex_spect: &Array2<Complex<f64>>,
        power: Option<NonZeroUsize>,
    ) -> (Array2<f64>, Array2<Complex<f64>>) {
        // Magnitude: elementwise absolute value

        let mut mag = complex_spect.mapv(num_complex::Complex::norm);

        // zeros_to_ones: 1.0 where mag == 0, else 0.0
        let zeros_to_ones = mag.mapv(|x| if x == 0.0 { 1.0 } else { 0.0 });

        // mag_nonzero = mag + zeros_to_ones
        let mag_nonzero = &mag + &zeros_to_ones;

        // Compute phase = D / mag_nonzero, but handle zeros separately
        let mut phase = complex_spect.clone();

        let power = power.map_or(1.0, |p| p.get() as f64);

        // Perform elementwise division for real and imaginary parts
        ndarray::Zip::from(&mut phase)
            .and(&mag_nonzero)
            .and(&zeros_to_ones)
            .for_each(|p, &m_nz, &z| {
                let div = Complex {
                    re: p.re / m_nz + z, // add 1.0 if originally zero
                    im: p.im / m_nz,
                };
                *p = div;
            });

        // Raise magnitude to the given power
        mag.mapv_inplace(|x| x.powf(power));

        (mag, phase)
    }
}

/// Extracts the single mono channel of `audio` as an owned `Vec<f64>`, erroring
/// on multi-channel input.
fn mono_f64<T: StandardSample>(audio: &AudioSamples<'_, T>) -> AudioSampleResult<Vec<f64>> {
    if audio.num_channels().get() != 1 {
        return Err(AudioSampleError::Layout(
            LayoutError::channel_count_unsupported(
                "convolve/deconvolve (mono only)",
                ChannelRequirement::Mono,
                audio.num_channels().get(),
            ),
        ));
    }
    let working = audio.to_format::<f64>();
    let ch = working.as_slice().expect("validated single channel");
    Ok(ch.to_vec())
}

/// Runs an FFT convolution-style op on the mono channels of `a` and `b`.
fn convolve_or_deconvolve_mono<T, F>(
    a: &AudioSamples<'_, T>,
    b: &AudioSamples<'_, T>,
    op: F,
) -> AudioSampleResult<Vec<f64>>
where
    T: StandardSample,
    F: FnOnce(
        &non_empty_slice::NonEmptySlice<f64>,
        &non_empty_slice::NonEmptySlice<f64>,
    ) -> AudioSampleResult<non_empty_slice::NonEmptyVec<f64>>,
{
    let a_vec = mono_f64(a)?;
    let b_vec = mono_f64(b)?;
    // AudioSamples::new_mono rejects empty data, so this error path is unreachable for well-formed inputs.
    let a_ne = non_empty_slice::NonEmptySlice::new(a_vec.as_slice()).ok_or_else(|| {
        AudioSampleError::Parameter(ParameterError::invalid_value("self", "empty signal"))
    })?;
    let b_ne = non_empty_slice::NonEmptySlice::new(b_vec.as_slice()).ok_or_else(|| {
        AudioSampleError::Parameter(ParameterError::invalid_value("other", "empty signal"))
    })?;
    Ok(op(a_ne, b_ne)?.into_vec())
}

/// Wraps a `Vec<f64>` back into a mono `AudioSamples<T>` at `sample_rate`.
fn from_f64_mono<T: StandardSample>(
    data: Vec<f64>,
    sample_rate: std::num::NonZeroU32,
) -> AudioSampleResult<AudioSamples<'static, T>> {
    use non_empty_slice::NonEmptyVec;
    let ne_vec = NonEmptyVec::new(data).map_err(|_| {
        AudioSampleError::Parameter(ParameterError::invalid_value("data", "empty result"))
    })?;
    Ok(AudioSamples::<T>::from_mono_vec::<f64>(ne_vec, sample_rate))
}

#[cfg(test)]
mod convolution_tests {
    use crate::{AudioSampleError, AudioSamples, AudioTransforms, LayoutError, sample_rate};
    use ndarray::array;

    #[test]
    fn convolve_with_unit_impulse_is_identity() {
        let signal =
            AudioSamples::new_mono(array![1.0_f64, 2.0, 3.0], sample_rate!(48000)).unwrap();
        let impulse = AudioSamples::new_mono(array![1.0_f64], sample_rate!(48000)).unwrap();
        let out = signal.convolve(&impulse).unwrap();
        assert_eq!(out.sample_rate().get(), 48000);
        let ch = out.as_slice().unwrap();
        assert!((ch[0] - 1.0).abs() < 1e-9);
        assert!((ch[1] - 2.0).abs() < 1e-9);
        assert!((ch[2] - 3.0).abs() < 1e-9);
    }

    #[test]
    fn deconvolve_inverts_convolve() {
        let excitation = AudioSamples::new_mono(
            array![1.0_f64, 0.7, -0.3, 0.2, 0.9, -0.5, 0.1, 0.4],
            sample_rate!(48000),
        )
        .unwrap();
        let system =
            AudioSamples::new_mono(array![0.0_f64, 0.0, 1.0, 0.5], sample_rate!(48000)).unwrap();
        let recorded = excitation.convolve(&system).unwrap();
        let recovered = recorded.deconvolve(&excitation, 0.0).unwrap();
        let r = recovered.as_slice().unwrap();
        for (i, want) in [0.0, 0.0, 1.0, 0.5].iter().enumerate() {
            assert!(
                (r[i] - want).abs() < 1e-6,
                "tap {i}: got {}, want {want}",
                r[i]
            );
        }
    }

    #[test]
    fn multichannel_input_is_rejected() {
        // 2 channels, 4 samples each
        let stereo = AudioSamples::new_multi_channel(
            array![[1.0_f64, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0]],
            sample_rate!(48000),
        )
        .unwrap();
        let mono = AudioSamples::new_mono(array![1.0_f64], sample_rate!(48000)).unwrap();
        let err = stereo.convolve(&mono).unwrap_err();
        assert!(
            matches!(
                err,
                AudioSampleError::Layout(LayoutError::ChannelCountUnsupported { .. })
            ),
            "expected Layout(ChannelCountUnsupported), got: {err:?}"
        );
    }
}

#[cfg(test)]
mod psd_tests {
    use super::Psd;

    #[test]
    fn psd_accessors_and_into_parts_round_trip() {
        let psd = Psd::new(vec![0.0, 10.0, 20.0], vec![1.0, 0.5, 0.25]);
        assert_eq!(psd.frequencies(), &[0.0, 10.0, 20.0]);
        assert_eq!(psd.density(), &[1.0, 0.5, 0.25]);

        let (freqs, density) = psd.into_parts();
        assert_eq!(freqs, vec![0.0, 10.0, 20.0]);
        assert_eq!(density, vec![1.0, 0.5, 0.25]);
    }
}