Chapter 2 described an LTI system by its impulse response $h[n]$ and the convolution sum. A second — and computationally far cheaper — description is the linear constant-coefficient difference equation（線性常係數差分方程, LCCDE）: the present output is computed from a finite number of past outputs, the present input, and past inputs.
Solving for the present output gives the recursive computational form:
| Non-recursive（非遞迴） | Recursive（遞迴） | |
|---|---|---|
| Condition | $N = 0$: no feedback terms | $N \ge 1$: output feeds back |
| Form | $y[n] = \sum_{k=0}^{M} b_k\,x[n-k]$ | $y[n] = -\sum_{k=1}^{N} a_k y[n-k] + \sum_{k=0}^{M} b_k x[n-k]$ |
| Impulse response | Finite length $M+1$ → FIR | Generally infinite → IIR |
| Memory needed | Past inputs only | Past inputs and past outputs |
| Example | moving average | accumulator $y[n]=y[n-1]+x[n]$ |
A classic non-recursive example is the $L$-point moving average（移動平均）, used everywhere for smoothing noisy data:
The direct method（直接法）solves an LCCDE in the time domain by splitting the total solution into two parts:
看你要什麼。直接法（$y_c+y_p$）和遞迴前推是兩種目的不同的工具：
只要「有限幾個數值」$y[0], y[1], \dots, y[n_0]$（$n_0$ 不大）→ 遞迴前推最保險。只要給了初始條件 $y[-1],\dots,y[-N]$，把 LCCDE 解成遞迴計算式
$$y[n] = \frac{1}{a_0}\!\left(\sum_{k=0}^{M} b_k\,x[n-k] \;-\; \sum_{k=1}^{N} a_k\,y[n-k]\right)$$從 $n=0$ 一格一格代下去即可。它不用猜特解型式、不用解特徵根、不用聯立求 $C_k$ —— 少掉所有代數出錯點，對遞迴／IIR 系統一樣有效。唯一要小心：初始條件代對、index 不要錯位。
要「任意 $n$ 的封閉解 $y[n]$」或看漸近／穩定行為 → 只能用直接法。遞迴前推只給你一串數值，給不出 symbolic 通式，也看不出 $n\to\infty$ 的趨勢。這時老實走 $y_c[n]+y_p[n]$：特徵方程求根 → 定特解 → 用初始條件解 $C_k$。
實務建議：需要通式就解 $y_c+y_p$，解完再用遞迴前幾項驗算（把 $n=0,1,2$ 代回封閉解，對照遞迴結果是否一致）—— 這是最穩的雙保險。考試若只問某個具體且 $n$ 小的 $y[n_0]$，直接遞迴最快、最不易錯。
The complementary solution（齊次解 / 互補解） $y_c[n]$ solves the homogeneous equation — the LCCDE with the input forced to zero:
Procedure. Assume an exponential trial solution $y[n] = \lambda^{\,n}$. Substituting and dividing out $\lambda^{\,n-N}$ yields the characteristic polynomial（特徵多項式）:
1. 代入試解。齊次式為 $\sum_{k=0}^{N} a_k\,y[n-k]=0$。把 $y[n]=\lambda^{\,n}$ 代進去，於是 $y[n-k]=\lambda^{\,n-k}$：
$$\sum_{k=0}^{N} a_k\,\lambda^{\,n-k} \;=\; a_0\lambda^{\,n} + a_1\lambda^{\,n-1} + \cdots + a_N\lambda^{\,n-N} \;=\; 0$$2. 提出最低次幂。每一項都含 $\lambda^{\,n-N}$（最小的指數），把它提出來：
$$\lambda^{\,n-N}\underbrace{\left(a_0\lambda^{N} + a_1\lambda^{N-1} + \cdots + a_N\right)}_{\text{特徵多項式}} \;=\; 0$$3. 約掉 $\lambda^{\,n-N}$。$\lambda=0$ 是平凡解（對應 $y[n]\equiv0$，沒意義），所以找非零解時 $\lambda^{\,n-N}\neq 0$，可以整條除掉 —— 指數項 $n$ 就此消失，只剩一個純代數方程：
$$a_0\lambda^{N} + a_1\lambda^{N-1} + \cdots + a_N = 0$$為什麼可以約掉 $\lambda^{\,n-N}$？因為它是「跟 $n$ 有關、但恆不為零」的因子，兩邊同除非零數是合法的。三層理由：
(i) 前提就排除 $\lambda=0$。$\lambda=0$ 給出 $y[n]=0^{\,n}\equiv0$，是沒資訊的平凡解；任何有意義的自然模態都對應 $\lambda\neq0$。
(ii) $\lambda\neq0 \Rightarrow \lambda^{\,n-N}\neq0$（對每個 $n$）。於是第 2 步的式子是「非零 × 某常數 $=0$」：
$$\lambda^{\,n-N}\cdot\underbrace{\big(a_0\lambda^{N}+\cdots+a_N\big)}_{=\,C\ \text{（不含 }n\text{ 的常數）}}=0 \;\Longrightarrow\; C=0$$(iii) 那個多項式因子完全不含 $n$。$n$ 只出現在被提出來的 $\lambda^{\,n-N}$ 裡；括號那項對所有 $n$ 是同一個常數。既然前面的因子恆不為零，要式子對所有 $n$ 成立，就只能令這個常數 $=0$ —— 這正是「約掉 $\lambda^{\,n-N}$」的真正意思：非零因子不背鍋，責任全丟給那個決定 $\lambda$ 的多項式。
關鍵：這一除，把「對每個 $n$ 都要成立」的差分方程，縮成「只需解出 $\lambda$」的 $N$ 次多項式。再把首項係數正規化（$a_0=1$）就得到上面的特徵方程。它的 $N$ 個根 $\lambda_k$ 就是系統的自然模態／極點；每個根貢獻一個 $\lambda_k^{\,n}$ 分量。
With $N$ distinct roots $\lambda_1, \dots, \lambda_N$ (the natural modes / poles of the system):
把 LCCDE 寫成狀態遞迴 $\mathbf{s}[n+1] = A\,\mathbf{s}[n]$（$A$ 為伴隨矩陣 companion matrix），則 $A$ 的特徵值恰好就是特徵根 $\lambda_k$，而 $y[n]$ 完全由 $A^n$ 決定。上面各種模態的形狀，其實就是 $A^n$ 的樣子。
相異根 → 可對角化。$N$ 個相異根時 $A$ 可對角化，$A = P\Lambda P^{-1}$，於是
$$A^n = P\,\Lambda^n P^{-1}, \qquad \Lambda^n = \operatorname{diag}(\lambda_k^{\,n}).$$得到純粹的模態 $\lambda_k^{\,n}$，沒有任何 $n$ 的多項式前綴——對應上表「distinct roots」那列。
重根 → 缺陷矩陣 defective → Jordan block。某根重複 $m$ 次時，伴隨矩陣少了特徵向量（幾何重數 < 代數重數），無法對角化，只能化為 Jordan block $J = \lambda I + N$（$N$ 為冪零矩陣）。此時
$$J^n = \sum_{k=0}^{m-1}\binom{n}{k}\lambda^{\,n-k} N^{k},$$其中二項式係數 $\binom{n}{k}$ 是 $n$ 的多項式（次數最高 $m-1$）。這正是
$$\lambda^{\,n},\; n\lambda^{\,n},\; n^2\lambda^{\,n},\;\dots,\; n^{m-1}\lambda^{\,n}$$這一整組模態的來源——你手動補上的 $n^{m-1}$ 因子，本質上就是 Jordan block 的 $\binom{n}{k}$。
共振陷阱其實是隱藏的 Jordan block。把輸入 $\alpha^n$ 看成由一個自己的外生狀態（特徵值 $\alpha$）產生，再與系統疊成一個更大的自治矩陣 $\tilde A$。當 $\alpha = \lambda_k$ 時，這個特徵值的代數重數比特徵向量多一個 → $\tilde A$ 變成缺陷矩陣，冒出一個 $2\times2$ Jordan block，其 $J^n$ 就產生 $n\,\alpha^{\,n}$ 項。所以「試解乘上 $n$」不是憑空規則，而是補上 Jordan 結構所要求的廣義特徵向量 generalized eigenvector 模態。
連續時間完全對應：Jordan block 的重根給出 $t^k e^{\lambda t}$，對應離散的 $n^k \lambda^n$。
The particular solution（特解） $y_p[n]$ is any one solution of the full equation with the input present, valid for $n \ge 0$. Assume a form that mimics the input:
| Input $x[n]$ (for $n \ge 0$) | Assumed particular solution $y_p[n]$ |
|---|---|
| $A$ (constant, e.g. step) | $K$ |
| $A\,\alpha^{\,n}$ | $K\,\alpha^{\,n}$ |
| $A\,n^{p}$ | $K_0 + K_1 n + \cdots + K_p n^{p}$ |
| $A\cos\omega_0 n$ or $A\sin\omega_0 n$ | $K_1\cos\omega_0 n + K_2 \sin\omega_0 n$ |
| $A\,\alpha^{\,n}$ where $\alpha$ is a characteristic root | $K\,n\,\alpha^{\,n}$ (multiply by $n$) |
Substitute the assumed form into the LCCDE and solve for $K$ by matching coefficients.
考慮一階系統，輸入的底數 $\alpha = 3$ 剛好等於特徵根 $\lambda = 3$：
$$y[n] - 3\,y[n-1] = 3^{\,n}, \qquad n \ge 0.$$先看錯誤的猜測 $y_p[n] = K\,3^{\,n}$。代入左式：
$$K\,3^{\,n} - 3\!\cdot\! K\,3^{\,n-1} = K\,3^{\,n} - K\,3^{\,n} = 0 .$$但右式是 $3^{\,n}\neq 0$，於是得到矛盾 $0 = 3^{\,n}$——這就是「$0 = \text{某值}$」。原因是 $3^{\,n}$ 本身就是自然模態，已滿足齊次方程 $y[n]-3y[n-1]=0$，永遠打不出非零的右式。
正確做法：乘上 $n$，試 $y_p[n] = K\,n\,3^{\,n}$。代入：
$$K\,n\,3^{\,n} - 3\!\cdot\! K\,(n\!-\!1)\,3^{\,n-1} = K\,3^{\,n}\big[\,n-(n\!-\!1)\,\big] = K\,3^{\,n}.$$令其等於右式 $3^{\,n}$，得 $K = 1$，故特解為
$$\boxed{\,y_p[n] = n\,3^{\,n}\,}.$$多出來的因子 $n$ 讓試解「跳出」齊次解的空間，代入後才會殘留非零項來對應輸入。若特徵根為雙重根（$\lambda$ 重複兩次）且輸入又是同一個模態，則要乘 $n^2$；一般而言重合 $m$ 重就乘 $n^m$。這與微分方程中遇到重根／共振時乘 $t$、$t^2$ 的規則完全對應。
Solve, for $n \ge 0$:
$$y[n] + y[n-1] - 12\,y[n-2] = x[n], \qquad x[n] = 10\,u[n], \qquad y[-1] = -1,\;\; y[-2] = -\tfrac12 .$$filter/filtic simulate: a = [1 1 -12]; b = 1; zi = filtic(b,a,[-1 -0.5]); y = filter(b,a,10*ones(1,5),zi) returns 5 -7 77 -151 1085 — matching the closed form. Exercise 3 repeats this workflow on the homework system.
The same total response can be split a different way, by cause instead of by mathematical form:
考試常考：給了非零初始條件的系統不是 LTI（線性性與時不變性會被初始狀態破壞），但它的零狀態部分 $y_{zs}[n]$ 才是 LTI——唯有在初始條件全為零（系統靜止）時，系統才滿足 $y_{zs}[n] = x[n]*h[n]$ 這個卷積關係。
The impulse response is the zero-state response to $x[n] = \delta[n]$. For a recursive LCCDE, the simplest route is to iterate the recursion with rest initial conditions; for low orders the closed form follows from the natural modes, since for $n > M$ the input terms vanish and $h[n]$ must be a pure combination of modes:
$y[n] = a\,y[n-1] + x[n]$, initially at rest, $x[n]=\delta[n]$.
scipy.signal.lfilter. The step response is the running sum of $h[n]$ and settles to $\sum_n h[n] = \frac{1}{1-0.8} = 5$.（一階遞迴系統：脈衝響應呈幾何衰減，步階響應收斂到直流增益 5。）The structure of the LCCDE determines the impulse-response length:
| FIR — Finite Impulse Response（有限脈衝響應） | IIR — Infinite Impulse Response（無限脈衝響應） | |
|---|---|---|
| LCCDE | $N = 0$ (non-recursive) | $N \ge 1$ (recursive) |
| Impulse response | $h[n] = b_n$ for $0 \le n \le M$, else $0$ — length $M+1$ | combination of modes $\lambda_k^n$, lasts forever |
| Stability | always BIBO stable | stable iff all $|\lambda_k| < 1$ |
| Implementation | convolution = the LCCDE itself | must use recursion (convolution would never end) |
| Memory of the past | exactly $M$ samples | fading but unlimited |
A system is BIBO stable（有界輸入有界輸出穩定） if every bounded input produces a bounded output: $|x[n]| \le B_x < \infty$ for all $n$ implies $|y[n]| \le B_y < \infty$ for all $n$.
For an LTI system, $|y[n]| = \big|\sum_k h[k]x[n-k]\big| \le B_x \sum_k |h[k]|$, which yields the complete characterization:
The condition is also necessary: if $S = \infty$, the bounded "worst-case" input $x[n] = \operatorname{sgn}\!\big(h[-n]\big)$ makes $y[0] = \sum_k |h[k]|$ blow up.
這一整章其實就是數位控制器的數學本體。對應關係：
1. LCCDE ＝ 數位控制器本身。數位控制器 $D(z)$ 落到硬體（MCU / DSP / FPGA）上執行時，就是一條差分方程的遞迴計算式：
$$u[n] = -\sum_{k=1}^{N} a_k\,u[n-k] + \sum_{k=0}^{M} b_k\,e[n-k]$$每個取樣時刻拿過去的 $u$、$e$ 算出當前輸出 $u[n]$ —— 正是 §3.1 的 recursive computational form。數位 PID（position / velocity form）就是這個。
2. 特徵根 $\lambda_k$ ＝ z 平面的極點。本節「特徵根判穩」的 $\lambda_k$，換到第 5 章 z-transform 語言就寫成極點 $z_k$。連續系統看「極點在左半平面（$\operatorname{Re}<0$）」；離散／數位系統看「極點在單位圓內（$|\lambda_k|<1$）」 —— 兩者由 $z=e^{sT}$ 對映。這裡的 BIBO 穩定條件，就是數位控制的閉迴路穩定條件。
用 s domain 一起看：模態的取樣就是 $z=e^{sT}$ 的由來。
連續系統的自然模態是 $e^{s_k t}$（$s_k$ 是 Laplace／s 平面的極點）。以週期 $T$ 取樣，令 $t=nT$：
$$e^{s_k t}\Big|_{t=nT} = \big(e^{s_k T}\big)^{n} = \lambda_k^{\,n}, \qquad \boxed{\;\lambda_k = e^{s_k T}\;}$$所以本章的離散模態 $\lambda_k^{\,n}$ 就是連續模態 $e^{s_k t}$ 的取樣版，$s$ 平面極點與 $z$ 平面極點由 $z=e^{sT}$ 一一對映。穩定區域也跟著搬：
| s 平面（連續） | $z=e^{sT}$ → z 平面（數位） | 模態行為 |
|---|---|---|
| 左半平面 $\operatorname{Re}(s)<0$ | 單位圓內 $|\lambda|<1$ | 衰減 → 穩定 |
| 虛軸 $\operatorname{Re}(s)=0$ | 單位圓上 $|\lambda|=1$ | 臨界 / 等幅振盪 |
| 右半平面 $\operatorname{Re}(s)>0$ | 單位圓外 $|\lambda|>1$ | 發散 → 不穩定 |
因此「s 平面看左半平面、z 平面看單位圓內」講的是同一件事 —— $e^{sT}$ 把「$\operatorname{Re}(s)<0$」這個左半平面直接彎成單位圓的內部。附帶一提：$s$ 的虛部（頻率 $\omega$）對映成 $\lambda$ 的輻角 $\omega T$，而 $\omega$ 每增加 $2\pi/T$ 就繞單位圓一圈回到原點 —— 這個「多對一」正是取樣的頻疊（aliasing），也是 z domain 只需看主圓、不需看整條 s 軸的原因。
3. 齊次解／特解 ＝ 暫態／穩態。$y_c[n]$（自然模態 $\lambda_k^{\,n}$）＝ transient response，$y_p[n]$＝ steady-state。控制器設計的極點配置（pole placement）本質就是「選 $\lambda_k$ 擺哪裡」，決定衰減快慢、有無振盪。
4. 閉迴路。受控體經 ZOH＋取樣得到 $G(z)$，配上控制器 $D(z)$，閉迴路特徵方程 $1+D(z)G(z)=0$ 的根＝閉迴路極點，決定整個系統響應。
一句話：這章教你「解一條差分方程 ＋ 判極點穩定」；數位控制課就是拿這套去「設計那條差分方程，把極點擺進單位圓」。
An LCCDE maps directly to hardware/software using three primitive elements（三種基本元件）:
Writing the second-order LCCDE in computational form,
and drawing one branch per term gives the Direct Form I（直接 I 型） realization — a feed-forward tapped delay line on the input plus a feedback tapped delay line on the output:
filter(b, a, x) uses the LCCDE written with all $y$ terms on the left: $\sum a_k y[n-k] = \sum b_k x[n-k]$. The multipliers in the feedback branches of the block diagram are therefore $-a_1, -a_2$ (negated!). Forgetting this sign flip is the single most common block-diagram error.
Correlation（相關） measures the degree of similarity between two sequences as a function of relative time shift. It is the standard tool for detecting a known waveform inside noise (radar, sonar, GPS), estimating time delays, and finding hidden periodicity.
The integer $\ell$ is the lag（延遲）: $r_{xy}[\ell]$ asks "how much does $x$ look like $y$ shifted right by $\ell$?" Reversing the roles only mirrors the lag axis:
Comparing with the convolution sum shows correlation is convolution with one sequence time-reversed — so all convolution machinery applies:
The autocorrelation（自相關） is the correlation of a signal with itself, $r_{xx}[\ell] = \sum_n x[n]\,x[n-\ell]$. Its key properties:
and for the cross-correlation the Cauchy–Schwarz bound gives $\big|r_{xy}[\ell]\big| \le \sqrt{r_{xx}[0]\,r_{yy}[0]} = \sqrt{E_x E_y}$. Dividing by the bound yields the normalized correlation（正規化相關係數）, handy because it is scale-free:
Periodic (power) signals have infinite energy, so the sum is replaced by a time average. For $x[n]$ periodic with period $N$:
The canonical application is radar/sonar delay estimation（雷達延遲估測）: transmit $x[n]$, receive an attenuated, delayed, noisy echo
Cross-correlating receive with transmit gives $r_{yx}[\ell] = \alpha\, r_{xx}[\ell - D] + r_{wx}[\ell]$. Since $r_{xx}$ peaks at zero lag and the noise term is small, $r_{yx}[\ell]$ peaks at $\ell = D$ — read off the round-trip delay, hence target range. The same idea underlies GPS code acquisition and pattern matching.
All four problems of ICE503 DSP Homework #3, fully worked.
Problem: For each of the following systems, determine whether the system is stable.
(a) $y[n] = \cos\big(x[n]\big)$ (b) $y[n] = r^{\,n} x[n],\; r > 1$ (c) $y[n] = n\,x[n]$ (d) $y[n] = \dfrac{1}{n}\,x[n]$
BIBO test: assume $|x[n]| \le B_x < \infty$ for all $n$ and check whether $|y[n]|$ must stay bounded. These are not LTI systems, so we argue directly from the definition (the $\sum|h|$ test does not apply).
(a) $y[n] = \cos(x[n])$ — STABLE. The cosine of any real number lies in $[-1,1]$, so for every input (bounded or not!) $|y[n]| \le 1$. Bounded output guaranteed → stable.
(b) $y[n] = r^n x[n]$, $r>1$ — UNSTABLE. Counterexample: the bounded input $x[n] = u[n]$ ($B_x = 1$) gives $y[n] = r^n$ for $n \ge 0$, and since $r > 1$, $r^n \to \infty$. No finite $B_y$ exists → unstable. (The time-varying gain $r^n$ is itself unbounded.)
(c) $y[n] = n\,x[n]$ — UNSTABLE. Counterexample: $x[n] = u[n]$ gives $y[n] = n \to \infty$. The gain $|n|$ grows without bound, so a constant input is amplified without limit → unstable.
(d) $y[n] = \frac{1}{n} x[n]$ — STABLE. For every $n \ne 0$, $\big|\tfrac{1}{n}\big| \le 1$, hence $$|y[n]| = \frac{|x[n]|}{|n|} \le |x[n]| \le B_x .$$ The output bound $B_y = B_x$ works for all bounded inputs → stable. (At $n = 0$ the rule $1/n$ is undefined; the system is understood as defined for $n \ne 0$, or with $y[0] := 0$ — either way the bound holds. Contrast with (c): here the gain shrinks with $|n|$ instead of growing.)
Problem: An LTI system is described by $$y[n] + 2y[n-1] - 8y[n-2] = 80\,x[n]$$ with initial conditions $y[-1] = -8$, $y[-2] = -9$. If the input is $x[n] = u[n]$, determine $y[n]$ for $n \ge 0$.
Step 1 — complementary solution. Characteristic equation from $y[n]+2y[n-1]-8y[n-2]=0$ with trial $y=\lambda^n$:
$$\lambda^2 + 2\lambda - 8 = 0 \;\Rightarrow\; (\lambda + 4)(\lambda - 2) = 0 \;\Rightarrow\; \lambda_1 = -4,\;\; \lambda_2 = 2$$so $y_c[n] = C_1(-4)^n + C_2\,2^n$.
Step 2 — particular solution. For $n \ge 0$ the input is the constant $x[n] = 1$, so try $y_p[n] = K$:
$$K + 2K - 8K = 80 \;\Rightarrow\; -5K = 80 \;\Rightarrow\; K = -16$$Step 3 — total solution form.
$$y[n] = C_1(-4)^n + C_2\,2^n - 16, \qquad n \ge 0$$Step 4 — get $y[0], y[1]$ from the recursion $y[n] = 80\,x[n] - 2y[n-1] + 8y[n-2]$:
$$\begin{aligned} y[0] &= 80(1) - 2y[-1] + 8y[-2] = 80 + 16 - 72 = 24\\ y[1] &= 80(1) - 2y[0] + 8y[-1] = 80 - 48 - 64 = -32 \end{aligned}$$Step 5 — solve for the constants.
$$\begin{aligned} n=0:&\quad C_1 + C_2 - 16 = 24 &&\Rightarrow\; C_1 + C_2 = 40\\ n=1:&\quad -4C_1 + 2C_2 - 16 = -32 &&\Rightarrow\; -2C_1 + C_2 = -8 \end{aligned}$$Subtracting: $3C_1 = 48 \Rightarrow C_1 = 16$, then $C_2 = 24$.
Verification by recursion: the formula gives $y[2] = 16(16)+24(4)-16 = 336$ and the recursion gives $y[2] = 80 - 2(-32) + 8(24) = 80+64+192 = 336$. ✓ Similarly $y[3] = -848$, $y[4] = 4464$ both ways.
Remark: the roots $-4$ and $2$ lie outside the unit circle, so this system is unstable — the response grows without bound even for the bounded step input, exactly as the numbers show.
Problem: MATLAB simulation of the system in Exercise 2, $y[n] + 2y[n-1] - 8y[n-2] = 80x[n]$, $x[n]=u[n]$, $y[-1]=-8$, $y[-2]=-9$:
(a) Implement the LCCDE (closed form) and determine $y[n]$ for $0 \le n \le 4$.
(b) Use a for loop to implement the recursion and determine $y[n]$ for $0 \le n \le 4$.
(c) Use filtic and filter to determine $y[n]$ for $0 \le n \le 4$. All three results must agree.
All three methods must reproduce $y[n] = 16(-4)^n + 24\,(2)^n - 16$ from Exercise 2, i.e.
$$y[0..4] = \{\,24,\; -32,\; 336,\; -848,\; 4464\,\}.$$clc; clear;
% LCCDE: y[n] + 2y[n-1] - 8y[n-2] = 80 x[n]
a = [1, 2, -8]; % output-side coefficients [a0 a1 a2]
b = 80; % input-side coefficient
x = ones(1, 5); % x[n] = u[n] for 0 <= n <= 4
y_m1 = -8; % y[-1]
y_m2 = -9; % y[-2]
% (a) closed-form solution from Exercise 2
n = 0:4;
y_closed = 16*(-4).^n + 24*2.^n - 16;
disp('(a) closed form:'); disp(y_closed);
% (b) direct recursion with a for loop
y = zeros(1, 5);
for k = 1:5
if k == 1 % n = 0: use y[-1], y[-2]
y(k) = b*x(k) - a(2)*y_m1 - a(3)*y_m2;
elseif k == 2 % n = 1: use y[0], y[-1]
y(k) = b*x(k) - a(2)*y(k-1) - a(3)*y_m1;
else % n >= 2
y(k) = b*x(k) - a(2)*y(k-1) - a(3)*y(k-2);
end
end
disp('(b) for loop:'); disp(y);
% (c) filtic + filter
zi = filtic(b, a, [-8, -9]); % past outputs [y(-1), y(-2)]
y_filt = filter(b, a, x, zi);
disp('(c) filtic + filter:'); disp(y_filt);
Output (identical for all three):
(a) closed form: 24 -32 336 -848 4464 (b) for loop: 24 -32 336 -848 4464 (c) filtic + filter: 24 -32 336 -848 4464
How each method works:
filtic(b, a, [y(-1) y(-2)]) converts the past outputs into the equivalent internal delay-line state zi of the Direct-Form filter, and filter(b, a, x, zi) then runs the same recursion in optimized form. With no zi, filter assumes a system at rest — it would compute only the zero-state response.Result regenerated numerically (Python scipy.signal.lfilter + lfiltic, the exact analogues of filter/filtic):
filter/filtic, and the closed form $16(-4)^n + 24\cdot 2^n - 16$ coincide exactly. The alternating, exploding values betray the unstable roots $-4$ and $2$.（三種算法完全一致；數值正負交替且快速放大，反映不穩定的特徵根。）Problem: MATLAB simulation:
(a) Generate the sinusoidal signal
$$x[n] = \begin{cases} \cos\!\big(\tfrac{1}{5}\pi n\big), & -5 \le n \le 5\\[2pt] 0, & -10 \le n < -5 \;\text{ and }\; 5 < n \le 10 \end{cases}$$
(b) Use the stem function to plot the autocorrelation of $x[n]$.
The autocorrelation of the finite-energy sequence is $r_{xx}[\ell] = \sum_n x[n]\,x[n-\ell]$, computed for lags $-20 \le \ell \le 20$ (a length-21 sequence has $2\cdot 21 - 1 = 41$ correlation lags).
clc; clear;
% (a) generate the windowed cosine on -10 <= n <= 10
n = -10:10;
x = zeros(1, length(n));
for i = 1:length(n)
if abs(n(i)) <= 5
x(i) = cos((1/5)*pi*n(i));
end % else: stays 0
end
% (b) autocorrelation r_xx[l] = sum_n x[n] x[n-l]
lags = -(length(n)-1) : (length(n)-1);
r = zeros(1, length(lags));
for k = 1:length(lags)
s = 0;
for i = 1:length(x)
j = i + lags(k); % index of x[n - l] (shifted copy)
if j >= 1 && j <= length(x) % keep index inside the array
s = s + x(i) * x(j);
end
end
r(k) = s;
end
% one-line check: r = xcorr(x); gives the same 41 values
figure;
subplot(2,1,1);
stem(n, x, 'filled'); grid on;
title('x[n] = cos(\pi n/5), |n| \leq 5'); xlabel('n'); ylabel('x[n]');
subplot(2,1,2);
stem(lags, r, 'filled'); grid on;
title('Autocorrelation r_{xx}[l]'); xlabel('lag l'); ylabel('r_{xx}[l]');
Result regenerated numerically (Python numpy.correlate(x, x, 'full')):
Checks against the theory of §3.7.2: