《数字信号处理》六、LTI 系统的变换域（z 域与频域）分析

在数字信号处理（DSP）的世界里，时域分析（如卷积和差分方程）虽然直观，但在面对复杂的系统设计（如滤波器设计、信道均衡）时，往往显得力不从心。

本章将带我们跨越时域的限制，进入更广阔的变换域（z 域与频域）。在这里，我们将解锁系统的“基因密码”——系统函数，并深入剖析全通系统、最小相位系统以及在实际工程中至关重要的广义线性相位系统。

6.1 LTI 系统的系统函数

6.1.1 LTI 系统的 z 域表示

对于一个线性时不变（LTI）系统，时域输入 $x(n)$ 与输出 $y(n)$ 的关系由单位脉冲响应 $h(n)$ 唯一确定：

$$y(n)=x(n)\ast h(n)=\sum _{k=-\infty }^{\infty }h(k)x(n-k)$$

根据 z 变换的时域卷积定理，时域的卷积在 z 域转化为简单的乘积：

$$Y(z)=H(z)X(z)$$

其中，$H(z)$ 称为系统的系统函数（System Function），它是单位脉冲响应 $h(n)$ 的 z 变换：

$$H(z)=\sum _{n=-\infty }^{\infty }h(n)z^{-n}$$

在实际工程中，LTI 系统通常由常系数线性差分方程描述：

$$\sum _{k=0}^N{a}_{k}y(n-k)=\sum _{r=0}^M{b}_{r}x(n-r)(a_0=1)$$

对两边同时取 z 变换，并利用位移性质 $Z\{x(n-k)\}=z^{-k}X(z)$（假设初始条件为零）：

$$Y(z)\sum _{k=0}^N{a}_k{z}^{-k}=X(z)\sum _{r=0}^M{b}_r{z}^{-r}$$

从而得到有理系统函数：

$$H(z)=\frac{Y(z)}{X(z)}=\frac{{\sum }_{r=0}^M{b}_r{z}^{-r}}{{\sum }_{k=0}^N{a}_k{z}^{-k}}$$

对其进行因式分解，可以写成零极点形式：

$$H(z)=A\frac{{\prod }_{r=1}^M(1-z_r{z}^{-1})}{{\prod }_{k=1}^N(1-p_k{z}^{-1})}=A{z}^{N-M}\frac{{\prod }_{r=1}^M(z-z_r)}{{\prod }_{k=1}^N(z-p_k)}$$

• $z_r$：系统的零点（Zeros），使 $H(z)=0$ 的点。
• $p_k$：系统的极点（Poles），使 $H(z)\to \infty$ 的点。

6.1.2 因果性和稳定性

一个 LTI 系统的时域特性（因果、稳定）在 z 域完全由其收敛域（ROC, Region of Convergence）及零极点分布决定。

1. 因果性（Causality）

• 时域条件：系统在 $n<0$ 时，单位脉冲响应 $h(n)=0$（即系统为右边序列）。
• z 域条件：$H(z)$ 的收敛域形如 $r_{max}<|z|\le \infty$，即收敛域在某个圆之外，且包含 $\infty$ 点。如果是因果有理系统，则系统函数的分子阶数 $M$ 不能大于分母阶数 $N$（在不包含无穷远极点的前提下）。

2. 稳定性（Stability）

• 时域条件：绝对可和，即 ${\sum }_{n=-\infty }^{\infty }|h(n)|<\infty$。
• z 域条件：系统函数 $H(z)$ 的收敛域必须包含单位圆（$|z|=1$）。

3. 因果稳定系统的充要条件

若系统既是因果的又是稳定的，则其收敛域必须同时满足上述两条：包含单位圆且单位圆在最大极点圆之外。这意味着：

因果稳定系统的所有极点必须严格位于单位圆以内，即 $|p_k|<1$。
零点的位置则不受此限制，可以在单位圆内、圆上或圆外。

6.1.3 单位脉冲响应计算

计算单位脉冲响应 $h(n)$，本质上是求解 $H(z)$ 的逆 z 变换（IZT）。常用的方法包括留数法、幂级数展开法（长除法）和部分分式展开法。

若 $H(z)$ 具有不同的单极点 $p_k$，可将其展开为：

$$H(z)=\sum _{r=0}^{M-N}{B}_r{z}^{-r}+\sum _{k=1}^N\frac{A_k}{1-p_k{z}^{-1}}$$

一旦求出系数 $A_k$，结合因果性条件，即可查表得到时域表达式：

$$h(n)=\sum _{r=0}^{M-N}{B}_r\delta (n-r)+\sum _{k=1}^N{A}_k(p_k{)}^{n}u(n)$$

例：求解$H(z)$时域表达$h(n)$
$$H(z)=\frac{1-0.5z^{-1}}{1-1.3z^{-1}+0.42z^{-2}}$$
分母可以因式分解为：
$$1-1.3z^{-1}+0.42z^{-2}=(1-0.7z^{-1})(1-0.6z^{-1})$$
由此得到系统的两个单极点分别为：$p_1=0.7$ 和 $p_2=0.6$。
因为分子阶数 $M=1$ 小于分母阶数 $N=2$（属于真分式），所以我们不需要进行长除法，前面的整式部分为 0。根据单极点展开公式，直接设：
$$H(z)=\frac{1-0.5z^{-1}}{(1-0.7z^{-1})(1-0.6z^{-1})}=\frac{A_1}{1-0.7z^{-1}}+\frac{A_2}{1-0.6z^{-1}}$$
我们可以使用留数法（留数定理）来快速求解这两个常数。
– 求解 $A_1$：方程两边同乘以 $A_1$ 的分母 $(1-0.7z^{-1})$，然后令 $z^{-1}=\frac{1}{0.7}$（即 $1-0.7z^{-1}=0$）：
A1=H(z)⋅(1−0.7z−1) ∣z−1=10.7=1−0.5z−11−0.6z−1 ∣z−1=10.7=2A_1 = H(z) \\cdot (1 – 0.7z^{-1}) \\ \\Big|_{z^{-1} = \\frac{1}{0.7}} = \\frac{1 – 0.5z^{-1}}{1 – 0.6z^{-1}} \\ \\Big|_{z^{-1} = \\frac{1}{0.7}}=2A1​=H(z)⋅(1−0.7z−1) ​z−1=0.71​​=1−0.6z−11−0.5z−1​ ​z−1=0.71​​=2
– 求解 $A_2$：方程两边同乘以 $A_2$ 的分母 $(1-0.6z^{-1})$，然后令 $z^{-1}=\frac{1}{0.6}$：
A2=H(z)⋅(1−0.6z−1) ∣z−1=10.6=1−0.5z−11−0.7z−1 ∣z−1=10.6=−1A_2 = H(z) \\cdot (1 – 0.6z^{-1}) \\ \\Big|_{z^{-1} = \\frac{1}{0.6}} = \\frac{1 – 0.5z^{-1}}{1 – 0.7z^{-1}} \\ \\Big|_{z^{-1} = \\frac{1}{0.6}}=-1A2​=H(z)⋅(1−0.6z−1) ​z−1=0.61​​=1−0.7z−11−0.5z−1​ ​z−1=0.61​​=−1
将求得的常数系数代回原式，得到最简部分分式形式：
$$H(z)=\frac{2}{1-0.7z^{-1}}-\frac{1}{1-0.6z^{-1}}$$
根据因果系统的常规假设，收敛域为外域 $|z|>0.7$。我们根据常用的经典 z 变换对：
$${\mathcal{Z}}^{-1}\left\{\frac{1}{1-p_k{z}^{-1}}\right\}=(p_k{)}^{n}u(n)$$
利用线性性质，直接逐项查表求出时域的单位脉冲响应 $h(n)$：
$$h(n)=2\cdot (0.7{)}^{n}u(n)-1\cdot (0.6{)}^{n}u(n)$$

下面用 MATLAB 代码演示一个因果系统的零极点分布，并计算、绘制其单位脉冲响应。

% ————————————————————————-
% Example: Analysis of Poles/Zeros and Impulse Response for a Causal LTI System
% H(z) = (1 – 0.5*z^-1) / (1 – 1.3*z^-1 + 0.42*z^-2)
% ————————————————————————-
clear; clc; close all;

% Define numerator and denominator coefficients
b = [1, –0.5];
a = [1, –1.3, 0.42];

% 1. Calculate zeros and poles
zeros_mat = roots(b);
poles_mat = roots(a);

fprintf('System Zero (z) locations: \\n'); disp(zeros_mat);
fprintf('System Pole (p) locations: \\n'); disp(poles_mat);

% 2. Plot Pole-Zero Map
figure('Position', [100, 100, 900, 400]);
subplot(1, 2, 1);
zplane(b, a);
grid on;
title('Pole-Zero Map in z-Domain');
legend('Zeros (z)', 'Poles (p)');

% 3. Compute and plot Impulse Response h(n)
N_samples = 30;
[h, t] = impz(b, a, N_samples);

subplot(1, 2, 2);
stem(t, h, 'filled', 'LineWidth', 1.5);
grid on;
xlabel('Time Index (n)');
ylabel('Amplitude h(n)');
title('Unit Impulse Response h(n)');

运行结果分析：
通过计算可以得到极点为 $p_1=0.7$, $p_2=0.6$。因为两个极点都绝对小于 1（在单位圆内部），所以该系统是因果且稳定的。从右侧的冲激响应图中也可以清楚地看到，$h(n)$ 随着 $n$ 的增大迅速衰减趋于 0。

6.2 LTI 系统的频域分析

6.2.1 频率响应的基本概念

如果系统是稳定的，其收敛域包含单位圆，那么我们就可以在单位圆上评估其系统函数。令 $z=e^{j\omega }$，此时的系统函数就蜕变为系统的频率响应（Frequency Response）：

H(ejω)=H(z)∣z=ejω=∑n=−∞∞h(n)e−jωnH(e^{j\\omega}) = H(z)\\Big|_{z=e^{j\\omega}} = \\sum_{n=-\\infty}^{\\infty} h(n)e^{-j\\omega n}H(ejω)=H(z)​z=ejω​=n=−∞∑∞​h(n)e−jωn

频率响应通常是一个复变函数，为了揭示其物理特性，我们常将其写成模与相角的形式（即极坐标形式）：

$$H(e^{j\omega })=|H(e^{j\omega })|e^{j\theta (\omega )}$$

• $|H(e^{j\omega })|$：幅度响应（Magnitude Response），表示系统对不同频率信号的放大或衰减倍数。
• $\theta (\omega )=\arg [H(e^{j\omega })]$：相位响应（Phase Response），表示不同频率信号通过系统后的相位滞后。

此外，为了评估信号通过系统后的时延失真，引入了群延时（Group Delay）的概念，定义为相位响应导数的负值：

$${\tau }_g(\omega )=-\frac{d\theta (\omega )}{d\omega }$$

例：假设一个数字滤波器的差分方程为：
$$y(n)=ay(n-1)+x(n)(0<a<1)$$
对差分方程两边取 $z$ 变换：
$$Y(z)=a{z}^{-1}Y(z)+X(z)⟹Y(z)(1-a{z}^{-1})=X(z)$$
得到系统函数为：
$$H(z)=\frac{1}{1-a{z}^{-1}}$$
该系统有一个极点 $p_1=a$。因为 $0<a<1$，极点在单位圆内，系统因果且稳定，收敛域（ROC）为 $|z|>a$，包含单位圆。因此，我们可以通过令 $z=e^{j\omega }$ 来计算它的频率响应。
将 $z=e^{j\omega }$ 代入系统函数：
$$H(e^{j\omega })=\frac{1}{1-a{e}^{-j\omega }}$$
利用欧拉公式 $e^{-j\omega }=\cos \omega -j\sin \omega$ 展开分母：
$$H(e^{j\omega })=\frac{1}{(1-a\cos \omega )+ja\sin \omega }$$
为了分别提取模（幅度）和相角（相位），我们需要对复数分母进行有理化（分子分母同乘以分母的共轭复数）：
$$H(e^{j\omega })=\frac{(1-a\cos \omega )-ja\sin \omega }{[(1-a\cos \omega )+ja\sin \omega ][(1-a\cos \omega )-ja\sin \omega ]}$$
$$H(e^{j\omega })=\frac{(1-a\cos \omega )-ja\sin \omega }{(1-a\cos \omega {)}^2+(a\sin \omega {)}^2}=\frac{(1-a\cos \omega )-ja\sin \omega }{1-2a\cos \omega +a^2}$$
现在，复数已被清晰地拆分为实部和虚部：
– 实部 $\text{Re}\{H(e^{j\omega })\}=\frac{1-a\cos \omega }{1-2a\cos \omega +a^2}$
– 虚部 $\text{Im}\{H(e^{j\omega })\}=\frac{-a\sin \omega }{1-2a\cos \omega +a^2}$
幅度响应是频率响应的模长，等于实部平方与虚部平方和的平方根。我们也可以直接从有理化前的分式直接求模：
$$|H(e^{j\omega })|=\frac{|1|}{|(1-a\cos \omega )+ja\sin \omega |}$$
$$|H(e^{j\omega })|=\frac{1}{\sqrt{(1-a\cos \omega {)}^2+(a\sin \omega {)}^2}}=\frac{1}{\sqrt{1-2a\cos \omega +a^2}}$$
物理意义分析：
– 当 $\omega =0$（直流/低频）时，$\cos (0)=1$，幅度为 $|H(e^{j0})|=\frac{1}{1-a}$。因为 $0<a<1$，这是一个较大的增益。
– 当 $\omega =\pi$（高频）时，$\cos (\pi )=-1$，幅度为 $|H(e^{j\pi })|=\frac{1}{1+a}$。这是一个较小的增益。
– 这表明该系统是一个典型的低通滤波器。
相位响应是复数的辐角，定义为 $\theta (\omega )=\arctan \left(\frac{\text{虚部}}{\text{实部}}\right)$：
$$\theta (\omega )=\arctan \left(\frac{\frac{-a\sin \omega }{1-2a\cos \omega +a^2}}{\frac{1-a\cos \omega }{1-2a\cos \omega +a^2}}\right)=\arctan \left(\frac{-a\sin \omega }{1-a\cos \omega }\right)$$
因为 $\arctan (-x)=-\arctan (x)$，可以提取负号表示相位滞后：
$$\theta (\omega )=-\arctan \left(\frac{a\sin \omega }{1-a\cos \omega }\right)$$
根据定义，群延时是对相位响应求导的负值：${\tau }_g(\omega )=-\frac{d\theta (\omega )}{d\omega }$。
由于 $\theta (\omega )=-\arctan (f(\omega ))$，利用复合函数求导法则以及 $(\arctan x{)}^{\prime }=\frac{1}{1+x^2}$：
$${\tau }_g(\omega )=\frac{d}{d\omega }\left[\arctan \left(\frac{a\sin \omega }{1-a\cos \omega }\right)\right]$$
令 $u=\frac{a\sin \omega }{1-a\cos \omega }$，先求 $\frac{du}{d\omega }$（使用商的求导公式）：
$$\frac{du}{d\omega }=\frac{(a\cos \omega )(1-a\cos \omega )-(a\sin \omega )(a\sin \omega )}{(1-a\cos \omega {)}^2}=\frac{a\cos \omega -a^2{\cos }^2\omega -a^2{\sin }^2\omega }{(1-a\cos \omega {)}^2}=\frac{a\cos \omega -a^2}{(1-a\cos \omega {)}^2}$$
然后代入总导数公式：
$${\tau }_g(\omega )=\frac{1}{1+u^2}\cdot \frac{du}{d\omega }=\frac{1}{1+{\left(\frac{a\sin \omega }{1-a\cos \omega }\right)}^2}\cdot \frac{a\cos \omega -a^2}{(1-a\cos \omega {)}^2}$$
$${\tau }_g(\omega )=\frac{(1-a\cos \omega {)}^2}{(1-a\cos \omega {)}^2+a^2{\sin }^2\omega }\cdot \frac{a\cos \omega -a^2}{(1-a\cos \omega {)}^2}$$
$${\tau }_g(\omega )=\frac{a\cos \omega -a^2}{1-2a\cos \omega +a^2{\cos }^2\omega +a^2{\sin }^2\omega }=\frac{a\cos \omega -a^2}{1-2a\cos \omega +a^2}$$
物理意义分析：
在低频段（$\omega \to 0$），${\tau }_g(0)=\frac{a-a^2}{1-2a+a^2}=\frac{a(1-a)}{(1-a{)}^2}=\frac{a}{1-a}$。若 $a=0.9$，则低频群延时为 $9$ 个采样点。由于群延时随 $\omega$ 变化，说明该系统存在相位失真（时延非线性）。

完整的 MATLAB 仿真验证代码

% ————————————————————————-
% Example: Multi-Frequency Analysis with Clean Subplot Zoom-ins
% H(z) = 1 / (1 – a*z^-1), with a = 0.5
% ————————————————————————-
clear; clc; close all;

a_coef = 0.5;
b = 1;
a = [1, –a_coef];

%% 1. Frequency Domain Analysis
w = linspace(0, pi, 512);
H_theory = 1 ./ (1 – a_coef * exp(–1i*w));
magnitude_theory = abs(H_theory);
phase_theory = –atan2(a_coef*sin(w), 1 – a_coef*cos(w));
gd_theory = (a_coef*cos(w) – a_coef^2) ./ (1 – 2*a_coef*cos(w) + a_coef^2);

%% 2. Time Domain Simulation with Multiple Frequencies
t_index = 0:150; % Discrete time index

w_low = 0.05 * pi;
w_med = 0.25 * pi;
w_high = 0.60 * pi;

envelope = exp(–((t_index–40)./15).^2);
x_low = envelope .* sin(w_low .* t_index);
x_med = envelope .* sin(w_med .* t_index);
x_high = envelope .* sin(w_high .* t_index);

y_low = filter(b, a, x_low);
y_med = filter(b, a, x_med);
y_high = filter(b, a, x_high);

gd_f = @(w_val) (a_coef*cos(w_val) – a_coef^2) / (1 – 2*a_coef*cos(w_val) + a_coef^2);

%% 3. Plotting 3-Row Comprehensive Visualization
figure('Position', [30, 30, 1400, 950]);

% ==========================================
% ROW 1: FREQUENCY DOMAIN CHARACTERISTICS
% ==========================================
subplot(3, 3, 1);
plot(w/pi, 20*log10(magnitude_theory), 'b', 'LineWidth', 2); grid on;
xlabel('Normalized Frequency (\\times\\pi rad/sample)'); ylabel('Magnitude (dB)');
title('Magnitude Response');

subplot(3, 3, 2);
plot(w/pi, phase_theory, 'g', 'LineWidth', 2); grid on;
xlabel('Normalized Frequency (\\times\\pi rad/sample)'); ylabel('Phase (rad)');
title('Phase Response');

subplot(3, 3, 3);
plot(w/pi, gd_theory, 'r', 'LineWidth', 2); grid on;
xlabel('Normalized Frequency (\\times\\pi rad/sample)'); ylabel('Group Delay (samples)');
title('Group Delay Response');

% ==========================================
% ROW 2: OVERVIEW OF TIME DOMAIN PULSES
% ==========================================
% Low Freq Overview
subplot(3, 3, 4);
plot(t_index, x_low, 'k–', 'LineWidth', 1.2); hold on;
plot(t_index, y_low, 'Color', [0.85 0.33 0.1], 'LineWidth', 2); grid on;
xlim([0, 120]); ylim([–1.2, 1.2]);
ylabel('Amplitude'); title('Low Freq Full View (\\omega = 0.05\\pi)');
legend('Input', 'Output', 'Location', 'northeast');
text(5, 0.9, sprintf('GD: %.3f', gd_f(w_low)), 'FontWeight', 'bold', 'BackgroundColor', [1 1 1 0.7]);

% Med Freq Overview
subplot(3, 3, 5);
plot(t_index, x_med, 'k–', 'LineWidth', 1.2); hold on;
plot(t_index, y_med, 'Color', [0.93 0.69 0.13], 'LineWidth', 2); grid on;
xlim([0, 120]); ylim([–1.2, 1.2]);
title('Medium Freq Full View (\\omega = 0.25\\pi)');
legend('Input', 'Output', 'Location', 'northeast');
text(5, 0.9, sprintf('GD: %.3f', gd_f(w_med)), 'FontWeight', 'bold', 'BackgroundColor', [1 1 1 0.7]);

% High Freq Overview
subplot(3, 3, 6);
plot(t_index, x_high, 'k–', 'LineWidth', 1.2); hold on;
plot(t_index, y_high, 'Color', [0 0.45 0.74], 'LineWidth', 2); grid on;
xlim([0, 120]); ylim([–1.2, 1.2]);
title('High Freq Full View (\\omega = 0.60\\pi)');
legend('Input', 'Output', 'Location', 'northeast');
text(5, 0.9, sprintf('GD: %.3f', gd_f(w_high)), 'FontWeight', 'bold', 'BackgroundColor', [1 1 1 0.7]);

% ==========================================
% ROW 3: ZOOMED-IN DETAILS (DELAY MARKS & TEXTS REMOVED)
% ==========================================
zoom_range = [32, 52]; % Focus precisely on the peak zone around sample 40

% 1. Low Freq Zoom
subplot(3, 3, 7);
plot(t_index, x_low, 'k–o', 'LineWidth', 1.2, 'MarkerIndices', 1:150); hold on;
plot(t_index, y_low, 'Color', [0.85 0.33 0.1], 'Marker', 'x', 'LineWidth', 2, 'MarkerIndices', 1:150);
grid on; xlim(zoom_range); ylabel('Amplitude'); xlabel('Time Index (n)');
title('Low Freq Zoom-in');
% All explicit peak vertical lines and delay text labels have been removed

% 2. Med Freq Zoom
subplot(3, 3, 8);
plot(t_index, x_med, 'k–o', 'LineWidth', 1.2, 'MarkerIndices', 1:150); hold on;
plot(t_index, y_med, 'Color', [0.93 0.69 0.13], 'Marker', 'x', 'LineWidth', 2, 'MarkerIndices', 1:150);
grid on; xlim(zoom_range); xlabel('Time Index (n)');
title('Medium Freq Zoom-in');
% All explicit peak vertical lines and delay text labels have been removed

% 3. High Freq Zoom
subplot(3, 3, 9);
plot(t_index, x_high, 'k–o', 'LineWidth', 1.2, 'MarkerIndices', 1:150); hold on;
plot(t_index, y_high, 'Color', [0 0.45 0.74], 'Marker', 'x', 'LineWidth', 2, 'MarkerIndices', 1:150);
grid on; xlim(zoom_range); xlabel('Time Index (n)');
title('High Freq Zoom-in');
% All explicit peak vertical lines and delay text labels have been removed

%% 4. Print Theoretical Analytical Delays to Command Window
fprintf('— Theoretical Group Delay Verification —\\n');
fprintf('Low Freq (0.05 pi) -> Delay: %.3f samples\\n', gd_f(w_low));
fprintf('Med Freq (0.25 pi) -> Delay: %.3f samples\\n', gd_f(w_med));
fprintf('High Freq (0.60 pi) -> Delay: %.3f samples\\n', gd_f(w_high));

该图完揭示了频域指标（频率响应）是如何直观映射到时域波形（延时与衰减）之中的：

6.2.2 LTI 系统的频率响应

利用几何方法，我们可以直观地从零极点向量看出系统的频率响应。将 $z=e^{j\omega }$ 代入零极点形式中：

$$H(e^{j\omega })=A{e}^{j\omega (N-M)}\frac{{\prod }_{r=1}^M(e^{j\omega }-z_r)}{{\prod }_{k=1}^N(e^{j\omega }-p_k)}$$

在 z 平面上，$e^{j\omega }$ 代表单位圆上的一个点，从零点 $z_r$ 指向该点的向量为 $V_r=e^{j\omega }-z_r$，从极点 $p_k$ 指向该点的向量为 $U_k=e^{j\omega }-p_k$。因此：

• 幅频特性：与零点向量长度的乘积成正比，与极点向量长度的乘积成反比。当单位圆上的点接近某个极点时，分母变小，幅频响应出现峰值；接近某个零点时，分子变小，幅频响应出现谷值。

6.3 全通系统

6.3.1 相同幅度响应的系统

在某些应用场景中（如全通均衡器），我们需要一种系统：它对所有频率成分的能量既不放大也不缩小，只改变它们的相位。这种系统被称为全通系统（Allpass System）。

其频率响应的幅度恒为常数（通常归一化为 1）：

$$|H_{ap}(e^{j\omega })|=1\forall \omega$$

为了在 z 域满足这一特性，全通系统的零极点必须呈共轭倒数关系。一个一阶全通系统的基本形式为：

$$H_{ap}(z)=\frac{z^{-1}-p^{\ast }}{1-p{z}^{-1}}=\frac{1-p^{\ast }z}{z-p}$$

若系统是因果稳定的，其极点 $p$ 必须在单位圆内（$|p|<1$），则对应的零点为 $1/p^{\ast }$，必然落在单位圆之外。

对于一个 $N$ 阶有理全通系统，其形式可以写成多个一阶或二阶全通子系统的级联：

$$H_{ap}(z)=\prod _{k=1}^N\frac{z^{-1}-p_k^{\ast }}{1-p_k{z}^{-1}}$$

6.3.2 全通系统的主要性质

下面的代码构建了一个二阶全通滤波器，验证其幅频特性是否为水平直线。

% ————————————————————————-
% 示例：二阶全通系统分析
% 极点(p)设在 0.6 + 0.4j 和 0.6 – 0.4j
% ————————————————————————-
p = [0.6 + 0.4i; 0.6 – 0.4i]; % 极点 p
z = 1 ./ conj(p); % 零点 z（共轭倒数）

% 从零极点构建多项式系数
a = poly(p);
b = poly(z);
% 保持全通增益归一化：全通分子系数是分母系数的倒序
b_ap = fliplr(a);

% 计算频响
[H_ap, w_ap] = freqz(b_ap, a, 512);

% 绘图
figure('Position', [100, 100, 900, 400]);
subplot(1, 2, 1);
zplane(b_ap, a);
title('全通系统零极点配对图');

subplot(1, 2, 2);
plot(w_ap/pi, abs(H_ap), 'LineWidth', 2);
ylim([0, 2]); grid on;
xlabel('归一化频率 (\\times\\pi rad/sample)');
ylabel('幅度响应');
title('全通系统幅度特性（恒等于1）');

6.4 最小相位系统

6.4.1 原系统及其逆系统

对于一个 LTI 系统 $H(z)$，如果我们想在接收端完全消除它对信号的影响，就需要构建一个逆系统（Inverse System） $H_i(z)$，使得：

$$H(z)H_i(z)=1⟹H_i(z)=\frac{1}{H(z)}$$

显然，原系统的零点变成了逆系统的极点，原系统的极点变成了逆系统的零点。
如果要求原系统和逆系统同时满足因果稳定，则两者的极点都必须严格位于单位圆以内。这意味着：

原系统 $H(z)$ 的所有极点 $p$ 和零点 $z$ 都必须在单位圆以内。
满足这一条件的系统，我们称之为最小相位系统（Minimum-Phase System）。

反之，若系统内部含有单位圆外部的零点，则称为最大相位系统或混合相位系统（统称为非最小相位系统）。

6.4.2 最小相位和全通分解

任何一个具有实质性物理意义的因果稳定系统 $H(z)$，如果其零点位于单位圆外，我们都可以利用全通系统将这些圆外的零点“反射”回单位圆内部，从而将其分解为：

$$H(z)=H_{min}(z)\cdot H_{ap}(z)$$

• $H_{min}(z)$：最小相位系统（包含了原系统的所有极点，以及被映射回圆内的零点）。
• $H_{ap}(z)$：全通系统（负责纠正相位差）。

这一性质在信道均衡与盲解卷积中具有极其重大的应用价值：它允许我们把不满足逆系统因果稳定条件的常规系统，拆解成一个可逆的骨架（最小相位）加一个不改动能量谱的全通滤片。

6.4.3 最小相位系统性质

在所有具有相同幅度响应的系统集合中：

我们设计两个系统，它们的幅频响应完全一样，但一个零点在圆内（最小相位），一个零点在圆外（非最小相位）。

% ————————————————————————-
% 示例：最小相位系统与非最小相位系统对比
% H1(z) = 1 – 0.5*z^-1 (零点在 0.5, 圆内 -> 最小相位)
% H2(z) = -0.5 + z^-1 (零点在 2.0, 圆外 -> 非最小相位)
% ————————————————————————-
b1 = [1, –0.5]; a1 = 1;
b2 = [–0.5, 1]; a2 = 1;

[H1, w] = freqz(b1, a1, 512);
[H2, w] = freqz(b2, a2, 512);

figure('Position', [100, 100, 1000, 400]);

% 1. 幅频对比
subplot(1, 2, 1);
plot(w/pi, 20*log10(abs(H1)), 'b-', 'LineWidth', 2); hold on;
plot(w/pi, 20*log10(abs(H2)), 'r–', 'LineWidth', 1.5);
grid on;
legend('H1 (最小相位)', 'H2 (非最小相位)');
xlabel('归一化频率 (\\times\\pi rad/sample)');
ylabel('幅度响应 (dB)');
title('幅度响应对比（完全相同）');

% 2. 相位响应对比
subplot(1, 2, 2);
plot(w/pi, unwrap(angle(H1)), 'b-', 'LineWidth', 2); hold on;
plot(w/pi, unwrap(angle(H2)), 'r–', 'LineWidth', 1.5);
grid on;
legend('H1 (最小相位)', 'H2 (非最小相位)');
xlabel('归一化频率 (\\times\\pi rad/sample)');
ylabel('相位响应 (rad)');
title('相位响应对比（H1 的相位延迟更小）');

6.5 广义线性相位系统

在现代数字通信和图像处理中，线性相位（Linear Phase）是一个极为苛刻但关键的指标。若系统相位非线性，不同频率的波形成分在通过系统后会产生不同的时延，从而导致时域波形严重畸变（如图像边缘模糊、数据码间串扰）。

6.5.1 线性相位系统特点

所谓第一类严格线性相位，是指系统的相位响应满足：

$$\theta (\omega )=-\alpha \omega$$

这意味着群延时 ${\tau }_g(\omega )=\alpha$ 是一个常数。

而广义线性相位（Generalized Linear Phase）则允许在原点处存在一个固定的恒定相位偏置 $\beta$：

$$\theta (\omega )=\beta -\alpha \omega$$

为了使 FIR 滤波器满足广义线性相位条件，其时域的单位脉冲响应 $h(n)$ 必须展现出严格的对称性或反对称性。设 FIR 滤波器的长度为 $N$（时域范围 $0\le n\le N-1$），对称中心为 $\alpha =(N-1)/2$：

• 对称（第一类广义线性相位, $\beta =0$ 或 $\pi$）：

$$h(n)=h(N-1-n)$$

• 反对称（第二类广义线性相位, $\beta =\pi /2$ 或 $-\pi /2$）：

$$h(n)=-h(N-1-n)$$

证：一个长度为 $N$ 的 FIR 滤波器的频率响应（DTFT）定义为：
$$H(\omega )=\sum _{n=0}^{N-1}h(n)e^{-j\omega n}$$
现在，我们找到这个序列的时间几何中心，令 $\alpha =\frac{N-1}{2}$。注意，不管 $N$ 是奇数还是偶数，$\alpha$ 永远是这个序列的正中心。
我们强行从求和公式中提取出一个整体的延迟因子 $e^{-j\omega \alpha }$，将公式重写为：
$$H(\omega )=e^{-j\omega \alpha }\sum _{n=0}^{N-1}h(n)e^{-j\omega (n-\alpha )}$$
为了看清求和号内部发生了什么，我们利用欧拉公式把 $e^{-j\omega (n-\alpha )}$ 展开为实部和虚部：
$$H(\omega )=e^{-j\omega \alpha }\sum _{n=0}^{N-1}h(n)\left[\cos \left(\omega (n-\alpha )\right)-j\sin \left(\omega (n-\alpha )\right)\right]$$
1. 当 $h(n)$ 严格偶对称时
偶对称意味着关于中心点对称的两侧系数相等，即 $h(n)=h(N-1-n)$。
– 观察虚部（正弦项）：由于 $\sin (-x)=-\sin (x)$ 是奇函数，在关于 $\alpha$ 对称的位置上，$\sin \left(\omega (n-\alpha )\right)$ 大小相等、符号相反。乘以相等的 $h(n)$ 后，所有的虚部在累加时两两完美抵消，求和结果为 0！
– 观察实部（余弦项）：由于 $\cos (-x)=\cos (x)$ 是偶函数，对称位置的项符号相同，互相叠加。
最终，系统函数精简为：
$$H(\omega )=e^{-j\omega \alpha }\cdot \underset{\text{这是一个纯实数函数 }A(\omega )}{\sum _{n=0}^{N-1}h(n)\cos \left(\omega (n-\alpha )\right)}$$
此时，整个系统的总相位就是前面那个提取出来的因子：
$$\theta (\omega )=-\alpha \omega$$
对照直线方程，此时 $\beta =0$，相位与频率 $\omega$ 成完美严格的正比（线性相位）！
2. 当 $h(n)$ 严格奇对称时
奇对称意味着关于中心点对称的两侧系数相反，即 $h(n)=-h(N-1-n)$。
– 这一次，实部（余弦项）因为正负相反，在累加时两两完全抵消，求和结果为 0。
– 而虚部（正弦项）由于“负负得正”，反而保留了下来。
最终，系统函数精简为：
$$H(\omega )=e^{-j\omega \alpha }\cdot \left[-j\sum _{n=0}^{N-1}h(n)\sin \left(\omega (n-\alpha )\right)\right]$$
虚数单位 $-j$ 可以写成指数形式 $-j=e^{-j\frac{\pi }{2}}$。我们把它和前面的项合并：
$$H(\omega )=e^{j\left(-\frac{\pi }{2}-\omega \alpha \right)}\cdot \underset{\text{纯实数函数 }A(\omega )}{\sum _{n=0}^{N-1}h(n)\sin \left(\omega (n-\alpha )\right)}$$
此时，整个系统的总相位为：
$$\theta (\omega )=-\frac{\pi }{2}-\alpha \omega$$
对照直线方程，此时 $\beta =-\frac{\pi }{2}$，同样满足广义线性相位的条件！

6.5.2 四类线性相位系统

根据 FIR 滤波器长度 $N$ 的奇偶性以及 $h(n)$ 的对称性，我们将广义线性相位系统划分为著名的四种类型。不同类型的系统具有不同的频域约束，直接决定了它们能用于设计哪种类型的滤波器。

下面通过表格全面总结这四类系统：

6.5.3 线性相位系统零点分布

由于 $h(n)$ 具有对称或反对称性，我们可以证明其系统函数满足：

$$H(z)=\pm z^{-(N-1)}H(z^{-1})$$

这导致线性相位 FIR 系统的零点(z)分布具有极其独特的互为倒数且共轭的镜像对称规律。

因此，除非零点刚好落在单位圆上（此时零点关于单位圆成对出现）或落在实轴上，否则一般的线性相位系统零点都是四个一组成非对称镜像分布的。

下面的脚本将一举生成四种不同类型的线性相位 FIR 系统，并输出它们的时域对称图形及对应的零点分布，供读者直观对比。

% ————————————————————————-
% MATLAB Script: Visual Comparison of 4 Types of Linear-Phase FIR Filters
% Demonstrating Time-Domain Symmetry and Z-Plane Zero Locations
% ————————————————————————-
clear; clc; close all;

%% 1. Define Filter Coefficients for the 4 Types

% Type I: Odd Length (N=7), Symmetric
h1 = [1, 2, 3, 4, 3, 2, 1];

% Type II: Even Length (N=6), Symmetric
h2 = [1, 2, 3, 3, 2, 1];

% Type III: Odd Length (N=7), Anti-Symmetric (Center must be 0)
h3 = [1, 2, 3, 0, –3, –2, –1];

% Type IV: Even Length (N=6), Anti-Symmetric
h4 = [1, 2, 3, –3, –2, –1];

%% 2. Plotting Layout (4 Rows x 2 Columns)
figure('Position', [100, 50, 950, 900], 'Name', '4 Types of Linear Phase FIR Systems');

% =========================================================================
% ROW 1: TYPE I
% =========================================================================
subplot(4, 2, 1);
stem(0:length(h1)–1, h1, 'b', 'LineWidth', 2, 'MarkerFaceColor', 'b'); grid on;
title('Type I: Odd Length (N=7), Symmetric');
xlabel('Time Index (n)'); ylabel('h[n]');

subplot(4, 2, 2);
zplane(h1, 1); grid on;
title('Type I Z-Plane (No Constraints)');

% =========================================================================
% ROW 2: TYPE II
% =========================================================================
subplot(4, 2, 3);
stem(0:length(h2)–1, h2, 'r', 'LineWidth', 2, 'MarkerFaceColor', 'r'); grid on;
title('Type II: Even Length (N=6), Symmetric');
xlabel('Time Index (n)'); ylabel('h[n]');

subplot(4, 2, 4);
zplane(h2, 1); grid on;
title('Type II Z-Plane (Zero at \\omega = \\pi)');

% =========================================================================
% ROW 3: TYPE III
% =========================================================================
subplot(4, 2, 5);
stem(0:length(h3)–1, h3, 'g', 'LineWidth', 2, 'MarkerFaceColor', 'g'); grid on;
title('Type III: Odd Length (N=7), Anti-Symmetric');
xlabel('Time Index (n)'); ylabel('h[n]');

subplot(4, 2, 6);
zplane(h3, 1); grid on;
title('Type III Z-Plane (Zeros at \\omega = 0, \\pi)');

% =========================================================================
% ROW 4: TYPE IV
% =========================================================================
subplot(4, 2, 7);
stem(0:length(h4)–1, h4, 'm', 'LineWidth', 2, 'MarkerFaceColor', 'm'); grid on;
title('Type IV: Even Length (N=6), Anti-Symmetric');
xlabel('Time Index (n)'); ylabel('h[n]');

subplot(4, 2, 8);
zplane(h4, 1); grid on;
title('Type IV Z-Plane (Zero at \\omega = 0)');

% Adjust layout spacing
sgtitle('Comparison of 4 Types of Linear-Phase FIR Filters', 'FontSize', 14, 'FontWeight', 'bold');

关键技术细节复盘：

• 在观察 Type II 的零点图时，你会发现 $\omega =\pi$（即 $z=-1$ 处）有一个固定死掉的零点。这就解释了为什么它无法通过高频信号，从而绝对不能用来做高通滤波器。
• 在观察 Type III 和 IV 的零点图时，注意 $z=1$ （即 $\omega =0$ 处）都有一个零点，说明它们天生就具备压制直流成分（DC）的特性。

本章小结与工程启示

第六章的内容是串联整个数字信号处理中“分析”与“设计”的黄金桥梁。

• 我们用 $H(z)$ 的极点 $p$ 去框定系统的因果稳定性，确保系统在工程上可控、可用。
• 我们用 全通系统 作为相位校正的利器，在不改变信号能量分布的前提下微调相位。
• 我们用 最小相位系统 解决因果求逆的问题，为信道盲均衡和系统辨识提供了底层数学方案。
• 我们用 广义线性相位系统 指导 FIR 滤波器的宏观设计，通过时域的对称性换取频域的无延迟失真。

掌握了变换域的这四种方法，在接下来面对具体的 FIR 和 IIR 滤波器工程设计时，你将拥有一种高屋建瓴的全局视角。

*欢迎在评论区留下您的思考！