# calibrate_weight_sine

## Overview

`calibrate_weight_sine` provides comprehensive, production-ready ADC foreground calibration using sinewave input. This is the full-featured version supporting automatic frequency search, rank deficiency handling, harmonic nuisance fitting, and multi-dataset calibration.

**Key Features:**
- ✅ Automatic frequency search with selectable coarse-estimator policies and fine refinement
- ✅ Rank deficiency handling for redundant ADC architectures
- ✅ Harmonic nuisance fitting for source/test-chain harmonic terms
- ✅ Multi-dataset calibration for time-interleaved ADCs
- ✅ Numerical conditioning for robust convergence
- ✅ Comprehensive diagnostics (weights, fitted residual metrics, error signals)

## Syntax

```python
from adctoolbox import analyze_spectrum
from adctoolbox.calibration import calibrate_weight_sine, scale_calibration_output

# Basic usage (auto frequency search)
result = calibrate_weight_sine(bits)

# MATLAB wcalsin-compatible coarse frequency search
result = calibrate_weight_sine(bits, freq=0, frequency_policy="matlab")

# With known frequency
result = calibrate_weight_sine(bits, freq=0.001587)

# With redundant weights
result = calibrate_weight_sine(bits, freq=0.001587,
                               nominal_weights=[2048, 1024, 512, 256, 128, 128, ...])

# With harmonic nuisance fitting
result = calibrate_weight_sine(bits, freq=0.001587, harmonic_order=3)

# Multi-dataset (time-interleaved ADC)
result = calibrate_weight_sine([bits_ch0, bits_ch1, bits_ch2, bits_ch3])
```

## Parameters

- **`bits`** (ndarray or list of ndarrays) — Binary data matrix
  - Single dataset: (N samples × M bits) ndarray
  - Multi-dataset: list of ndarrays for time-interleaved ADCs
  - Each row is one sample, each column is a bit (MSB first)

- **`freq`** (float, array, or None, optional) — Normalized frequency (f_in / f_s)
  - `None`: Automatic frequency search (default)
  - `float`: Single frequency for all datasets
  - `array`: Per-dataset frequencies for multi-dataset mode

- **`force_search`** (bool or None, optional) — Fine frequency search policy
  - Default: `None`
  - `None`: refine automatically estimated frequencies; keep explicitly provided nonzero frequencies fixed
  - `True`: refine provided frequencies too
  - `False`: disable fine search unless a zero frequency placeholder remains

- **`frequency_policy`** (`"python"` or `"matlab"`, optional) — Coarse estimator used for automatic frequency search
  - Default: `"python"` preserves the historical Python bit-activity/SNDR coarse estimator
  - `"matlab"` uses the MATLAB `wcalsin(freq=0)` compatible estimator: rank-patched nominal bit-prefix reconstructions, one frequency fit per prefix, then the median candidate
  - Explicit nonzero `freq` values remain fixed unless `force_search=True`

- **`nominal_weights`** (array, optional) — Nominal bit weights (only effective when rank is deficient)
  - Default: `[2^(M-1), 2^(M-2), ..., 2, 1]`
  - Required for redundant ADCs (e.g., `[128, 128, 64, ...]`)

- **`harmonic_order`** (int, optional) — Number of harmonic terms included in the fitted reference
  - Default: `1` (fundamental only)
  - Higher values include H2/H3/... as fitted nuisance terms in `ideal`
  - Fitted harmonic terms are excluded from `error`
  - Example: `harmonic_order=3` fits the fundamental plus H2/H3 in the reference
  - Use this to keep source/test-chain harmonics from contaminating weight estimation; it does not prove ADC harmonic distortion was removed from `calibrated_signal`

- **`learning_rate`** (float, optional) — Adaptive learning rate for frequency updates
  - Default: `0.5`
  - Range: `0 < learning_rate < 1`
  - Lower values = more conservative convergence

- **`reltol`** (float, optional) — Relative error tolerance for convergence
  - Default: `1e-12`

- **`max_iter`** (int, optional) — Maximum iterations for fine frequency search
  - Default: `100`

- **`verbose`** (int, optional) — Print verbosity level
  - `0`: Silent (default)
  - `1`: Basic progress
  - `2`: Detailed diagnostics

## Returns

Dictionary with keys:

- **`weight`** (ndarray) — Calibrated bit weights, normalized by sinewave magnitude
  - Length: M (bit width)
  - Normalized so max weight ≈ 1.0

- **`offset`** (float or ndarray) — Calibrated DC offset(s)
  - Single value for single dataset
  - Array for multi-dataset

- **`calibrated_signal`** (ndarray or list) — Signal after calibration
  - Single array for single dataset
  - List of arrays for multi-dataset

- **`ideal`** (ndarray or list) — Best fitted reference waveform
  - Includes the fundamental and fitted harmonics up to `harmonic_order`
  - Single array for single dataset
  - List of arrays for multi-dataset

- **`error`** (ndarray or list) — Residue errors after calibration
  - `error = calibrated_signal - offset - ideal`
  - Excludes the harmonic terms included in `ideal`

- **`refined_frequency`** (float or ndarray) — Refined frequency from calibration
  - Single value for single dataset
  - Array for multi-dataset

- **`initial_frequency`** (float or ndarray) — Coarse frequency used before fine search
  - Useful for diagnosing automatic frequency-search policy differences

- **`frequency_policy`** (str) — Coarse-estimator policy used for this result

- **`snr_db`** (float or ndarray) — Fitted-reference to residual ratio in dB
  - When `harmonic_order > 1`, this is not standard ADC dynamic SNDR

- **`enob`** (float or ndarray) — Effective number of bits
  - Calculated as: `(snr_db - 1.76) / 6.02`
  - Interpret as a fitted-residual ENOB estimate, not FFT ENOB, when harmonics are included in `ideal`

- **`rank_patch`** (dict) — Rank-deficiency diagnostics
  - `applied`: whether rank-deficiency patching was used
  - `bit_width_effective`: number of effective columns used by the solver
  - `bit_to_col_map`: original bit column to effective column map (`-1` means unmapped)
  - `bit_weight_ratios`: nominal-weight ratios used for merged columns
  - `dropped_constant_bits`: bit columns that were constant in this capture
  - `unmapped_bits`: bit columns returned as zero for this fitted model

## Output Scale Convention

`calibrate_weight_sine` returns calibrated waveforms in **solver-unit-sine**
scale. This is a mathematical gauge used to make the least-squares problem
identifiable: the fitted fundamental sine magnitude is fixed to one. It is not
automatically the ADC voltage or code full-scale.

The result dictionary includes `scale_convention = "solver_unit_sine"`. Ratio
metrics in the calibration result, such as `snr_db` and `enob`, are unchanged
by a later linear rescale. Absolute full-scale-referenced spectrum metrics,
such as `sig_pwr_dbfs`, `noise_floor_dbfs`, and `nsd_dbfs_hz`, require an
explicit ADC/code reference scale before calling `analyze_spectrum`.
When `max_scale_range=None`, the spectrum analyzer uses the waveform's own
range as a self-reference; that keeps ratio metrics valid but does not recover
the physical ADC full-scale. If an explicit `max_scale_range` is supplied and
the waveform exceeds it, the analyzer emits an overrange warning without
clipping or changing the data.

Use `scale_calibration_output` to map the solver result to a chosen convention:

```python
result = calibrate_weight_sine(bits, freq=freq_true, nominal_weights=w_nominal)

# Map back to the same code/voltage convention as the nominal reconstruction
# weights before interpreting dBFS or NSD.
result_adc = scale_calibration_output(result, target_weights=w_nominal)
metrics = analyze_spectrum(
    result_adc["calibrated_signal"][0],
    max_scale_range=(0.0, 1.0),
    create_plot=False,
)
```

If the training sine peak is known directly in ADC units, use
`target_sine_peak=A` instead. Do not treat `1/A` as a universal weight-scale
rule; the correct bridge depends on the bit encoding and ADC convention.

## Algorithm

### Modular Pipeline (8 Stages)

```
Input → Prepare → Patch Rank → Scale Columns → Estimate Freq
   ↓
Solve (Freq Search or Static) → Recover Scaling → Recover Rank → Post-process → Output
```

Each stage is implemented as a separate helper function.

### Mathematical Model

The weighted bit output is fitted to a reference waveform with optional
harmonic nuisance terms:

```
y(n) = Σ w_i · b_i(n) = Σ [A_k·cos(2πkfn) + B_k·sin(2πkfn)] + C
                        k=1 to H
```

where:
- `y(n)` = reconstructed analog signal
- `w_i` = weight of bit i (unknown)
- `b_i(n) ∈ {0, 1}` = binary value of bit i
- `f` = normalized frequency
- `H` = harmonic order
- `A_k, B_k` = fitted reference coefficients for harmonic k
- `C` = DC offset

For `H > 1`, the harmonics are part of the fitted reference. They are useful
when the training source or test chain contains harmonic content that should not
be absorbed into ADC bit weights. They should not be read as a claim that ADC
harmonic distortion has been removed from `calibrated_signal`.

### Dual-Basis Least Squares

Unlike the lite version, the full algorithm tries **both** basis assumptions:

**Cosine Basis** (A_1 = 1):
```
[B | 1 | sin(2πfn) | ... | sin(2πHfn) | cos(2π·2fn) | ...] × [w; C; B_1; ...; B_H; A_2; ...] = -cos(2πfn)
```

**Sine Basis** (B_1 = 1):
```
[B | 1 | cos(2πfn) | ... | cos(2πHfn) | sin(2π·2fn) | ...] × [w; C; A_1; ...; A_H; B_2; ...] = -sin(2πfn)
```

The algorithm solves both and selects the one with **lower residual error**.

### Rank Deficiency Handling

For redundant ADCs (e.g., weights `[128, 128, 64, ...]`), the bit matrix has rank deficiency. The algorithm:

1. **Identifies redundant columns** — Bits with identical patterns
2. **Groups identical bits** — `{b_i, b_j}` if `b_i(n) = b_j(n)` for all n
3. **Solves reduced system** — One representative per group
4. **Distributes weights** using nominal weights:
   ```
   w_i = w_eff^(group) × (w_i^nom / Σ w_j^nom)
                              j in group
   ```

**Example:**
- Redundant bits: `[b_4, b_5]` both have pattern `[0,1,1,0,1,...]`
- Nominal weights: `w_4^nom = 128, w_5^nom = 128`
- Solved effective weight: `w_eff = 256`
- Recovered weights: `w_4 = w_5 = 256 × (128/(128+128)) = 128` ✅

**This preserves redundancy rather than collapsing it to zero!**

### Frequency Search

**Coarse Search** (when `freq=None` or `freq=0`):
1. `frequency_policy="python"` preserves the historical Python estimator:
   it ranks bit columns by toggle activity, builds low-toggle and high-toggle
   surrogate waveforms, chooses the cleaner surrogate by spectrum quality, and
   fits its tone frequency.
2. `frequency_policy="matlab"` matches MATLAB `wcalsin(freq=0)` coarse search:
   after rank patching, reconstruct prefixes of the first `min(M, 5)`
   effective bit columns using nominal representative weights, fit one
   frequency per prefix, and use the median candidate.
3. Explicit nonzero `freq` values skip coarse search. They remain fixed by
   default and are refined only when `force_search=True`.

**Fine Search** (iterative refinement):
```
for iteration = 1 to max_iter:
    1. Solve weights at current frequency f^(t)
    2. Append the derivative of the fitted harmonic reference with respect to frequency
    3. Solve the augmented least-squares system for a frequency correction
    4. Update frequency: f^(t+1) = f^(t) + learning_rate * delta_f
    5. Check convergence: |f^(t+1) - f^(t)| / f^(t) < reltol
    6. If converged, break
```

## Examples

### Example 1: Basic Calibration (Known Frequency)

```python
import numpy as np
from adctoolbox.calibration import calibrate_weight_sine

# Generate test data
n_samples = 8192
bit_width = 12
freq_true = 13 / n_samples

# Create sinewave and quantize
t = np.arange(n_samples)
signal = 0.5 * np.sin(2 * np.pi * freq_true * t) + 0.5
quantized = np.clip(np.floor(signal * (2**bit_width)), 0, 2**bit_width - 1).astype(int)

# Extract bits
bits = (quantized[:, None] >> np.arange(bit_width - 1, -1, -1)) & 1

# Calibrate with known frequency
result = calibrate_weight_sine(bits, freq=freq_true, verbose=2)

# Access results
print(f"Recovered weights: {result['weight']}")
print(f"Refined frequency: {result['refined_frequency']:.8f}")
print(f"Fitted residual SNR: {result['snr_db']:.2f} dB")
print(f"Fitted residual ENOB: {result['enob']:.2f} bits")
```

### Example 2: Automatic Frequency Search

```python
# Calibrate without knowing frequency
result = calibrate_weight_sine(bits, freq=None, verbose=2)

# Algorithm will:
# 1. Estimate frequency from FFT
# 2. Refine using iterative search
# 3. Return refined frequency
print(f"Estimated frequency: {result['refined_frequency']:.8f}")
print(f"True frequency:      {freq_true:.8f}")
print(f"Error:               {abs(result['refined_frequency'] - freq_true):.2e}")
```

### Example 3: Redundant ADC Calibration

```python
# Redundant weights: [2048, 1024, 512, 256, 128, 128, 64, ...]
true_weights = np.array([2048, 1024, 512, 256, 128, 128, 64, 32, 16, 8, 4, 2])

# Generate redundant bit data using greedy decomposition
def decompose_to_redundant_bits(codes, weights):
    n_samples = len(codes)
    bit_width = len(weights)
    bits = np.zeros((n_samples, bit_width), dtype=int)

    for i, code in enumerate(codes):
        remaining = float(code)
        for j, weight in enumerate(weights):
            if remaining >= weight - 0.5:
                bits[i, j] = 1
                remaining -= weight

    return bits

bits_redundant = decompose_to_redundant_bits(quantized, true_weights)

# Provide nominal weights as hints
result = calibrate_weight_sine(bits_redundant, freq=freq_true,
                               nominal_weights=true_weights)

# Both redundant bits will be recovered correctly
print(f"True weights:      {true_weights}")
print(f"Recovered weights: {result['weight'] * np.max(true_weights)}")
# Expected: [2048, 1024, 512, 256, 128, 128, 64, 32, 16, 8, 4, 2] ✅
# NOT:      [2048, 1024, 512, 256, 128,   0, 64, 32, 16, 8, 4, 2] ❌
```

### Example 4: Harmonic Nuisance Fitting

```python
from adctoolbox import analyze_spectrum

# Fit source/test-chain harmonic nuisance terms up to H3.
result = calibrate_weight_sine(
    bits,
    freq=freq_true,
    harmonic_order=3,  # Model fundamental + 2nd + 3rd harmonics
    verbose=2
)

# The returned snr_db is a fitted-residual metric. Fitted H2/H3 terms are
# included in result["ideal"] and excluded from result["error"].
print(f"Fitted residual SNR: {result['snr_db']:.2f} dB")

# To evaluate ADC dynamic performance, analyze the calibrated output itself.
spec = analyze_spectrum(result["calibrated_signal"])
print(f"FFT SNDR: {spec['sndr_dbc']:.2f} dBc")
```

### Example 5: Multi-Dataset Calibration (Time-Interleaved ADC)

```python
# Time-interleaved ADC with 4 channels
bits_ch0 = ...  # Channel 0 data (N0 samples × M bits)
bits_ch1 = ...  # Channel 1 data (N1 samples × M bits)
bits_ch2 = ...  # Channel 2 data (N2 samples × M bits)
bits_ch3 = ...  # Channel 3 data (N3 samples × M bits)

# Calibrate all channels together (shared weights)
result = calibrate_weight_sine(
    bits=[bits_ch0, bits_ch1, bits_ch2, bits_ch3],
    freq=None,  # Auto frequency search for each channel
    verbose=2
)

# Results are lists for multi-dataset
print(f"Shared weights: {result['weight']}")
print(f"Per-channel offsets: {result['offset']}")
print(f"Per-channel frequencies: {result['refined_frequency']}")
print(f"Per-channel fitted residual ENOB: {result['enob']}")
```

### Example 6: Forced Frequency Refinement

```python
# Even with provided frequency, force refinement
result = calibrate_weight_sine(
    bits,
    freq=13/8192,  # Approximate frequency
    force_search=True,  # Refine it anyway
    learning_rate=0.3,  # Conservative learning rate
    reltol=1e-14,  # Tight convergence
    max_iter=200,  # Allow more iterations
    verbose=2
)

print(f"Initial frequency: {13/8192:.8f}")
print(f"Refined frequency: {result['refined_frequency']:.8f}")
```

## Performance

### Computational Complexity

**Time Complexity**: `O(I·N·M² + M³)`
- I = number of frequency search iterations
- N = number of samples
- M = bit width

**Space Complexity**: `O(N·M)`

**Typical Performance** (12-bit ADC, N=8192, Intel i7):
- Known frequency: ~20 ms
- Frequency search: ~50-100 ms (depends on convergence)
- Multi-dataset (4 channels): ~80-150 ms

### Accuracy

**Weight Recovery**:
- Ideal conditions: Error < 10⁻⁶ LSB
- With noise (SNR > 60 dB): Error < 10⁻³ LSB
- Redundant weights: Fully preserved (both bits active)

**Frequency Estimation**:
- Coarse FFT: Accuracy ≈ 1/N bins
- Fine search: Accuracy < 10⁻¹² (relative)

**Fitted-Residual ENOB Improvement**:
- Binary ADC (no INL): +0.5 to +2 ENOB
- Redundant ADC: +2 to +4 ENOB (error correction)
- Time-interleaved: +3 to +6 ENOB (timing skew correction)

These values are based on the fitted residual returned by calibration. When
`harmonic_order > 1`, fitted harmonics are excluded from that residual. Use
`analyze_spectrum(result["calibrated_signal"])` for standard FFT SNDR, THD, HDx,
and dynamic ENOB.

## Limitations

### Input Signal Requirements

1. **Amplitude**: Input should be > -6 dBFS for stable weight recovery
2. **Purity**: Input signal should have low distortion (THD < -60 dB) for best results. If source/test-chain harmonics are unavoidable, `harmonic_order > 1` can model them as nuisance terms for weight estimation, but the calibrated output spectrum remains the source of truth for ADC distortion.
3. **Coherency**: For FFT-based frequency estimation, use coherent sampling when possible
4. **Bit activity**: Each weight can only be estimated if its bit column has AC activity in the capture

### Frequency Constraints

1. **Nyquist**: `0 < f < 0.5` (normalized)
2. **Avoid DC and Nyquist**: `0.01 < f < 0.49` recommended
3. **Multi-tone interference**: Avoid frequencies near `k/M` where k is small integer

### Convergence Issues

Frequency search may fail to converge if:
1. Initial frequency estimate is very poor (> 10% error)
2. Input amplitude is too low (< -10 dBFS)
3. Learning rate is too aggressive (α > 0.8)

**Solutions:**
- Provide better initial frequency estimate
- Reduce `learning_rate` to 0.1-0.3
- Increase `max_iter` to 200-500
- Check `verbose=2` output for convergence diagnostics

### Multi-Dataset Assumptions

Multi-dataset calibration assumes:
1. All channels share same bit weights (valid for time-interleaved ADCs)
2. Independent offset and gain per channel
3. Small timing skew between channels

For ADCs with **per-channel weight variation**, use separate single-dataset calibrations.

### Rank-Deficient Captures

Constant bit columns have no AC information and are not identifiable separately
from the fitted offset. If all bit columns are constant, calibration raises a
`ValueError` because no effective bit columns remain. If only some columns are
constant or otherwise unmapped, the solver continues with observable columns,
sets the unmapped fitted weights to zero for the current model, emits a
`UserWarning`, and reports the affected columns in `result["rank_patch"]`.

Those zero fitted weights do **not** mean the physical ADC weights are zero.
They mean the current capture did not observe those columns. Increase input
amplitude, improve code coverage, use a different common-mode/window, or combine
multiple captures if those weights must be calibrated.

## Comparison with Lite Version

| Feature | Lite Version | Full Version |
|---------|-------------|--------------|
| **Frequency handling** | Known frequency only | Auto search + refinement |
| **Rank deficiency** | ❌ Collapses redundancy | ✅ Preserves all weights |
| **Harmonic nuisance fitting** | ❌ No | ✅ Configurable order |
| **Multi-dataset** | ❌ Single only | ✅ Multiple datasets |
| **Numerical stability** | Basic lstsq | Column scaling + conditioning |
| **Output** | Weights only (ndarray) | Full diagnostics (dict) |
| **Return values** | 1 (weights) | weights, offset, signals, fitted residual metrics, rank patch diagnostics, etc. |
| **Code complexity** | ~40 lines | ~600+ lines (modular) |
| **Typical runtime** | 5 ms | 20-100 ms |

### When to Use Full Version

✅ **Use full version when:**
- Unknown or imprecise frequency
- Redundant ADC architecture
- Need comprehensive diagnostics
- Multi-dataset or time-interleaved ADCs
- Harmonic nuisance fitting required for source/test-chain distortion
- Production/research environments

❌ **Use lite version when:**
- Frequency precisely known
- Binary weighted ADC (no redundancy)
- Speed critical (embedded systems)
- Minimal code footprint needed

## Pipeline Details

### Stage 1: Input Preparation (`_prepare_input`)

**Purpose**: Normalize input to unified format and validate

**Operations**:
1. Convert single array to list: `[bits]`
2. Validate all arrays have same bit width M
3. Stack: `bits_stacked = [bits_1; bits_2; ...]`
4. Generate nominal weights if not provided: `w_i^nom = 2^(M-1-i)`

### Stage 2: Rank Deficiency Patching (`_patch_rank_deficiency`)

**Purpose**: Detect and handle redundant bits

**Algorithm**:
1. Compute column norms: `||b_i|| = sqrt(Σ b_i(n)²)`
2. Find identical columns: `||b_i - b_j|| < ε`
3. Group identical bits: `G_1, G_2, ..., G_K` where K = M_eff
4. Select representatives for each group
5. Compute weight ratios: `r_i = w_i^nom / Σ w_j^nom` (j in group)

### Stage 3: Column Scaling (`_scale_columns_for_conditioning`)

**Purpose**: Normalize bit columns for numerical stability

**Operations**:
```
s_i = ||b_i||_2
b̃_i = b_i / s_i
```

### Stage 4: Frequency Estimation (`_estimate_frequencies`)

**Purpose**: Estimate or validate input frequencies

**Logic**:
- If `freq=None` or `freq=0`: coarse estimation for each dataset using
  `frequency_policy`
- If `freq=scalar`: broadcast to all datasets
- If `freq=array`: validate length matches number of datasets

### Stage 5: Least Squares Solve

**Purpose**: Solve for weights and sinewave parameters

**Known Frequency Path**:
1. Build design matrix A with all harmonics
2. Try both cosine and sine basis
3. Solve both and select basis with lower residual

**Frequency Search Path**:
1. Initialize with coarse frequency
2. Loop until convergence:
   - Solve weights at current frequency
   - Append the frequency-derivative column
   - Solve the augmented least-squares system for a correction
   - Check convergence

### Stage 6-8: Recovery and Post-Processing

- **Recover column scaling**: Undo normalization
- **Recover rank deficiency**: Distribute effective weights to all physical bits
- **Post-process**: Assemble results dictionary with fitted residual metrics and error signals

## References

1. **IEEE Standard 1057-2017**: "IEEE Standard for Digitizing Waveform Recorders"

2. Vogel, C., & Johansson, H. (2006). "Time-interleaved analog-to-digital converters: Status and future directions." *IEEE Transactions on Circuits and Systems I*, 53(11), 2386-2394.

3. Jin, H., & Lee, E. K. F. (1992). "A digital-background calibration technique for minimizing timing-error effects in time-interleaved ADCs." *IEEE Transactions on Circuits and Systems II*, 47(7), 603-613.

4. Le Dortz, N., et al. (2014). "A 1.62GS/s time-interleaved SAR ADC with digital background mismatch calibration achieving interleaving spurs below 70dBFS." *ISSCC Digest of Technical Papers*, 386-387.

5. MATLAB Signal Processing Toolbox: `sinefitweights` documentation

## See Also

- [calibrate_weight_sine_lite](calibrate_weight_sine_lite.md) - Lightweight version for embedded systems
- [analyze_spectrum](analyze_spectrum.md) - Spectral analysis for SNDR/ENOB
- [fit_sine_4param](fit_sine_4param.md) - IEEE 1057 sinewave fitting