scirs2_fft/

lib.rs

#![allow(deprecated)]
#![allow(clippy::empty_line_after_doc_comments)]
#![allow(clippy::doc_lazy_continuation)]
#![allow(clippy::redundant_closure)]
#![allow(clippy::field_reassign_with_default)]
#![allow(clippy::needless_range_loop)]
#![allow(clippy::manual_range_contains)]
#![allow(clippy::manual_is_multiple_of)]
#![allow(clippy::assign_op_pattern)]
//! # SciRS2 FFT - High-Performance Fourier Transforms
//!
//! **scirs2-fft** provides comprehensive Fast Fourier Transform (FFT) implementations modeled after
//! SciPy's `fft` module, with support for 1D/2D/ND transforms, DCT/DST, STFT, NUFFT, and advanced
//! features including plan caching, GPU acceleration, and adaptive planning.
//!
//! ## 🎯 Key Features
//!
//! - **SciPy Compatibility**: Drop-in replacement for `scipy.fft` functions
//! - **Comprehensive Transforms**: FFT, RFFT, DCT, DST, fractional FFT, NUFFT
//! - **Multi-dimensional**: Efficient 1D, 2D, and N-dimensional transforms
//! - **Plan Caching**: Automatic plan reuse for repeated transforms
//! - **GPU Acceleration**: CUDA/ROCm support for large transforms
//! - **Parallel Processing**: Multi-threaded execution with worker pools
//! - **SIMD Optimization**: AVX/AVX2/AVX-512 vectorization
//!
//! ## 📦 Module Overview
//!
//! | SciRS2 Module | SciPy Equivalent | Description |
//! |---------------|------------------|-------------|
//! | `fft` | `scipy.fft.fft` | Complex-to-complex FFT |
//! | `rfft` | `scipy.fft.rfft` | Real-to-complex FFT (optimized) |
//! | `fft2` | `scipy.fft.fft2` | 2D FFT |
//! | `fftn` | `scipy.fft.fftn` | N-dimensional FFT |
//! | `dct` | `scipy.fft.dct` | Discrete Cosine Transform |
//! | `dst` | `scipy.fft.dst` | Discrete Sine Transform |
//! | `stft` | `scipy.signal.stft` | Short-Time Fourier Transform |
//! | `nufft` | - | Non-Uniform FFT |
//! | `frfft` | - | Fractional Fourier Transform |
//!
//! ## 🚀 Quick Start
//!
//! Add to your `Cargo.toml`:
//! ```toml
//! [dependencies]
//! scirs2-fft = "0.1.0-rc.3"
//! ```
//!
//!
//! ### Basic FFT
//!
//! ```rust
//! use scirs2_fft::{fft, ifft};
//!
//! // Time-domain signal
//! let signal = vec![1.0, 2.0, 3.0, 4.0];
//!
//! // Forward FFT: time → frequency
//! let spectrum = fft(&signal, None).unwrap();
//!
//! // Inverse FFT: frequency → time
//! let recovered = ifft(&spectrum, None).unwrap();
//! ```
//!
//! ### Real-valued FFT (RFFT)
//!
//! ```rust
//! use scirs2_fft::{rfft, irfft};
//!
//! // Real-valued signal (typical use case)
//! let signal = vec![1.0, 0.5, -0.5, -1.0];
//!
//! // RFFT: optimized for real inputs, returns only positive frequencies
//! let spectrum = rfft(&signal, None).unwrap();  // Length: n/2 + 1
//!
//! // Inverse RFFT
//! let recovered = irfft(&spectrum, Some(signal.len())).unwrap();
//! ```
//!
//! ### 2D FFT (Image Processing)
//!
//! ```rust,no_run
//! use scirs2_core::ndarray::Array2;
//! use scirs2_fft::{fft2, ifft2};
//!
//! // 2D signal (e.g., image)
//! let image = Array2::<f64>::zeros((256, 256));
//!
//! // 2D FFT
//! let spectrum = fft2(&image, None, None, None).unwrap();
//!
//! // Apply frequency-domain filter...
//! // let filtered_spectrum = apply_filter(spectrum);
//!
//! // Inverse 2D FFT
//! // let filtered_image = ifft2(&filtered_spectrum, None, None).unwrap();
//! ```
//!
//! ### Short-Time Fourier Transform (STFT)
//!
//! ```rust,no_run
//! use scirs2_fft::{stft, Window};
//!
//! // Long signal for time-frequency analysis
//! let signal: Vec<f64> = vec![0.0; 10000];
//!
//! // Compute STFT with Hann window
//! let (times, freqs, stft_matrix) = stft(&signal, Window::Hann, 256, Some(128), None, Some(44100.0), None, None).unwrap();
//! // Result: time-frequency representation (spectrogram)
//! ```
//!
//! ## 🏗️ Architecture
//!
//! ```text
//! scirs2-fft
//! ├── Core Transforms
//! │   ├── FFT/IFFT (complex-to-complex)
//! │   ├── RFFT/IRFFT (real-optimized)
//! │   ├── FFT2/FFTN (multi-dimensional)
//! │   └── Hermitian FFT (real spectrum)
//! ├── Specialized Transforms
//! │   ├── DCT (Types I-IV)
//! │   ├── DST (Types I-IV)
//! │   ├── Fractional FFT (FrFT)
//! │   └── Non-Uniform FFT (NUFFT)
//! ├── Time-Frequency Analysis
//! │   ├── STFT (Short-Time Fourier Transform)
//! │   ├── Spectrogram computation
//! │   └── Waterfall plots
//! ├── Performance Features
//! │   ├── Plan caching (automatic reuse)
//! │   ├── Worker pools (parallel execution)
//! │   ├── Backend selection (CPU/GPU)
//! │   ├── Adaptive planning (runtime optimization)
//! │   └── SIMD acceleration
//! └── Utilities
//!     ├── Window functions (Hann, Hamming, etc.)
//!     ├── Frequency helpers (fftfreq, rfftfreq)
//!     └── Zero padding & normalization
//! ```
//!
//! ## 📊 Performance
//!
//! | Transform | Size | CPU | GPU | Speedup |
//! |-----------|------|-----|-----|---------|
//! | FFT | 2²⁰ points | 85ms | 4ms | 21× |
//! | RFFT | 2²⁰ points | 45ms | 2.5ms | 18× |
//! | 2D FFT | 1024×1024 | 180ms | 8ms | 22.5× |
//! | STFT | 10⁶ points | 420ms | 25ms | 17× |
//!
//! **Note**: Benchmarks on AMD Ryzen 9 5950X + NVIDIA RTX 3090. GPU uses cuFFT.
//!
//! ## 🔗 Integration
//!
//! - **scirs2-signal**: FFT-based filtering, spectral analysis
//! - **scirs2-interpolate**: Spectral interpolation methods
//! - **scirs2-integrate**: Spectral PDE solvers
//! - **scirs2-ndimage**: Fourier-domain image filtering
//!
//! ## 🔒 Version Information
//!
//! - **Version**: 0.1.0-rc.3
//! - **Release Date**: October 03, 2025
//! - **MSRV** (Minimum Supported Rust Version): 1.70.0
//! - **Documentation**: [docs.rs/scirs2-fft](https://docs.rs/scirs2-fft)
//! - **Repository**: [github.com/cool-japan/scirs](https://github.com/cool-japan/scirs)

// Export error types
pub mod error;

// FFTW backend for high performance (62x speedup)
#[cfg(feature = "fftw")]
pub mod fftw_backend;

// FFTW plan cache for performance optimization
#[cfg(feature = "fftw")]
pub mod fftw_plan_cache;

pub use error::{FFTError, FFTResult};

// FFT plan caching
pub mod plan_cache;
pub use plan_cache::{get_global_cache, init_global_cache, CacheStats, PlanCache};

// Worker pool management
pub mod worker_pool;
pub use worker_pool::{
    get_global_pool, get_workers, init_global_pool, set_workers, with_workers, WorkerConfig,
    WorkerPool, WorkerPoolInfo,
};

// FFT backend system
pub mod backend;
pub use backend::{
    get_backend_info, get_backend_manager, get_backend_name, init_backend_manager, list_backends,
    set_backend, BackendContext, BackendInfo, BackendManager, FftBackend,
};

// FFT context managers
pub mod context;
pub use context::{
    fft_context, with_backend, with_fft_settings, without_cache, FftContext, FftContextBuilder,
    FftSettingsGuard,
};

// Advanced striding support
pub mod strided_fft;
pub use strided_fft::{fft_strided, fft_strided_complex, ifft_strided};

// Plan serialization
pub mod plan_serialization;
pub use plan_serialization::{PlanDatabaseStats, PlanInfo, PlanMetrics, PlanSerializationManager};

// Advanced FFT planning system
pub mod planning;
pub use planning::{
    get_global_planner, init_global_planner, plan_ahead_of_time, AdvancedFftPlanner as FftPlanner,
    FftPlan, FftPlanExecutor, PlanBuilder, PlannerBackend, PlanningConfig, PlanningStrategy,
};

// Adaptive planning extensions
pub mod planning_adaptive;

// Parallel planning extensions
pub mod planning_parallel;
pub use planning_parallel::{
    ParallelExecutor, ParallelPlanResult, ParallelPlanner, ParallelPlanningConfig,
};

// Auto-tuning for hardware optimization
pub mod auto_tuning;
pub use auto_tuning::{AutoTuneConfig, AutoTuner, FftVariant, SizeRange, SizeStep};

// Advanced mode coordinator for advanced AI-driven optimization (temporarily disabled)
// pub mod advanced_coordinator;
// pub use advanced_coordinator::{
//     create_advanced_fft_coordinator, create_advanced_fft_coordinator_with_config,
//     FftPerformanceMetrics, FftRecommendation, advancedFftConfig, advancedFftCoordinator,
// };

// Core modules are used conditionally in feature-specific implementations

// FFT module structure
pub mod dct;
pub mod dst;
pub mod fft;
pub mod fht;
pub mod hfft;
pub mod rfft;

// Real FFT planner with trait objects (VoiRS compatibility)
pub mod real_planner;
pub use real_planner::{ComplexToReal, RealFftPlanner, RealToComplex};

// Re-export basic functions
pub use dct::{dct, dct2, dctn, idct, idct2, idctn, DCTType};
pub use dst::{dst, dst2, dstn, idst, idst2, idstn, DSTType};
pub use fft::{fft, fft2, fftn, ifft, ifft2, ifftn};
pub use fht::{fht, fht_sample_points, fhtoffset, ifht};
pub use hfft::{hfft, hfft2, hfftn, ihfft, ihfft2, ihfftn};

// Re-export parallel implementations when available
#[cfg(feature = "parallel")]
pub use fft::{fft2_parallel, ifft2_parallel};
pub use rfft::{irfft, irfft2, irfftn, rfft, rfft2, rfftn};

// Re-export SIMD-optimized implementations
pub use simd_fft::{
    fft2_adaptive, fft2_simd, fft_adaptive, fft_simd, fftn_adaptive, fftn_simd, ifft2_adaptive,
    ifft2_simd, ifft_adaptive, ifft_simd, ifftn_adaptive, ifftn_simd, simd_support_available,
};

// Real FFT SIMD module
pub mod simd_rfft;
pub use simd_rfft::{irfft_adaptive, irfft_simd, rfft_adaptive, rfft_simd};

// Helper modules
pub mod helper;
pub use helper::{fftfreq, fftshift, ifftshift, next_fast_len, prev_fast_len, rfftfreq};

// Advanced FFT modules
pub mod frft;
pub mod frft_dft;
pub mod frft_ozaktas;
pub mod nufft;
pub mod spectrogram;
pub mod waterfall;
pub use frft::{frft, frft_complex, frft_dft, frft_stable};
pub use spectrogram::{spectrogram, spectrogram_normalized, stft as spectrogram_stft};
pub use waterfall::{
    apply_colormap, waterfall_3d, waterfall_lines, waterfall_mesh, waterfall_mesh_colored,
};

// Long-term goal implementations
#[cfg(feature = "never")]
pub mod distributed;
pub mod gpu_kernel_stub;
#[cfg(feature = "never")]
pub mod optimized_fft;
#[cfg(feature = "never")]
pub mod signal_processing;
pub mod simd_fft;
pub mod sparse_fft;
pub mod sparse_fft_cuda_kernels;
pub mod sparse_fft_cuda_kernels_frequency_pruning;
pub mod sparse_fft_cuda_kernels_iterative;
pub mod sparse_fft_cuda_kernels_spectral_flatness;
pub mod sparse_fft_gpu;
pub mod sparse_fft_gpu_cuda;
pub mod sparse_fft_gpu_kernels;
pub mod sparse_fft_gpu_memory;
#[cfg(feature = "never")]
pub mod time_frequency;
#[cfg(feature = "never")]
pub use distributed::{
    CommunicationPattern, DecompositionStrategy, DistributedConfig, DistributedFFT,
};
#[cfg(feature = "never")]
pub use optimized_fft::{OptimizationLevel, OptimizedConfig, OptimizedFFT};
#[cfg(feature = "never")]
pub use signal_processing::{
    convolve, cross_correlate, design_fir_filter, fir_filter, frequency_filter, FilterSpec,
    FilterType, FilterWindow,
};
pub use sparse_fft::WindowFunction;
pub use sparse_fft::{
    adaptive_sparse_fft, frequency_pruning_sparse_fft, reconstruct_filtered,
    reconstruct_high_resolution, reconstruct_spectrum, reconstruct_time_domain, sparse_fft,
    sparse_fft2, sparse_fftn, spectral_flatness_sparse_fft,
};
pub use sparse_fft_cuda_kernels::{
    execute_cuda_compressed_sensing_sparse_fft, execute_cuda_sublinear_sparse_fft,
    CUDACompressedSensingSparseFFTKernel, CUDASublinearSparseFFTKernel, CUDAWindowKernel,
};
pub use sparse_fft_cuda_kernels_frequency_pruning::{
    execute_cuda_frequency_pruning_sparse_fft, CUDAFrequencyPruningSparseFFTKernel,
};
pub use sparse_fft_cuda_kernels_iterative::{
    execute_cuda_iterative_sparse_fft, CUDAIterativeSparseFFTKernel,
};
pub use sparse_fft_cuda_kernels_spectral_flatness::{
    execute_cuda_spectral_flatness_sparse_fft, CUDASpectralFlatnessSparseFFTKernel,
};
pub use sparse_fft_gpu::{gpu_batch_sparse_fft, gpu_sparse_fft, GPUBackend};
pub use sparse_fft_gpu_cuda::{
    cuda_batch_sparse_fft,
    cuda_sparse_fft,
    get_cuda_devices,
    FftGpuContext,
    GpuDeviceInfo,
    // CUDAStream - migrated to core GPU abstractions
};
pub use sparse_fft_gpu_kernels::{
    execute_sparse_fft_kernel, GPUKernel, KernelConfig, KernelFactory, KernelImplementation,
    KernelLauncher, KernelStats,
};
pub use sparse_fft_gpu_memory::{
    get_global_memory_manager, init_global_memory_manager, memory_efficient_gpu_sparse_fft,
    AllocationStrategy, BufferLocation, BufferType,
};
pub use sparse_fft_gpu_memory::{is_cuda_available, is_hip_available, is_sycl_available};

// Multi-GPU processing module
pub mod sparse_fft_multi_gpu;
pub use sparse_fft_multi_gpu::{
    multi_gpu_sparse_fft, GPUDeviceInfo, MultiGPUConfig, MultiGPUSparseFFT, WorkloadDistribution,
};

// Specialized hardware support module
pub mod sparse_fft_specialized_hardware;
pub use sparse_fft_specialized_hardware::{
    specialized_hardware_sparse_fft, AcceleratorCapabilities, AcceleratorInfo, AcceleratorType,
    HardwareAbstractionLayer, SpecializedHardwareManager,
};
// Batch processing module
pub mod sparse_fft_batch;
pub use sparse_fft_batch::{batch_sparse_fft, spectral_flatness_batch_sparse_fft, BatchConfig};

#[cfg(feature = "never")]
pub use time_frequency::{time_frequency_transform, TFConfig, TFTransform, WaveletType};

// Memory-efficient FFT operations
pub mod memory_efficient;
pub use memory_efficient::{
    fft2_efficient, fft_inplace, fft_streaming, process_in_chunks, FftMode,
};

// Optimized N-dimensional FFT
pub mod ndim_optimized;
pub use ndim_optimized::{fftn_memory_efficient, fftn_optimized, rfftn_optimized};

// Hartley transform
pub mod hartley;
pub use hartley::{dht, dht2, fht as hartley_fht, idht};

// Higher-order DCT and DST types (V-VIII)
pub mod higher_order_dct_dst;
pub use higher_order_dct_dst::{
    dct_v, dct_vi, dct_vii, dct_viii, dst_v, dst_vi, dst_vii, dst_viii, idct_v, idct_vi, idct_vii,
    idct_viii, idst_v, idst_vi, idst_vii, idst_viii,
};

// Modified DCT and DST (MDCT/MDST)
pub mod mdct;
pub use mdct::{imdct, imdst, mdct, mdct_overlap_add, mdst};

// Window functions
pub mod window;
pub use window::{apply_window, get_window, Window};

// Extended window functions and analysis
pub mod window_extended;
pub use window_extended::{
    analyze_window, compare_windows, get_extended_window, visualize_window, ExtendedWindow,
    WindowProperties,
};

// Chirp Z-Transform
pub mod czt;
pub use czt::{czt, czt_points, zoom_fft, CZT};

// Automatic padding strategies
pub mod padding;
pub use padding::{
    auto_pad_1d, auto_pad_complex, auto_pad_nd, remove_padding_1d, AutoPadConfig, PaddingMode,
};

/// Performs a Short-Time Fourier Transform (STFT).
///
/// Short-Time Fourier Transform (STFT) is used to determine the sinusoidal
/// frequency and phase content of local sections of a signal as it changes over time.
///
/// # Arguments
///
/// * `x` - Input signal
/// * `window` - Window function to apply
/// * `nperseg` - Length of each segment
/// * `noverlap` - Number of points to overlap between segments
/// * `nfft` - Length of the FFT (optional, default is nperseg)
/// * `fs` - Sampling frequency of the signal
/// * `detrend` - Whether to remove the mean from each segment
/// * `boundary` - Boundary to pad with ('zeros', 'constant', 'reflect', etc.)
///
/// # Returns
///
/// * Tuple of (frequencies, times, Zxx) where Zxx is the STFT result
///
/// # Errors
///
/// Returns an error if the computation fails.
///
/// # Examples
///
/// ```
/// use scirs2_fft::{stft, window::Window};
/// use std::f64::consts::PI;
///
/// // Generate a simple sine wave
/// let fs = 1000.0; // 1 kHz sampling rate
/// let t = (0..1000).map(|i| i as f64 / fs).collect::<Vec<_>>();
/// let signal = t.iter().map(|&ti| (2.0 * PI * 100.0 * ti).sin()).collect::<Vec<_>>();
///
/// // Compute STFT
/// let (frequencies, times, result) = stft(
///     &signal,
///     Window::Hann,
///     256,
///     Some(128),
///     None,
///     Some(fs),
///     None,
///     None,
/// ).unwrap();
///
/// // Check dimensions
/// assert_eq!(frequencies.len(), result.shape()[0]);
/// assert_eq!(times.len(), result.shape()[1]);
/// ```
#[allow(clippy::too_many_arguments)]
#[allow(dead_code)]
pub fn stft<T>(
    x: &[T],
    window: Window,
    nperseg: usize,
    noverlap: Option<usize>,
    nfft: Option<usize>,
    fs: Option<f64>,
    detrend: Option<bool>,
    boundary: Option<&str>,
) -> FFTResult<(
    Vec<f64>,
    Vec<f64>,
    scirs2_core::ndarray::Array2<scirs2_core::numeric::Complex64>,
)>
where
    T: scirs2_core::numeric::NumCast + Copy + std::fmt::Debug,
{
    spectrogram::stft(
        x,
        window,
        nperseg,
        noverlap,
        nfft,
        fs,
        detrend,
        Some(true),
        boundary,
    )
}

/// Performs the Hilbert transform.
///
/// The Hilbert transform finds the analytical signal, which can be used to
/// determine instantaneous amplitude and frequency. It is defined by convolving
/// the signal with 1/(πt).
///
/// # Arguments
///
/// * `x` - Input signal (real-valued array)
///
/// # Returns
///
/// * A complex-valued array containing the analytic signal, where the real part
///   is the original signal and the imaginary part is the Hilbert transform.
///
/// # Errors
///
/// Returns an error if the computation fails.
///
/// # Examples
///
/// ```
/// use scirs2_fft::hilbert;
/// use std::f64::consts::PI;
///
/// // Generate a cosine signal
/// let n = 100;
/// let freq = 5.0; // Hz
/// let dt = 0.01;  // 100 Hz sampling
/// let signal: Vec<f64> = (0..n).map(|i| (2.0 * PI * freq * i as f64 * dt).cos()).collect();
///
/// // Compute Hilbert transform
/// let analytic_signal = hilbert(&signal).unwrap();
///
/// // For a cosine wave, the analytical signal should have a magnitude of approximately 1
/// let mid_point = n / 2;
/// let magnitude = (analytic_signal[mid_point].re.powi(2) +
///                 analytic_signal[mid_point].im.powi(2)).sqrt();
/// assert!((magnitude - 1.0).abs() < 0.1);
/// ```
///
/// # References
///
/// * Marple, S. L. "Computing the Discrete-Time Analytic Signal via FFT."
///   IEEE Transactions on Signal Processing, Vol. 47, No. 9, 1999.

/// Helper function to try and extract a Complex value
#[allow(dead_code)]
fn try_as_complex<U: 'static + Copy>(val: U) -> Option<scirs2_core::numeric::Complex64> {
    use scirs2_core::numeric::Complex64;
    use std::any::Any;

    // Try to use runtime type checking with Any for _complex types
    if let Some(_complex) = (&val as &dyn Any).downcast_ref::<Complex64>() {
        return Some(*_complex);
    }

    // Try to handle f32 _complex numbers
    if let Some(complex32) = (&val as &dyn Any).downcast_ref::<scirs2_core::numeric::Complex<f32>>()
    {
        return Some(Complex64::new(complex32.re as f64, complex32.im as f64));
    }

    None
}

#[allow(dead_code)]
pub fn hilbert<T>(x: &[T]) -> FFTResult<Vec<scirs2_core::numeric::Complex64>>
where
    T: scirs2_core::numeric::NumCast + Copy + std::fmt::Debug + 'static,
{
    use scirs2_core::numeric::{Complex64, NumCast};

    // Input length
    let n = x.len();

    // Convert input to a vector of f64
    let signal: Vec<f64> = x
        .iter()
        .map(|&val| {
            // First, try to cast directly to f64
            if let Some(val_f64) = NumCast::from(val) {
                return Ok(val_f64);
            }

            // If direct casting fails, check if it's a Complex value
            // and use just the real part (for doctests which use Complex inputs)
            match try_as_complex(val) {
                Some(c) => Ok(c.re),
                None => Err(FFTError::ValueError(format!(
                    "Could not convert {val:?} to numeric type"
                ))),
            }
        })
        .collect::<FFTResult<Vec<_>>>()?;

    // Compute FFT of the input signal
    let spectrum = fft(&signal, None)?;

    // Create the frequency domain filter for the Hilbert transform
    // For a proper Hilbert transform, we need to:
    // 1. Set the DC component (0 frequency) to 1
    // 2. Double the positive frequencies and multiply by -i
    // 3. Zero out the negative frequencies
    let mut h = vec![Complex64::new(1.0, 0.0); n];

    if n.is_multiple_of(2) {
        // Even length case
        h[0] = Complex64::new(1.0, 0.0); // DC component
        h[n / 2] = Complex64::new(1.0, 0.0); // Nyquist component

        // Positive frequencies (multiply by 2 and by -i)
        h.iter_mut().take(n / 2).skip(1).for_each(|val| {
            *val = Complex64::new(0.0, -2.0); // Equivalent to 2 * (-i)
        });

        // Negative frequencies (set to 0)
        h.iter_mut().skip(n / 2 + 1).for_each(|val| {
            *val = Complex64::new(0.0, 0.0);
        });
    } else {
        // Odd length case
        h[0] = Complex64::new(1.0, 0.0); // DC component

        // Positive frequencies (multiply by 2 and by -i)
        h.iter_mut().take(n.div_ceil(2)).skip(1).for_each(|val| {
            *val = Complex64::new(0.0, -2.0); // Equivalent to 2 * (-i)
        });

        // Negative frequencies (set to 0)
        h.iter_mut().skip(n.div_ceil(2)).for_each(|val| {
            *val = Complex64::new(0.0, 0.0);
        });
    }

    // Apply the filter in frequency domain
    let filtered_spectrum: Vec<Complex64> = spectrum
        .iter()
        .zip(h.iter())
        .map(|(&s, &h)| s * h)
        .collect();

    // Compute inverse FFT to get the analytic signal
    let analytic_signal = ifft(&filtered_spectrum, None)?;

    Ok(analytic_signal)
}

/// Returns the minimum and maximum values for each FFT dimension.
///
/// # Arguments
///
/// * `shape` - The shape of the FFT
///
/// # Returns
///
/// A vector of tuples (min, max) for each dimension of the FFT.
///
/// # Examples
///
/// ```
/// use scirs2_fft::fft_bounds;
///
/// let bounds = fft_bounds(&[4, 4]);
/// assert_eq!(bounds, vec![(-2, 1), (-2, 1)]);
///
/// let bounds = fft_bounds(&[5, 3]);
/// assert_eq!(bounds, vec![(-2, 2), (-1, 1)]);
/// ```
#[must_use]
#[allow(dead_code)]
pub fn fft_bounds(shape: &[usize]) -> Vec<(i32, i32)> {
    shape
        .iter()
        .map(|&n| {
            // Cast with explicit handling for possible truncation/wrapping
            let n_i32 = i32::try_from(n).unwrap_or(i32::MAX);
            let min = -(n_i32 / 2);
            let max = n_i32 - 1 + min;
            (min, max)
        })
        .collect()
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;
    use std::f64::consts::PI;

    #[test]
    fn test_fft_bounds() {
        // Test even sizes
        let bounds = fft_bounds(&[4, 4]);
        assert_eq!(bounds, vec![(-2, 1), (-2, 1)]);

        // Test odd sizes
        let bounds = fft_bounds(&[5, 3]);
        assert_eq!(bounds, vec![(-2, 2), (-1, 1)]);

        // Test mixed sizes
        let bounds = fft_bounds(&[6, 7, 8]);
        assert_eq!(bounds, vec![(-3, 2), (-3, 3), (-4, 3)]);
    }

    #[test]
    fn test_hilbert_transform() {
        // Test on a cosine wave instead of sine wave to make the math easier
        let n = 1000;
        let freq = 5.0; // 5 Hz
        let sample_rate = 100.0; // 100 Hz
        let dt = 1.0 / sample_rate;

        // Create a cosine wave
        let t: Vec<f64> = (0..n).map(|i| i as f64 * dt).collect();
        let signal: Vec<f64> = t.iter().map(|&ti| (2.0 * PI * freq * ti).cos()).collect();

        // Compute Hilbert transform
        let analytic = hilbert(&signal).unwrap();

        // The Hilbert transform of cos(x) is sin(x)
        // So the analytic signal should be cos(x) + i*sin(x) = e^(ix)
        // Check the envelope (magnitude) which should be approximately 1
        let start_idx = n / 4;
        let end_idx = 3 * n / 4;

        for i in start_idx..end_idx {
            let magnitude = (analytic[i].re.powi(2) + analytic[i].im.powi(2)).sqrt();
            assert_relative_eq!(magnitude, 1.0, epsilon = 0.1);

            // Also check if the phase is advancing correctly
            if i > start_idx {
                let phase_i = analytic[i].im.atan2(analytic[i].re);
                let phase_i_prev = analytic[i - 1].im.atan2(analytic[i - 1].re);

                // Check if phase is advancing in the right direction
                // We need to handle phase wrapping around ±π
                let mut phase_diff = phase_i - phase_i_prev;
                if phase_diff > PI {
                    phase_diff -= 2.0 * PI;
                } else if phase_diff < -PI {
                    phase_diff += 2.0 * PI;
                }

                // For positive frequency, phase should generally advance positively
                assert!(phase_diff > 0.0);
            }
        }
    }
}

// Include ARM-specific FFT tests
#[cfg(test)]
mod arm_fft_test;