Metadata-Version: 2.2
Name: sgnl-cpu-interp
Version: 0.2.0
Summary: Fast polyphase resampling with multi-architecture SIMD support
Keywords: signal processing,interpolation,upsampling,polyphase,resampling,SIMD
Author-Email: Chad Hanna <crh184@psu.edu>
License: MIT
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: C
Classifier: Topic :: Scientific/Engineering
Classifier: Operating System :: OS Independent
Project-URL: Homepage, https://git.ligo.org/greg/fast-resample-cpu
Project-URL: Documentation, https://greg.docs.ligo.org/fast-resample-cpu
Project-URL: Issues, https://git.ligo.org/greg/fast-resample-cpu/issues
Requires-Python: >=3.9
Requires-Dist: numpy>=1.19.0
Provides-Extra: bench
Requires-Dist: scipy; extra == "bench"
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Requires-Dist: pytest-cov; extra == "test"
Provides-Extra: lint
Requires-Dist: black; extra == "lint"
Requires-Dist: flake8; extra == "lint"
Requires-Dist: flake8-bandit; extra == "lint"
Requires-Dist: flake8-black; extra == "lint"
Requires-Dist: flake8-bugbear; extra == "lint"
Requires-Dist: flake8-isort; extra == "lint"
Requires-Dist: flake8-pyproject; extra == "lint"
Requires-Dist: isort; extra == "lint"
Requires-Dist: mypy; extra == "lint"
Requires-Dist: mypy-extensions; extra == "lint"
Requires-Dist: typing_extensions; extra == "lint"
Provides-Extra: dev
Requires-Dist: sgnl-cpu-interp[lint]; extra == "dev"
Requires-Dist: sgnl-cpu-interp[test]; extra == "dev"
Description-Content-Type: text/markdown

<h1 align="center">sgnl-cpu-interp</h1>

<p align="center">Fast polyphase resampling for multichannel data with multi-architecture SIMD support</p>

<p align="center">
  <a href="https://git.ligo.org/greg/fast-resample-cpu/-/pipelines/latest">
    <img alt="pipeline" src="https://git.ligo.org/greg/fast-resample-cpu/badges/main/pipeline.svg" />
  </a>
  <a href="https://git.ligo.org/greg/fast-resample-cpu/-/pipelines/latest">
    <img alt="coverage" src="https://git.ligo.org/greg/fast-resample-cpu/badges/main/coverage.svg" />
  </a>
  <a href="https://pypi.org/project/sgnl-cpu-interp/">
    <img alt="pypi version" src="https://img.shields.io/pypi/v/sgnl-cpu-interp.svg" />
  </a>
</p>

---

## Features

- **Multi-Architecture SIMD**: Automatic runtime CPU detection with optimized kernels for:
  - x86_64: AVX-512, AVX2+FMA, AVX, SSE4.1, SSE2
  - ARM64: NEON (Apple Silicon, AWS Graviton, etc.)
  - Fallback: Optimized scalar implementation
- **No External Dependencies**: Only requires NumPy (removed GSL and FFTW dependencies)
- **High Performance**: ~5x faster than GSL-based implementations
- **Multichannel**: Optimized for processing many channels simultaneously (tested with 1024+ channels)
- **Quality**: Lanczos-windowed sinc interpolation for high-quality upsampling
- **Two Memory Layouts**: Supports both (time, channels) and (channels, time) layouts
- **Simple API**: Easy-to-use NumPy-based interface

## Installation

### From PyPI (when available)

```bash
pip install sgnl-cpu-interp
```

Binary wheels are built in CI for:

| Platform | Architecture | Minimum version |
|----------|--------------|-----------------|
| Linux | x86_64, aarch64 | glibc 2.28 / `manylinux_2_28` (e.g. RHEL/Rocky 8, Debian 10, Ubuntu 18.10) |
| macOS | x86_64 | macOS 13 (Ventura) |
| macOS | arm64 | macOS 15 (Sequoia) |

### From source

```bash
git clone https://git.ligo.org/greg/fast-resample-cpu.git
cd fast-resample-cpu
pip install .
```

No external dependencies are required beyond NumPy and a C compiler; the
build backend (scikit-build-core) provisions CMake and Ninja automatically
if they are not already installed. The build detects your CPU architecture
and compiles the appropriate SIMD kernels; the best implementation is
selected at runtime.

## Quick Start

```python
import numpy as np
from sgnl_cpu_interp import upsample, get_simd_info

# Check which SIMD implementation is being used
print(get_simd_info())
# {'implementation': 'NEON', 'available': ['NEON', 'Scalar'], 'cpu_features': 'NEON+FMA'}

# Upsample a 50 Hz sine wave from 128 Hz to 2048 Hz
fs_in = 128
fs_out = 2048
factor = fs_out // fs_in  # 16x upsampling

# Generate test signal
t = np.arange(0, 0.5, 1/fs_in)
signal = np.sin(2 * np.pi * 50 * t).astype(np.float32)

# Upsample
upsampled = upsample(signal, factor=factor, half_length=8)
print(f"Input: {len(signal)} samples at {fs_in} Hz")
print(f"Output: {len(upsampled)} samples at {fs_out} Hz")
```

## Usage Examples

### Single channel upsampling
```python
import numpy as np
from sgnl_cpu_interp import upsample

# 1D signal (single channel)
signal = np.random.randn(1024).astype(np.float32)
upsampled = upsample(signal, factor=2)
```

### Multichannel upsampling
```python
# 2D array: (n_samples, n_channels)
n_samples, n_channels = 1024, 128
data = np.random.randn(n_samples, n_channels).astype(np.float32)

# Upsample by factor of 2
upsampled = upsample(data, factor=2)
print(upsampled.shape)  # (2016, 128) - note: loses 2*half_length samples

# Upsample by factor of 4 with longer kernel for better quality
upsampled = upsample(data, factor=4, half_length=16)
print(upsampled.shape)  # (3972, 128)
```

### Transposed layout
```python
from sgnl_cpu_interp import upsample_transposed

# Transposed layout: (n_channels, n_samples)
data = np.random.randn(128, 1024).astype(np.float32)
upsampled = upsample_transposed(data, factor=2)
print(upsampled.shape)  # (128, 2016)
```

## API Reference

### `upsample(data, factor=2, half_length=8)`

Upsample multichannel data using polyphase filtering (standard layout).

**Parameters:**
- `data` (ndarray): Input array of shape `(n_samples,)` for single channel or `(n_samples, n_channels)` for multichannel. float32, float64, complex64, and complex128 are processed natively and preserved in the output; any other dtype is converted to float32. Complex data is resampled by filtering the real and imaginary parts independently with the same real kernel.
- `factor` (int, optional): Upsampling factor (default: 2). Must be >= 2.
- `half_length` (int, optional): Half-length of the sinc kernel (default: 8). Larger values provide better quality but are slower. Total kernel length = `2 * half_length + 1`.

**Returns:**
- `output` (ndarray): Upsampled array of shape `((n_samples - kernel_len + 1) * factor,)` or `((n_samples - kernel_len + 1) * factor, n_channels)` where `kernel_len = 2 * half_length + 1`. Same dtype as the (possibly converted) input.

### `upsample_transposed(data, factor=2, half_length=8)`

Upsample multichannel data using polyphase filtering (transposed layout).

Same as `upsample()` but expects input in `(n_channels, n_samples)` layout. Use this when your data is already in channels-first format to avoid transpose overhead.

**Parameters:**
- `data` (ndarray): Input array of shape `(n_channels, n_samples)`. Must be 2D; float32, float64, complex64, and complex128 are processed natively and preserved, other dtypes are converted to float32.
- `factor` (int, optional): Upsampling factor (default: 2). Must be >= 2.
- `half_length` (int, optional): Half-length of the sinc kernel (default: 8).

**Returns:**
- `output` (ndarray): Upsampled array of shape `(n_channels, (n_samples - kernel_len + 1) * factor)`.

### `downsample(data, factor=2, half_length=8)`

Downsample multichannel data with an anti-aliasing FIR filter (standard layout).

The signal is filtered with a Lanczos-windowed sinc lowpass (cutoff at the output Nyquist, unit DC gain) and decimated in a single pass — only every `factor`-th output is computed. Any integer factor >= 2 is supported (e.g. `factor=1024` for 16384 Hz -> 16 Hz).

**Parameters:**
- `data` (ndarray): Input array of shape `(n_samples,)` or `(n_samples, n_channels)`. float32, float64, complex64, and complex128 are processed natively and preserved; other dtypes are converted to float32.
- `factor` (int, optional): Downsampling factor (default: 2). Must be >= 2.
- `half_length` (int, optional): Filter half-length parameter (default: 8). The filter has `2 * half_length * factor + 1` taps.

**Returns:**
- `output` (ndarray): Downsampled array of length `(n_samples - kernel_len) // factor + 1` per channel, where `kernel_len = 2 * half_length * factor + 1`. The group delay is exactly `half_length * factor` input samples.

### `downsample_transposed(data, factor=2, half_length=8)`

Same as `downsample()` but expects `(n_channels, n_samples)` layout.

### `get_simd_info()`

Get information about the current SIMD implementation.

**Returns:**
- `dict` with keys:
  - `implementation`: Name of current implementation (e.g., 'AVX2+FMA', 'NEON', 'Scalar')
  - `available`: List of all available implementations for this CPU
  - `cpu_features`: Detected CPU SIMD features

### `set_implementation(name)`

Manually select a SIMD implementation. Useful for testing and benchmarking.

**Parameters:**
- `name` (str): Implementation name from `get_simd_info()['available']`

Can also be set via the `SGNL_CPU_IMPL` environment variable:
```bash
SGNL_CPU_IMPL=Scalar python my_script.py
```

## Important Notes

- **Edge loss**: The convolution loses `kernel_len - 1` samples from the edges. For `half_length=8`, you lose 16 input samples.
- **Time alignment**: The output has a delay of `(kernel_len - 1) / 2` samples at the input sample rate.
- **Minimum length**: Input must have at least `kernel_len` samples.
- Uses Lanczos-windowed sinc kernel: `h(x) = sinc(x/factor) * sinc(x/kernel_length)`

## Performance

Benchmark on 1024 channels, 1024 samples:

| Platform | Implementation | Time | Notes |
|----------|---------------|------|-------|
| Apple Silicon (M-series) | NEON | ~0.4 ms | Auto-selected |
| Apple Silicon | Scalar | ~0.4 ms | Compiler auto-vectorizes well |
| x86_64 (Haswell+) | AVX2+FMA | ~0.3 ms | Expected |
| x86_64 (older) | SSE2 | ~0.8 ms | Baseline x86_64 |

Comparison with previous GSL-based implementation:

| Implementation | Time | Speedup |
|----------------|------|---------|
| GSL BLAS (old) | 9.8 ms | 1.0x |
| This package | 0.4 ms | **~25x** |

## Architecture

The package automatically detects CPU features at module load time and selects the best available implementation:

```
┌─────────────────────────────────────────────────────────┐
│                    Python API                           │
│         upsample() / upsample_transposed()              │
└─────────────────────────────────────────────────────────┘
                           │
                           ▼
┌─────────────────────────────────────────────────────────┐
│                  Runtime Dispatch                        │
│         cpu_detect() → select best implementation        │
└─────────────────────────────────────────────────────────┘
                           │
          ┌────────────────┼────────────────┐
          ▼                ▼                ▼
    ┌──────────┐    ┌──────────┐    ┌──────────┐
    │ AVX-512  │    │   NEON   │    │  Scalar  │
    │  AVX2    │    │  (ARM)   │    │(fallback)│
    │   AVX    │    └──────────┘    └──────────┘
    │  SSE4.1  │
    │  SSE2    │
    │  (x86)   │
    └──────────┘
```

## Development

### Environment

Any C toolchain and Python >= 3.10 will do:

```bash
python -m venv .venv
source .venv/bin/activate
pip install -e ".[dev]"
```

Alternatively, a [nix flake](flake.nix) provides a self-contained
development environment (C toolchain, CMake, Ninja, Python, uv):

```bash
nix develop
uv venv --seed
uv pip install -e ".[dev]"
```

### Building and testing

```bash
make build         # rebuild the C extension and reinstall (editable, incremental)
make test          # run the test suite
make lint          # linter
make format        # code formatter
make type-check    # static type checker
make bench         # downsample benchmark vs scipy (needs pip install ".[bench]")
```

The C extension is built with scikit-build-core + CMake; per-kernel SIMD
flags live in [CMakeLists.txt](CMakeLists.txt).

### Project structure

```
fast-resample-cpu/
├── csrc/
│   ├── cpu_detect.c      # Runtime CPU feature detection
│   ├── dispatch.c        # Function pointer dispatch table
│   ├── resample_ext_simd.c  # Python extension wrapper
│   └── kernels/
│       ├── convolve_scalar.c   # Baseline implementation
│       ├── convolve_sse2.c     # x86 SSE2
│       ├── convolve_sse4.c     # x86 SSE4.1
│       ├── convolve_avx.c      # x86 AVX
│       ├── convolve_avx2.c     # x86 AVX2+FMA
│       ├── convolve_avx512.c   # x86 AVX-512
│       └── convolve_neon.c     # ARM NEON
├── src/
│   └── sgnl_cpu_interp/  # Python API
├── CMakeLists.txt        # C extension build configuration
├── pyproject.toml        # Package metadata and build backend
└── tests/                # Test suite
```

### Adding a new SIMD implementation

1. Create `csrc/kernels/convolve_<name>.c` implementing `convolve_<name>()` and `convolve_transposed_<name>()`
2. Add the implementation to the dispatch table in `csrc/dispatch.c`
3. Add per-file compile flags in `CMakeLists.txt` (see the existing `set_source_files_properties` calls)
4. Add CPU feature detection if needed in `csrc/cpu_detect.c`

## Algorithm

This implementation uses polyphase filtering for efficient upsampling:

1. **Kernel generation**: Creates a Lanczos-windowed sinc kernel and splits it into `factor` polyphase components
2. **SIMD convolution**: Vectorized dot product across channels (standard layout) or time samples (transposed layout)
3. **Phase-blocked upsampling**: Processes all output samples with the same phase together to maximize kernel data reuse in cache

The approach is specifically optimized for:
- Many channels (100+)
- Small to moderate upsampling factors (2-16x)
- Short to medium input lengths (100s to 1000s of samples)

## License

MIT License - see LICENSE file for details.

## Contributing

Contributions welcome! Please open an issue or merge request on [git.ligo.org](https://git.ligo.org/greg/fast-resample-cpu).

## Citation

If you use this in research, please cite:
```
@software{sgnl_cpu_interp,
  title = {sgnl-cpu-interp: Fast polyphase resampling with multi-architecture SIMD},
  url = {https://git.ligo.org/greg/fast-resample-cpu},
  year = {2025}
}
```
