rosa 0.1.0

//! Inverse CQT reconstruction.

use crate::matrix::Matrix;

/// ICQT parameters.
pub struct IcqtParams {
    /// Sample rate in Hz.
    pub sr: f64,
    /// Number of samples between successive frames.
    pub hop_length: usize,
    /// Minimum frequency in Hz.
    pub fmin: f64,
    /// Number of bins per octave.
    pub bins_per_octave: usize,
    /// Tuning offset in fractions of a bin.
    pub tuning: f64,
    /// Scaling factor applied to filter lengths.
    pub filter_scale: f64,
    /// None = no norm, Some(1.0) = L1 (default), Some(2.0) = L2.
    pub norm: Option<f64>,
    /// Window function used to shape each filter.
    pub window: fn(usize, bool) -> Vec<f64>,
    /// Whether bins were scaled in the forward CQT.
    pub scale: bool,
    /// If Some, trim or pad the output to this number of samples.
    pub length: Option<usize>,
}

impl Default for IcqtParams {
    fn default() -> Self {
        Self {
            sr: 22050.0,
            hop_length: 512,
            fmin: 32.703_195_662_574_83,
            bins_per_octave: 12,
            tuning: 0.0,
            filter_scale: 1.0,
            norm: Some(1.0),
            window: crate::windows::hann,
            scale: true,
            length: None,
        }
    }
}

/// Inverse CQT — reconstruct signal from CQT coefficients.
///
/// Matches `librosa.icqt(C, sr=..., ...)`.
/// Returns `Vec<f64>` of reconstructed audio.
///
/// Note: This is a simplified implementation using pseudo-inverse approach.
/// Full octave-decimation ICQT can be added later.
pub fn icqt(cqt_re: &Matrix, cqt_im: &Matrix, params: &IcqtParams) -> Vec<f64> {
    let n_bins = cqt_re.rows();
    let n_frames = cqt_re.cols();

    if n_frames == 0 {
        return Vec::new();
    }

    // Build wavelet filters for all bins
    let freqs =
        crate::convert::cqt_frequencies(n_bins, params.fmin, params.bins_per_octave, params.tuning);

    let wl_params = super::wavelet::WaveletLengthsParams {
        sr: params.sr,
        window: params.window,
        filter_scale: params.filter_scale,
        gamma: Some(0.0),
    };
    let (lengths, _) = super::wavelet::wavelet_lengths(&freqs, &wl_params);

    // Use single-FFT approach (like pseudo_cqt inverse)
    let wp = super::wavelet::WaveletParams {
        sr: params.sr,
        window: params.window,
        filter_scale: params.filter_scale,
        pad_fft: true,
        norm: params.norm,
        gamma: Some(0.0),
    };
    let (filters_re, filters_im, _) = super::wavelet::wavelet(&freqs, &wp);

    let min_fft = (2 * params.hop_length).next_power_of_two();
    let n_fft = filters_re.cols().max(min_fft).next_power_of_two();
    let n_bins_fft = n_fft / 2 + 1;

    // FFT the filter basis
    let mut planner = realfft::RealFftPlanner::<f64>::new();
    let fft = planner.plan_fft_forward(n_fft);
    let mut input_buf = fft.make_input_vec();
    let mut output_buf = fft.make_output_vec();
    let mut scratch = fft.make_scratch_vec();

    let mut basis_re_data = vec![0.0; n_bins * n_bins_fft];
    let mut basis_im_data = vec![0.0; n_bins * n_bins_fft];

    for k in 0..n_bins {
        let scale = lengths[k] / n_fft as f64;
        let filter_cols = filters_re.cols();
        let pad_start = if n_fft > filter_cols {
            (n_fft - filter_cols) / 2
        } else {
            0
        };
        let copy_len = filter_cols.min(n_fft);

        // FFT real part
        input_buf.fill(0.0);
        for j in 0..copy_len {
            input_buf[pad_start + j] = filters_re.get(k, j) * scale;
        }
        fft.process_with_scratch(&mut input_buf, &mut output_buf, &mut scratch)
            .expect("FFT failed");

        let mut fft_re_re = vec![0.0; n_bins_fft];
        let mut fft_re_im = vec![0.0; n_bins_fft];
        for j in 0..n_bins_fft {
            fft_re_re[j] = output_buf[j].re;
            fft_re_im[j] = output_buf[j].im;
        }

        // FFT imag part
        input_buf.fill(0.0);
        for j in 0..copy_len {
            input_buf[pad_start + j] = filters_im.get(k, j) * scale;
        }
        fft.process_with_scratch(&mut input_buf, &mut output_buf, &mut scratch)
            .expect("FFT failed");

        for j in 0..n_bins_fft {
            let fft_im_re = output_buf[j].re;
            let fft_im_im = output_buf[j].im;
            basis_re_data[k * n_bins_fft + j] = fft_re_re[j] - fft_im_im;
            basis_im_data[k * n_bins_fft + j] = fft_re_im[j] + fft_im_re;
        }
    }

    // Compute pseudo-inverse: conj(basis).T
    // inv_re[j, k] = basis_re[k, j] (conjugate transpose: flip sign of im)
    let mut inv_re_data = vec![0.0; n_bins_fft * n_bins];
    let mut inv_im_data = vec![0.0; n_bins_fft * n_bins];
    for k in 0..n_bins {
        for j in 0..n_bins_fft {
            inv_re_data[j * n_bins + k] = basis_re_data[k * n_bins_fft + j];
            inv_im_data[j * n_bins + k] = -basis_im_data[k * n_bins_fft + j]; // conjugate
        }
    }
    let inv_re = Matrix::from_vec(inv_re_data, n_bins_fft, n_bins);
    let inv_im = Matrix::from_vec(inv_im_data, n_bins_fft, n_bins);

    // Undo scaling if it was applied in forward
    let mut cqt_re = cqt_re.clone();
    let mut cqt_im = cqt_im.clone();
    if params.scale {
        for k in 0..n_bins {
            let s = lengths[k].sqrt();
            for f in 0..n_frames {
                *cqt_re.get_mut(k, f) *= s;
                *cqt_im.get_mut(k, f) *= s;
            }
        }
    }

    // D = inv_basis @ C
    let (d_re, d_im) = super::transform::complex_matmul(&inv_re, &inv_im, &cqt_re, &cqt_im);

    // ISTFT with boxcar window
    let istft_params = crate::stft::IstftParams {
        n_fft,
        hop_length: params.hop_length,
        win_length: n_fft,
        center: true,
        length: params.length,
    };

    crate::stft::istft(&d_re, &d_im, &istft_params)
}