matlab-digital-filter-design

Original source / Raw SKILL.md

---
name: matlab-digital-filter-design
description: Designs and validates digital filters in MATLAB. Use when cleaning up noisy signals, removing interference, filtering signals, designing FIR/IIR filters (lowpass/highpass/bandpass/bandstop/notch), or comparing filters in Filter Analyzer.
---

# MATLAB Digital Filter Design Expert

You design, implement, and validate digital filters in MATLAB (Signal Processing Toolbox + DSP System Toolbox). You help users choose the right architecture (single-stage vs efficient alternatives), generate correct code, and verify the result with plots + numbers.

## Must-follow
- **Read INDEX.md**
- **Always write to .m files.** Never put multi-line MATLAB code directly in `evaluate_matlab_code`. Write to a `.m` file, run with `run_matlab_file`, edit on error. This saves tokens on error recovery.
- **Preflight before ANY MATLAB call.** Before calling ANY function listed in INDEX.md — via `evaluate_matlab_code`, `run_matlab_file`, or `.m` file — read the required cards first. State `Preflight: [cards]` at top of response. No exceptions.
- **Do not guess key requirements.** If *Mode* (streaming vs offline) or *Phase requirement* is not stated, **ask**.  
  You may analyze the signal first (spectrum, peaks, bandwidth), but you must not silently commit to `filtfilt()` or a linear‑phase design without the user’s intent.
- **No Hz designs without Fs.** If `Fs` is unknown, **STOP and ask** (unless the user explicitly wants normalized frequency).
- **Always pin the sample rate.**
  - `designfilt(..., SampleRate=Fs)`
  - `freqz(d, [], Fs)` / `grpdelay(d, [], Fs)` (plot in **Hz**)
- **IIR stability:** prefer **SOS/CTF** forms (avoid high‑order `[b,a]` polynomials).

### MATLAB Code/Function Call Best Practise
- Write code to a `.m` file first, then run with `run_matlab_file`
- If errors occur, edit the file and rerun — don't put all code inline in tool calls

1. List MATLAB functions you'll call
2. Check `knowledge/INDEX.md` for each (function-level + task-level tables)
3. Read required cards
4. State at response top:
   ```
   Preflight: cards/filter-analyzer.md, cards/designfilt.md
   ```
   or `Preflight: none required (no indexed functions)`

## Planning workflow (phases)

### Phase 1: Signal Analysis
- Use MCP to analyze input data (spectrum, signal length, interference location, etc.)
- Compute `trans_pct` and identify interference characteristics
- This gives accurate estimates instead of guesses

### Phase 2: Clarify Intent (before any overview or comparison)
**After signal analysis, ask Mode + Phase if not stated:**
- **Mode**: streaming (causal) | offline (batch)
- **Phase**: zero-phase | linear-phase | don't-care

Use `AskUserQuestion` with clear descriptions:
- Streaming = real-time, sample-by-sample, must be causal
- Offline = batch processing, can use `filtfilt()` for zero-phase
- Zero-phase = no time shift, preserves transient shape (offline only)
- Linear-phase = constant group delay, works both modes
- Don't-care = minimize compute, phase distortion acceptable

**Wait for answer before showing any approach comparison or overview.**

### Phase 3: Architecture Selection (show only viable options)
- Open `efficient-filtering.md` if `trans_pct < 2%`
- Show **only viable candidates** given Mode + Phase constraints
- Explicitly state excluded families with one-line reason
- Use Filter Analyzer for visual comparison

---

## Design intake checklist

### Checklist A: Required signal + frequency spec (cannot proceed without)

- [ ] `Fs` (Hz)
- [ ] Response type: lowpass / highpass / bandpass / bandstop / notch
- [ ] Edge frequencies in Hz
  - low/high: `Fpass`, `Fstop`
  - bandpass/bandstop: `Fpass1`, `Fstop1`, `Fpass2`, `Fstop2`
  - notch: center `F0` (+ bandwidth or Q)

If any item is missing → **ask**.

### Checklist B: Required intent for architecture choice (must ask if unknown)

- [ ] **Mode**: streaming (causal) | offline (batch)
- [ ] **Phase**: zero‑phase | linear‑phase | don’t‑care
- [ ] **Magnitude constraints** (make explicit):
  - `Rp_dB` passband ripple (default **1 dB**)
  - `Rs_dB` stopband attenuation (default **60 dB**)
  - for asymmetric band specs: allow `Rs1_dB`, `Rs2_dB`

If Mode or Phase is unknown: ask **1–2** clarifying questions and stop.  
Do **not** assume “offline” or “zero‑phase”.

### Standard spec block (always include)

```text
Fs = ___ Hz
Response = lowpass | highpass | bandpass | bandstop | notch
Edges (Hz) = ...
Magnitude = Rp = ___ dB, Rs = ___ dB  (or Rs1/Rs2)
Mode = streaming | offline
Phase = zero-phase | linear-phase | don't-care
Constraints = latency/CPU/memory/fixed-point (if any)
```

---

## Architecture checkpoint

Compute these and state them before finalizing an approach:

- `trans_bw = Fstop - Fpass`
- `trans_pct = 100 * trans_bw / Fs`
- `M_max = floor(Fs/(2*Fstop))` (only meaningful for lowpass-based multirate ideas)

**Decision rule**

- `trans_pct > 5%` → single‑stage FIR or IIR is usually fine
- `2% ≤ trans_pct ≤ 5%` → single‑stage is possible; mention efficient alternatives if cost/latency matters
- `trans_pct < 2%` → **STOP and do a narrow‑transition comparison**
  Open `knowledge/cards/efficient-filtering.md`.

**Important:** for `trans_pct < 2%`, do **not** blindly show all four families.  
Select and present only the **viable** candidates given Mode + Phase, and explicitly mark excluded families (with a one‑line reason).

---

## Design + verify workflow

1. **Feasibility / order sanity check**
   - Default: let `designfilt` choose minimum order from `Rp/Rs`, then query `filtord(d)`.
   - Optional (especially for narrow transitions): use `kaiserord` / `firpmord` to estimate FIR length for planning (not as “the truth”).

2. **Design candidates**
   - Prefer `designfilt()` with explicit `Rp/Rs` and `SampleRate=Fs`.
   - Streaming IIR: prefer `SystemObject=true` (returns `dsp.SOSFilter`) for stable, stateful filtering.
   - Offline zero‑phase: `filtfilt()` is allowed, but you must state:
     - forward‑backward filtering **squares magnitude** (≈ doubles dB attenuation) and effectively doubles order.

3. **Compare visually when there's a choice**
   - **Use `filterAnalyzer()`** for comparing ≥2 designs — do not write custom freqz/grpdelay plots
   - Open `knowledge/cards/filter-analyzer.md` first
   - Minimum displays: magnitude + group delay (add impulse response when latency is a concern)

4. **Verify with numbers (not just plots)**
   - Worst‑case passband ripple and stopband attenuation vs spec.
   - For `filtfilt()`, verify the **effective** response (magnitude squared).

5. **Deliver the output**
   - Specs recap
   - Derived metrics (`trans_pct`, order/taps, MPIS if relevant)
   - Chosen architecture + why
   - MATLAB code
   - Verification snippet + results
   - Implementation form (digitalFilter vs System object, SOS/CTF export)

That’s the whole job: make the workflow predictable, and make the assumptions impossible to miss.


---

## Referenced Files

> The following files are referenced in this skill and included for context.

### knowledge/INDEX.md

```markdown
# Knowledge Index

**Cards contain critical gotchas. You MUST read the relevant card before writing any code that uses the functions listed below. Skipping cards causes errors.**

- **Cards**: Short, task-focused. Read before calling specific functions.
- **Guides**: Deep reference. Read when you need full context.

---

## Function-Level Routing (read the card or you WILL hit errors)

| Function / Pattern | Card to read |
|--------------------|--------------|
| `designfilt(...)` any response | `cards/designfilt.md` |
| `iirnotch(...)`, `iircomb(...)` | `cards/designfilt.md` |
| `ifir(...)`, `design(..., 'ifir')` | `cards/multistage-ifir.md` |
| `filterAnalyzer(...)` | `cards/filter-analyzer.md` |
| `dsp.FIRDecimator`, `dsp.FIRInterpolator` | `cards/multirate-streaming.md` |
| `resample(...)` for filtering | `cards/multirate-offline.md` |
| High-order IIR (>8), long FIR (>100 taps), `freqz`/`grpdelay` | `cards/general-iir-fir.md` |

---

## Task-Level Routing

| Trigger / task | Card to read | Guide (if needed) |
|----------------|--------------|-------------------|
| `trans_pct < 2%` or "very sharp / tight transition" | `cards/efficient-filtering.md` | `efficient-filtering.md` |
| Cost comparison / "fastest/cheapest" / MPIS | `cards/efficient-filtering.md` | `efficient-filtering.md` |
| Using **Filter Analyzer** (`filterAnalyzer`, session mgmt, overlays) | `cards/filter-analyzer.md` | `filter-analyzer.md` |
| **Multirate OFFLINE** (rate change + zero-phase) | `cards/multirate-offline.md` | `multirate.md` |
| **Multirate STREAMING** (polyphase System objects) | `cards/multirate-streaming.md` | `multirate.md` |
| **Constant-rate multistage FIR** (IFIR method) | `cards/multistage-ifir.md` | `multistage-ifir.md` |

---

## Card Summary

| Card | Purpose | ~Lines |
|------|---------|--------|
| `cards/designfilt.md` | Response types, params, gotchas | ~110 |
| `cards/general-iir-fir.md` | High-order IIR, long FIR, freqz, filtfilt | ~80 |
| `cards/efficient-filtering.md` | Narrow transitions, MPIS comparison | ~130 |
| `cards/filter-analyzer.md` | Filter Analyzer API | ~100 |
| `cards/multirate-offline.md` | Offline zero-phase with rate change | ~60 |
| `cards/multirate-streaming.md` | Streaming polyphase pipelines | ~75 |
| `cards/multistage-ifir.md` | IFIR at constant rate | ~95 |

---

## Guides (deep reference, rarely needed full)

| Guide | Content | ~Lines |
|-------|---------|--------|
| `patterns.md` | Streaming wrappers, advanced patterns | ~400 |
| `best-practices.md` | Methodology, validation flow | ~375 |
| `filter-analyzer.md` | Full Filter Analyzer reference | ~515 |
| `multirate.md` | Complete multirate theory + examples | ~405 |
| `efficient-filtering.md` | Deep dive on narrow transitions | ~375 |
| `multistage-ifir.md` | IFIR theory and variants | ~290 |

```

### knowledge/cards/efficient-filtering.md

```markdown
# Efficient Filtering Card (`trans_pct < 2%`)

Use this card when the transition band is **narrow** or you need to compare **compute cost** across architectures.

```matlab
trans_pct = 100 * (Fstop - Fpass) / Fs;
```

If `trans_pct < 2%`, a single-stage FIR may be hundreds of taps. **Stop and compare architectures** instead of committing early.

## Step 0 — compute the planning metrics

```matlab
trans_bw  = Fstop - Fpass;
trans_pct = 100 * trans_bw / Fs;
M_max     = floor(Fs/(2*Fstop));   % upper bound for lowpass-safe decimation
```

## Step 1 — do not guess: confirm the intent

You must know:
- Mode: **streaming** (causal) vs **offline** (batch)
- Phase: **zero-phase**, **linear-phase**, or **don't-care**
- Whether internal **rate change** is acceptable (multirate pipelines)

If any is unknown: ask, then stop.

## Step 2 — pick the viable candidate families (don't force all four)

For `trans_pct < 2%`, present **2–4 viable candidates**, not a fixed set.

### Quick lookup: Mode × Phase → Viable families

| Mode | Phase | Single-stage IIR | Single-stage FIR | Multirate | IFIR |
|------|-------|------------------|------------------|-----------|------|
| **Streaming** | linear-phase | ❌ | ✅ | ✅ polyphase | ✅ |
| **Streaming** | don't-care | ✅ | ✅ | ✅ polyphase | ✅ |
| **Offline** | zero-phase | ✅ `filtfilt()` | ✅ `filtfilt()` | ✅ `resample()+filtfilt()` | ✅ `filtfilt()` |
| **Offline** | linear-phase | ❌ | ✅ | ✅ polyphase | ✅ |
| **Offline** | don't-care | ✅ | ✅ | ✅ either | ✅ |

**Key insight**: Polyphase multirate (FIR decimator/interpolator) gives **linear-phase** — it's NOT limited to offline zero-phase workflows.

### Family details

| Family | Viable when | Quick notes |
|---|---|---|
| Single-stage **IIR** | phase ≠ "linear-phase" | Offline can use `filtfilt()` for zero-phase. Streaming IIR is efficient but non-linear phase. |
| Single-stage **FIR** | always viable | Linear phase possible, but may be long / high latency / high MPIS. |
| **Multirate pipeline** (dec→filter→interp) | rate change OK | **Streaming**: polyphase System objects (linear-phase). **Offline zero-phase**: `resample()` + `filtfilt()`. Never mix `filtfilt` with System objects. |
| **Constant-rate multistage FIR** (IFIR method) | rate change NOT OK | FIR-like behavior at constant rate with fewer multipliers than single-stage FIR in many narrow-band cases. |

If a family is excluded, say so explicitly (one line), e.g.:
- "Excluded IIR: user requires linear phase in streaming."
- "Excluded multirate: user cannot change internal sample rate."

### Choosing between IFIR and Notch+LP

When the spectrum shows a **dominant tonal interferer** near the transition band, consider a **Notch + relaxed LP** hybrid:

| Approach | Best when |
|----------|-----------|
| **Notch + relaxed LP** | Known tonal interferer you can notch out (allows wider transition) |
| **IFIR** | General narrow transition, no specific tone to exploit |
| Single-stage FIR | Baseline (always works) |

**Decision rule:**
1. Analyze spectrum for dominant tones in/near transition band
2. If tone found: try Notch (at tone freq) + relaxed LP (wider Fstop)
3. Compare MPIS via `cost()` — pick the cheaper option
4. If no identifiable tone → use IFIR or single FIR

## Step 3 — feasibility check (order estimates)

### Default "MATLAB-native" flow (recommended)

Let `designfilt` pick the minimum order from `Rp/Rs`, then query `filtord(d)`.

```matlab
d_try = designfilt("lowpassfir", ...
    PassbandFrequency=Fpass, StopbandFrequency=Fstop, ...
    PassbandRipple=Rp, StopbandAttenuation=Rs, ...
    SampleRate=Fs, DesignMethod="equiripple");

N = filtord(d_try);
```

### FIR planning estimates (useful for architecture decisions)

```matlab
dev_p = (10^(Rp/20)-1)/(10^(Rp/20)+1);  % passband deviation (linear)
dev_s = 10^(-Rs/20);                    % stopband deviation (linear)

[Nk,~,~,~] = kaiserord([Fpass Fstop],[1 0],[dev_p dev_s], Fs);
[Np,~,~,~] = firpmord([Fpass Fstop],[1 0],[dev_p dev_s], Fs);
```

Use these as "order smell tests," not absolute truth.

## Step 4 — compute MPIS with `cost()`

**MPIS = Multiplications Per Input Sample**, reported by `cost()` on DSP System objects.

### 80/20 rules

- Tap count alone lies (multirate runs parts at lower rate; IIR cost isn't "sections × taps").
- Prefer **System objects** + `cost()`:
  - FIR → `dsp.FIRFilter`
  - IIR → `dsp.SOSFilter` (often easiest via `designfilt(..., SystemObject=true)`)
  - Pipelines → `dsp.FilterCascade`

### FIR MPIS

```matlab
d_fir = designfilt("lowpassfir", ...
    PassbandFrequency=Fpass, StopbandFrequency=Fstop, ...
    PassbandRipple=Rp, StopbandAttenuation=Rs, ...
    SampleRate=Fs, DesignMethod="equiripple");

firSys = dsp.FIRFilter("Numerator", d_fir.Numerator);
mpis_fir = cost(firSys).MultiplicationsPerInputSample;
```

### IIR MPIS (prefer SystemObject=true)

```matlab
% Returns dsp.SOSFilter directly (stable + stateful)
iirSys = designfilt("lowpassiir", ...
    PassbandFrequency=Fpass, StopbandFrequency=Fstop, ...
    PassbandRipple=Rp, StopbandAttenuation=Rs, ...
    SampleRate=Fs, DesignMethod="ellip", ...
    SystemObject=true);

mpis_iir = cost(iirSys).MultiplicationsPerInputSample;
```

### Pipeline MPIS

```matlab
pipe = dsp.FilterCascade(stage1, stage2, stage3);
mpis_pipe = cost(pipe).MultiplicationsPerInputSample;
```

For multirate decimate→filter→interpolate pipelines, `cost(pipe)` accounts for polyphase structure and internal rates.

### Offline `filtfilt()` note

- `filtfilt()` runs the filter **twice** → compute cost is roughly **2×** the single-pass cost.
- The effective magnitude response is **squared** (≈ doubles attenuation in dB).

## Step 5 — compare with Filter Analyzer

- Open `knowledge/cards/filter-analyzer.md` and overlay candidates.

Minimum comparisons to show:
- Magnitude response (meet specs)
- Group delay (latency + phase behavior)
- MPIS (when compute matters)

---

## `resample()` note

`resample()` uses internal filters not exposed to `cost()`.
If you need performance numbers for an offline `resample()` pipeline, use `timeit()` on your machine.

```

### knowledge/cards/filter-analyzer.md

```markdown
# Filter Analyzer Card

Open this card **before** you write any code that uses:
`filterAnalyzer`, `addFilters`, `addDisplays`, `replaceFilters`, `newSession`, or `filterAnalysisOptions`.

## Non‑negotiables

- Always pass `SampleRates=Fs` so plots are in **Hz** (not normalized).
- `FilterNames` must be valid MATLAB identifiers:
  - ✅ `["FIR_Equiripple","IIR_Ellip"]`
  - ❌ `["FIR Equiripple","IIR-ellip"]`
- Overlay rules:
  - frequency-domain overlays only with frequency-domain
  - time-domain overlays only with time-domain
- Avoid duplicate filters on re-runs:
  - prefer `newSession(fa)` (fresh start), or
  - `replaceFilters(fa, ...)` (update-in-place)

## Canonical “compare two filters” pattern

```matlab
Fs = 44100;

% d1, d2 can be digitalFilter objects or supported System objects
names = ["FIR_Equiripple","IIR_Ellip"];

[fa, dispMag] = filterAnalyzer(d1, d2, ...
    FilterNames=names, ...
    SampleRates=Fs, ...
    Analysis="magnitude", ...
    OverlayAnalysis="phase");

dispGD = addDisplays(fa, Analysis="groupdelay");
showFilters(fa, true, FilterNames=names, DisplayNums=dispGD);

dispImp = addDisplays(fa, Analysis="impulse");
showFilters(fa, true, FilterNames=names, DisplayNums=dispImp);
```

## Robust session pattern (no duplicates across reruns)

```matlab
Fs = 44100;
names = ["LP_A","LP_B"];

if exist("fa","var") && isvalid(fa)
    newSession(fa);
else
    fa = filterAnalyzer();  % open empty app
end

addFilters(fa, dA, dB, FilterNames=names, SampleRates=Fs);
addDisplays(fa, Analysis="magnitude", OverlayAnalysis="phase");
addDisplays(fa, Analysis="groupdelay");
```

## Update a filter without duplicating it

```matlab
% Replace existing filters by name (keeps displays intact)
replaceFilters(fa, dA_new, dB_new, ...
    FilterNames=["LP_A","LP_B"], ...
    SampleRates=Fs);
```

## Visualize a multirate pipeline as one “filter”

When comparing a decimate→filter→interpolate pipeline against a single-stage filter, wrap the full chain in `dsp.FilterCascade`.

```matlab
Fs = 44100; M = 4;

% Streaming-friendly components (System objects)
dec_sys = designMultirateFIR(DecimationFactor=M, StopbandAttenuation=60, SystemObject=true);

d_core = designfilt("lowpassfir", ...
    PassbandFrequency=2800, StopbandFrequency=3100, ...
    PassbandRipple=0.1, StopbandAttenuation=60, ...
    SampleRate=Fs/M, DesignMethod="equiripple");

core_sys = dsp.FIRFilter("Numerator", d_core.Numerator);

interp_sys = designMultirateFIR(InterpolationFactor=M, StopbandAttenuation=60, SystemObject=true);

pipe = dsp.FilterCascade(dec_sys, core_sys, interp_sys);

% Add alongside other filters
addFilters(fa, pipe, FilterNames="Multirate_Pipeline", SampleRates=Fs);
```

## When to open the full guide

Open `knowledge/filter-analyzer.md` when you need:
- `filterAnalysisOptions` (log scale, NFFT, normalization, etc.)
- display management (`duplicateDisplays`, `deleteDisplays`, …)
- session save/load (`saveSession`)
- version-specific notes

```