matlab-fit-rational-model

# Rational Fitting and Time-Domain Analysis

Fit S-parameter data to rational function models using the `rational` object, then compute time-domain responses (TDR, impulse, step) and export to SPICE/Verilog-A.

## When to Use

- Fitting S-parameter data to rational function models (vector fitting, AAA)
- Computing time-domain responses: TDR, impulse, step, arbitrary waveform
- Exporting broadband circuit models to SPICE or Verilog-A
- Verifying fit accuracy in frequency domain with `freqresp`
- Enforcing causality or passivity on fitted models
- Extracting poles, zeros, and state-space matrices from fits

## When NOT to Use

- Loading or plotting S-parameter data -- use `matlab-manage-sparameters`
- Converting between network parameter types -- use `matlab-convert-network-parameters`
- RF budget or cascade analysis -- use `matlab-analyze-rf-budget` or `matlab-deembed-rf-cascade`

## Workflow

1. **Load S-parameters** — See `matlab-manage-sparameters` skill
2. **Fit** — Create a `rational` model from the S-parameter data
3. **Verify** — Compare fit against original data in frequency domain
4. **Analyze** — Compute time-domain responses (TDR, impulse, step)
5. **Export** — Generate SPICE subcircuit or Verilog-A model

## The `rational` Object (R2020a+)

Prefer the `rational` object (AAA algorithm) over `rationalfit` (vector-fitting algorithm). Both provide the same analysis methods (`impulse`, `stepresp`, `timeresp`, `zpk`, `generateSPICE`, `pwlresp`, `iscausal`, `ispassive`, `makepassive`, `passivity`, `freqresp`, `abcd`, `writeva`). The advantage of `rational` is that it tends to converge more robustly. The advantage of `rationalfit` is that it enables delay removal.

**Note:** `rationalfit` still ships and uses the vector-fitting algorithm. It returns the legacy `rfmodel.rational` object. Use `rational` (without the "fit" suffix) for new code unless you specifically need vector fitting.

### Fitting S-Parameters

```matlab
% Fit entire S-parameter object (all ports)
[fit, errdb] = rational(s);
fprintf('Ports: %d, Poles: %d, Error: %.1f dB\n', fit.NumPorts, fit.NumPoles, errdb);

% Fit a single parameter from (freq, data)
s21 = rfparam(s, 2, 1);
[fit21, err21] = rational(s.Frequencies, s21);
```

### Tuning the Fit

```matlab
[fit, errdb] = rational(s, ...
    'Tolerance', -50, ...              % Target accuracy in dB
    'MaxPoles', 30);                   % Limit pole count
```

| Option | Default | Purpose |
|--------|---------|---------|
| `Tolerance` | -40 | Target error in dB |
| `MaxPoles` | 1000 | Upper bound on poles |
| `Causal` | true | Enforce causal (stable) poles |
| `TendsToZero` | true | Force fit to decay to zero as frequency → ∞ |
| `NoiseFloor` | -Inf | Ignore data below this level (dB) |
| `Display` | `'off'` | Set to `'plot'` to visualize fit progress |

**`TendsToZero`:** When `true` (default), the rational model has strictly proper form (numerator order < denominator order), so the response decays to zero at high frequencies. This is appropriate for bandpass or lowpass responses (e.g., insertion loss S21 of a lossy channel). Set `TendsToZero=false` for responses that do not vanish at infinity (e.g., reflection coefficients, allpass networks, or broadband matched loads where S11 plateaus). The `false` setting produces a more general model with a nonzero direct-feedthrough term (D ≠ 0).

## Verify the Fit

```matlab
% Frequency response of the fit
resp = freqresp(fit, s.Frequencies);   % Returns N×N×K for N-port fit

% For multi-port fit, extract S21 element
figure;
plot(s.Frequencies/1e9, 20*log10(abs(rfparam(s, 2, 1))), ...
     s.Frequencies/1e9, 20*log10(abs(squeeze(resp(2,1,:)))), '--');
xlabel('Frequency (GHz)'); ylabel('Magnitude (dB)');
legend('Original', 'Rational Fit'); title('S21 Fit Verification');
```

### Extract Poles, Zeros, and State-Space

```matlab
% Direct property access (preferred)
poles = fit.Poles;                     % Pole vector
nPoles = fit.NumPoles;                 % Number of poles

% zpk form
[z, p, k] = zpk(fit);                 % Zeros, poles, gain
[A, B, C, D] = abcd(fit);             % State-space matrices
```

**Gotcha:** Use `fit.NumPoles` and `fit.Poles` for direct access. The `zpk(fit)` form returns a struct with cell arrays (`.Z`, `.P`, `.K`), which requires indexing like `z.P{1}` — the direct properties are simpler.

### Causality and Passivity Enforcement

```matlab
% Step 1: Fit with causal poles enforced (default Causal=true)
fit = rational(s, Causal=true);

% Step 2: Verify causality and passivity
iscausal(fit)                          % logical: true if all poles in left half-plane
ispassive(fit)                         % logical: true if passive at all frequencies

% For numeric quality metrics per IEEE P370 (R2023b+), use on sparameters:
% [cm, rm, pm] = ieee370QualityCheckFrequencyDomain(s)

% Step 3: Visualize passivity across frequency
figure;
passivity(fit);                        % Plots H-infinity norm vs frequency
                                       % Passive where norm <= 1

% Step 4: Enforce passivity if needed
if ~ispassive(fit)
    fit = makepassive(fit);            % Perturbs residues to enforce passivity
    passivity(fit);                    % Re-plot to confirm
end
```

**`Causal=true`** (the default) constrains all poles to the left half-plane during fitting, ensuring a stable/causal model. Set `Causal=false` only when fitting active devices or when stability is not required.

**`passivity(fit)`** plots the 2-norm of the transfer function vs frequency — the fit is passive wherever the norm is at or below 1. With output arguments, `pNorm = passivity(fit)` returns the scalar maximum passivity norm. Call without arguments for reliable results — `passivity(fit, freqVector)` with a custom frequency vector is not supported in all releases.

**`makepassive(fitObj)`** enforces passivity by running an optimization that perturbs the fit residues. This is different from `makepassive(sparametersObj)` on raw S-par