数字信号处理

• 杂七杂八
  • 卷积
    • 线性卷积
    • 圆周卷积
  • 模拟信号采样
  • 数学不好
• 信号变换
  • 傅里叶变换
    • 连续时间
      • CFS
      • CFT
    • 离散时间
      • DTFS
      • DTFT
    • 特点
  • z 变换
    • 双边
      • 定义
      • 常见 z 变换对
      • 性质
      • 逆变换
    • 单边
    • 与傅里叶变换的关系
• DFT & FFT
  • DFT & IDFT 的计算
  • 基2 FFT 算法
    • 按时间抽取
      • 蝶形运算
    • 按频率抽取
  • 基4 FFT 算法
  • N 为复合数的 FFT 算法
  • 调频 z 变换算法
• 离散时间系统
  • 系统函数
  • FIR 系统
    • 直接型结构
    • 级联型结构
  • IIR 系统
    • 直接型结构
      • 直接 I 型
      • 直接 II 型
    • 级联型结构
    • 并联型结构
• 滤波器设计
  • 数字滤波器
    • 技术指标
    • 频率响应
  • IIR 滤波器
    • 模拟滤波器
      • 巴特沃斯低通滤波器
      • 切比雪夫 I 型低通滤波器
      • 模拟频域频带变换
    • 间接法设计数字滤波器
      • 冲激响应不变法
      • 双线性变换法
      • 数字频域频带变换
  • FIR 滤波器
    • 滤波器特点
    • 窗函数设计法
      • 计算步骤
    • 频率抽样设计法
      • 基本思路
      • 设计公式
    • 与 IIR 比较

杂七杂八

卷积

线性卷积

\[ x(n) * h(n) = \sum\limits_{k = -\infty}^{\infty} x(k)h(n-k) \]

• 恒等性：\(x(n)*\delta(n) = x(n)\)
• 位移性：\(x(n)*\delta(n-k) = x(n-k)\)
• 交换律：\(x(n)*h(n) = h(n)*x(n)\)
• 结合律：\([x_1(n)*h_1(n)]*h_2(n) = x_1(n)*[h_1(n)*h_2(n)]\)
• 分配律：\(x(n)*[h_1(n) + h_2(n)] = x_1(n)*h_1(n) + x_1(n)*h_2(n)\)

圆周卷积

对于两个长度为 N 的序列 \(x(n)\) 和 \(h(n)\)：
\[ x(n) \otimes h(n) = \sum\limits_{k=0}^{N-1} x(k) [h(n-k) \mod N] \]
即将 \(h(n)\) 扩充为周期为 N 的序列。

可以通过 DFT 和 IDFT 的计算获得圆周卷积的结果： - 对输入序列 x(n) 和 h(n) 分别进行 0 填充，使二者长度均变为 2N-1。 - 对填充后的序列做 2N-1 点 DFT，得到 X(k) 和 H(k)。 - 对 X(k)H(k) 这个乘积做 IDFT，即可得到圆周卷积的结果。

模拟信号采样

\[ x(n) = x_a(nT) \]

• 模拟信号频率：\(-\infty < F < \infty,\ -\infty < \Omega < \infty\)

• 采样周期：\(T\)

• 采样频率：\(F_s = \dfrac{1}{T}\)

• 离散信号频率：\(-\dfrac{1}{2} < f = \dfrac{F}{F_s} < \dfrac{1}{2},\ -\pi < \omega = \Omega T < \pi\)

• 对某一采样率 \(F_s\)：\(F_{max} = \dfrac{F_s}{2} = \dfrac{1}{2T},\ \Omega_{max} = \pi F_s = \dfrac{\pi}{T}\)
  超过的部分会发生混叠。

数学不好

• \(\sin \theta = \dfrac{1}{2}(e^{j\theta} - e^{-j\theta}),\ \cos \theta = \dfrac{1}{2}(e^{j\theta} + e^{-j\theta})\)

信号变换

傅里叶变换

连续时间

连续时间信号具有非周期频谱。

CFS

连续时间周期信号的傅里叶极数。

对周期信号 \(x(t)\)：

变换	公式
逆变换	\(x(t) = \sum\limits_{k = -\infty}^{\infty} c_k e^{j 2\pi k F_0 t}\)
正变换	\(c_k = \dfrac{1}{T_p}\displaystyle \int_{T_p} x(t) e^{-j 2\pi k F_0 t} dt\)

其中，\(T_p = \dfrac{1}{F_0}\) 为信号 \(x(t)\) 的基本周期。

一个周期信号具有无限能量，但具有一个有限的平均功率。
功率信号的帕塞瓦关系式：
\[ P_x = \dfrac{1}{T_p} \displaystyle \int_{T_p} |x(t)|^2 dt = \sum\limits_{k = -\infty}^{\infty} |c_k|^2 \]

功率密度谱：以频率 \(F = kF_0\) 为横轴，\(|c_k|^2\) 为纵轴，偶对称。

\(c_k\) 幅度谱偶对称，相位谱奇对称。

CFT

连续时间非周期信号的傅里叶变换。

对无限时长的非周期信号 \(x(t)\)：

变换	公式（频率）	公式（角频率）
逆变换	\(x(t) = \displaystyle \int_{-\infty}^{\infty} X(F) e^{j 2\pi F t} dF\)	\(x(t) = \dfrac{1}{2\pi} \displaystyle\int_{-\infty}^{\infty} X(\Omega) e^{j \Omega t} d\Omega\)
正变换	\(X(F) = \displaystyle \int_{-\infty}^{\infty} x(t) e^{-j 2\pi F t} dt\)	\(X(\Omega) = \displaystyle\int_{-\infty}^{\infty} x(t) e^{-j \Omega t} dt\)

如果是有限时长的信号，修改积分上下限即可。

非周期有限能量信号的帕塞瓦关系式：
\[ E_x = \displaystyle \int_{-\infty}^{\infty} |x(t)|^2 dt = \displaystyle \int_{-\infty}^{\infty} |X(F)|^2 dF \]

能量密度谱：以频率 \(F\) 为横轴，\(S_{xx} = |X(F)|^2\) 为纵轴。
实信号的能量密度谱偶对称，频域幅度谱偶对称，相位谱奇对称。

离散时间

DTFS

离散时间周期信号的傅里叶极数。

对周期为 \(N\) 的序列 \(x(n)\)：

变换	公式
逆变换	\(x(n) = \sum\limits_{k=0}^{N-1} = c_k e^{j 2\pi k n/N}\)
正变换	\(c_k = \dfrac{1}{N} \sum\limits_{n=0}^{N-1} x(n) e^{-j 2\pi k n/N}\)

其中 \(k = 0, 1, \cdots, N-1\)。
如果扩展到该范围之外，则 \(c_{k+N} = c_k\)，即 \(c_k\) 是周期为 \(N\) 的序列。

离散周期信号的帕塞瓦关系式：
\[ P_x = \sum\limits_{k=0}^{N-1} |c_k|^2 = \dfrac{1}{N} \sum\limits_{n=0}^{N-1} |x(n)|^2 \]

与连续时间周期信号一样，频域幅度谱偶对称，相位谱奇对称。

DTFT

离散时间非周期信号的傅里叶变换。

对能量有限离散时间信号 \(x(n)\)：

变换	公式
逆变换	\(x(n) = \dfrac{1}{2\pi} \displaystyle\int_{-\pi}^{\pi} X(\omega) e^{j \omega n} d\omega\)
正变换	\(X(\omega) = \sum\limits_{n = -\infty}^{\infty} x(n) e^{-j \omega n}\)

记 \(X_N(\omega) = \sum\limits_{n = -N}^{N} x(n) e^{-j \omega n}\)
Gibbs 现象：在 \(X(\omega)\) 的一个不连续点处，\(X_N(\omega)\) 逼近 \(X(\omega)\) 的振荡行为。

有限能量的离散时间非周期信号的帕塞瓦关系式：
\[ E_x = \sum\limits_{n = -\infty}^{\infty} |x(n)|^2 = \dfrac{1}{2\pi} \displaystyle\int_{-\pi}^{\pi} |X(\Omega)|^2 \]

记 \(S_{xx} = |X(\omega)|^2\)，能量密度谱。
对于实信号，频域幅度谱偶对称，相位谱奇对称，能量密度谱偶对称。

特点

• 连续时间信号具有非周期频谱。
• 离散时间信号具有周期频谱。
• 周期信号具有离散频谱。
• 非周期有限能量信号具有连续频谱。

z 变换

双边

定义

\[ X(z) = \sum\limits_{n = -\infty}^{\infty} x(n) z^{-n} \]

收敛域： - 因果序列：\(|z| > r_2\) - 非因果序列：\(|z| < r_1\) - 双边序列：\(r_1 < |z| < r_2\)

常见 z 变换对

信号	Z变换结果	收敛域
\(\delta[n]\)	\(1\)	全平面
\(u[n]\)	\(\dfrac{1}{1 - z^{-1}}\)	\(\|z\| > 1\)
\(a^n u[n]\)	\(\dfrac{a z^{-1}}{(1 - az^{-1})^2}\)	\(\|z\| > \|a\|\)
\(n a^n u[n]\)	\(\dfrac{1}{1 - az^{-1}}\)	\(\|z\| > \|a\|\)
\(\sin(\omega_0 n) u(n)\)	\(\dfrac{z^{-1} \sin(\omega_0)}{1 - 2 z^{-1}\cos(\omega_0) + z^{-2}}\)	\(\|z\| > 1\)
\(\cos(\omega_0 n) u(n)\)	\(\dfrac{1 - z^{-1}\cos(\omega_0)}{1 - 2 z^{-1}\cos(\omega_0) + z^{-2}}\)	\(\|z\| > 1\)
\(e^{j\omega_0 n}\)	\(\dfrac{1}{1 - z^{-1}e^{j\omega_0}}\)	全平面

性质

性质	时域	z 域	收敛域
记号	\(x(n)\)	\(X(z)\)	收敛域0： \(r_2 < \|z\| < r_1\)
	\(x_1(n)\)	\(X_1(z)\)	收敛域1
	\(x_2(n)\)	\(X_2(z)\)	收敛域2
线性	\(a_1x_1[n] + a_2x_2[n]\)	\(a_1X_1(z) + a_2X_2(z)\)	至少是 1 和 2 的交集
时移	\(x[n-n_0]\)	\(X(z)z_0^{-n}\)	一般同 0
z 域尺度变换	\(a^nx[n]\)	\(X(a^{-1}z)\)	\(\|a\|r_2 < \|z\| < \|a\|r_1\)
时间翻转	\(x[-n]\)	\(X(z^{-1})\)	\(\dfrac{1}{r_1} < \|z\| <\dfrac{1}{r_2}\)
共轭	\(X^*(z^*)\)	\(x^*[n]\)	0
实部	\(Re\{x(n)\}\)	\(\dfrac{1}{2}[X(z) + X^*(z^*)]\)	包括 0
虚部	\(Im\{x(n)\}\)	\(\dfrac{1}{2} j [X(z) - X^*(z^*)]\)	包括 0
z 域求导	\(nx[n]\)	\(-z\dfrac{dX(z)}{dz}\)	0
时域积分	\(\sum\limits_{n=0}^{\infty} x[n]z^{-n}\)	\(\dfrac{X(z)}{1-z^{-1}}\)
卷积	\(x_1[n] * x_2[n]\)	\(X_1(z) \cdot X_2(z)\)	至少是 1 和 2 的交集
相关	\(r_{x_1 x_2} = x_1[l] * x_2[-l]\)	\(R_{x_1 x_2} = X_1(z)X_2(z)\)
初值定理	\(x[0]\)	\(\lim\limits_{z \to \infty}X(z)\)
终值定理	\(\lim\limits_{n \to \infty} x[n]\)	\(\lim\limits_{z \to 1}(z-1)X(z)\)
频域平移	\(e^{j\omega_0n}x[n]\)	\(X(ze^{j\omega_0})\)
频域缩放性质	\(\dfrac{1}{\|a\|}x\left [\dfrac{n}{a} \right]\)	\(X(az)\)

逆变换

• 围线积分法 \[ x(n) = \sum\limits_{z_i}\left[X(z)z^{n-1} 在 z=z_i 的留数\right] \]

  对于 m 阶极点 a，f(z) 的留数为： \[ \dfrac{1}{(m-1)!} \lim\limits_{z \rightarrow a}\left[(z - a)^m f(z)\right]^{(m-1)} \]

• 幂级数展开法
  \[ X(z) = \sum\limits_{n=-\infty}^{\infty} c_n z^{-n} \rightarrow x(n) = c_n \]

  采用长除法

  • 因果序列：高次幂在前
  • 非因果序列：低次幂在前

• 部分分式展开法

  • \(X(z)\) 有两个相异极点 \(p_1,\ p_2\)
    \[ X(z) = \dfrac{A_1}{z-p_1} + \dfrac{A_2}{z-p_2} \]

  • \(X(z)\) 有两个相同极点 \(p\)
    \[ X(z) = \dfrac{A_1}{z-p} + \dfrac{A_2}{(z-p)^2} \]

单边

\[ X^+(z) = \sum\limits_{n=0}^{\infty} x(n)z^{-n} \]

性质	时域	z 域
时移	\(x(n-k)\)	\(\left[x(-k) + x(-k+1)z^{-1} + \cdots +x(-1)z^{-k+1}\right] + z^{-k}X^+(z)\)
时间超前	\(x(n+k)\)	\(z^k \left[ X^+(z) - \sum\limits_{n=0}^{k-1} x(n)z^{-n}\right]\)
终值定理	\(\lim\limits_{n \rightarrow \infty} x(n)\)	\(\lim\limits_{z \rightarrow 1} (z-1)X^+(z)\)

与傅里叶变换的关系

\[ X(z=e^{j\omega}) = X(\omega) = \sum\limits_{n = -\infty}^{\infty} x(n) e^{-j\omega n} \]

DFT & FFT

DFT & IDFT 的计算

记 \[ W_N = e^{-j 2\pi / N} \] \(W_N\) 具有对称性和周期性，即 \[ W_N^{k + N/2} = -W_N^k \] \[ W_N^{k + N} = W_N^k \]

DFT：对长为 N 的数据序列 x(n) 计算长度为 N 的复序列 X(k)。
IDFT：DFT 的逆变换。

变换	公式
DFT	\(X(k) = \sum\limits_{n=0}^{N-1} x(n) W_N^{kn}\quad (0 \leq k \leq N-1)\)
IDFT	\(x(n) = \dfrac{1}{N} \sum\limits_{k=0}^{N-1} X(k) W_N^{-kn}\quad (0 \leq n \leq N-1)\)

对每个 k 值，直接计算 X(k) 涉及 N 次复数乘法（4N 次实数乘法）和 N-1 次复数加法（4N-2 次实数加法）。
因此计算长度为 N 的 DFT 要求 \(N^2\) 次复数乘法和 \(N^2 - N\) 次复数加法。

基2 FFT 算法

适用于 \(N = 2^v\)。

按时间抽取

将 \(x(n)\) 划分为奇序列和偶序列（以 2 为因子抽取），即 \[ f_1(n) = x(2n),\quad f_2(n) = x(2n + 1)\quad (n = 0, 1, \cdots, \frac{N}{2}-1) \]

二者的 \(\frac{N}{2}\) 点 DFT 结果分别为 \(F_1(k)\)、\(F_2(k)\)。

N 点 DFT 算法用抽取序列的 DFT 可表示为： \[ \begin{aligned} X(k) &= \sum\limits_{n=0}^{N-1} x(n) W_N^{kn},\quad k = 0,1, \cdots, N-1 \\ &= \sum\limits_{n为偶数} x(n) W_N^{kn} + \sum\limits_{n为奇数} x(n) W_N^{kn} \\ &= \sum\limits_{m=0}^{(N/2)-1} x(2m) W_N^{2mk} + \sum\limits_{m=0}^{(N/2)-1} x(2m+1) W_N^{(2m+1)k} \\ &= \sum\limits_{m=0}^{(N/2)-1} f_1(m) W_{N/2}^{mk} + W_N^k \sum\limits_{m=0}^{(N/2)-1} f_2(m) W_{N/2}^{mk} \\ &= F_1(k) + W_N^k F_2(k),\quad k = 0,1, \cdots, \frac{N}{2}-1 \\ X(k + \frac{N}{2}) &= F_1(k) - W_N^k F_2(k) \end{aligned} \]

直接计算 \(F_1(k)\) 需要 \((\frac{N}{2})^2\) 次复数乘法，直接计算 \(F_2(k)\) 同理，计算 \(W_N^k F_2(k)\) 还需要 \(\frac{N}{2}\) 次复数乘法。
所以，计算 \(X(k)\) 一共需要 \(\frac{N^2 + N}{2}\) 次复数乘法。

同样的，对 \(F_1(k)\)和\(F_2(k)\) 的计算也可以通过计算 \(\frac{N}{4}\) 点 DFT 来完成，从而只需要 \(\frac{N^2 + 2N}{4}\) 次复数乘法。
这时计算 \(X(k)\) 一共需要 \(\frac{N^2}{4} + N\) 次复数乘法。

数据序列的抽取可以不断重复，直到最后序列变为 1 点序列。
对于\(N = 2^v\)，抽取次数为 \(v\)，复数乘法总次数减少为 \(\frac{N}{2} \log_2(N)\)，复数加法的次数为 \(N \log_2(N)\)。

eg. 如果 \(x(n)\) 为 N 点实信号
从上述推导中可以得出：
\[ \begin{aligned} F_1(k) &= \dfrac{1}{2} \left[ X(k) + X(k + \frac{N}{2})\right] \\ F_2(k) &= \dfrac{1}{2} \left[ X(k) - X(k + \frac{N}{2})\right] W_N^{-k} \end{aligned} \]
其中 \(k = 0, 1, \cdots, \frac{N}{2} - 1\)。
令 \(y(n) = f_1(n) + j f_2(n)\)，易得 \(Y(k) = F_1(k) + j F_2(k)\)。
已知 \(y(n)\)，则可以通过计算 \(Y(k)\)（\(\frac{N}{2}\) 点 DFT）得到 \(X(k)\)（\(N\) 点 DFT）；
已知 \(Y(k)\)，则可以通过计算 \(y(n)\)（\(\frac{N}{2}\) 点 IDFT）得到 \(x(n)\)（\(N\) 点 IDFT）。

蝶形运算

对一对复数 (a, b)，将 b 乘上 \(W_N^k\) 然后 a 加上或减去这个积，得到一对新的复数 (A, B)。

对于一个 2 点序列 \(x(0)\)、\(x(1)\)，它们各自做 1 点 DFT，得到的结果就是它们自身，因此按上述基 2 FFT 算法，它们的 2 点 DFT 结果为 \(X(0) = x(0) + W_2^0 x(1)\)、\(X(1) = x(0) - W_2^0 x(1)\)。
即上述蝶形运算中 \(a = x(0)\)、\(b = x(1)\)、\(A = X(0)\)、\(B = X(1)\)。

eg. 计算 8 点基 2 FFT

• 序列抽取
  • 第一次抽取：
    • \(x(0), x(2), x(4), x(6)\)
    • \(x(1), x(3), x(5), x(7)\)
      
  • 第二次抽取：
    • \(x(0), x(4)\)
    • \(x(2), x(6)\)
    • \(x(1), x(5)\)
    • \(x(3), x(7)\)
      
  • 第三次抽取：
    \(y(n)\) 表示当前位置存储的值。
    • \(y(0) = x(0)\)
      
    • \(y(1) = x(4)\)
    • \(y(2) = x(2)\)
    • \(y(3) = x(6)\)
    • \(y(4) = x(1)\)
    • \(y(5) = x(5)\)
    • \(y(6) = x(3)\)
    • \(y(6) = x(7)\)
• DFT 计算
  • 第一阶段：4 级 4 组蝶形运算
    • \(x(0), x(4)\) 的 2 点 DFT 结果：
      \(y(0) = y(0) + W_2^0 y(1) = y(0) + W_8^0 y(1)\)
      \(y(1) = y(0) - W_2^0 y(1) = y(0) - W_8^0 y(1)\)
    • 后续 \(y(n)\) 的算法相同
  • 第二阶段：2 级 4 组蝶形运算
    • \(x(0), x(2), x(4), x(6)\) 4 点 DFT 结果：
      \(y(0) = y(0) + W_4^0 y(2) = y(0) + W_8^0y(2)\)
      \(y(1) = y(1) + W_4^1 y(3) = y(1) + W_8^2y(3)\)
      \(y(2) = y(0) - W_4^0 y(2) = y(0) - W_8^0y(2)\)
      \(y(3) = y(1) - W_4^1 y(3) = y(1) - W_8^2y(3)\)
    • \(x(1), x(3), x(5), x(7)\) 4 点 DFT 运算方式相同
  • 第三阶段：1 级 4 组蝶形运算
    • \(y(0) = y(0) + W_8^0 y(4)\)
      \(y(1) = y(1) + W_8^1 y(5)\)
      \(y(2) = y(2) + W_8^2 y(6)\)
      \(y(3) = y(3) + W_8^3 y(7)\)
      \(y(4) = y(0) - W_8^0 y(4)\)
      \(y(5) = y(1) - W_8^1 y(5)\)
      \(y(6) = y(2) - W_8^2 y(6)\)
      \(y(7) = y(3) - W_8^3 y(7)\)
    • \(X(k) = y(k)\)

按频率抽取

将 DFT 公式分为前 \(\frac{N}{2}\) 个数据点的求和与后 \(\frac{N}{2}\) 个数据点的求和： \[ \begin{aligned} X(k) &= \sum\limits_{m=0}^{(N/2)-1} x(n) W_N^{kn} + \sum\limits_{m=N/2}^{N-1} x(n) W_N^{kn} \\ & = \sum\limits_{m=0}^{(N/2)-1} x(n) W_N^{kn} + W_N^{kn} \sum\limits_{m=0}^{(N/2)-1} x(n + \frac{N}{2}) W_N^{kn} \\ & = \sum\limits_{m=0}^{(N/2)-1} \left[ x(n) + (-1)^k x(n + \frac{N}{2})\right] W_N^{kn} \end{aligned} \]

将 \(X(k)\) 分拆成奇序列和偶序列：
\[ \begin{aligned} X(2k) &= \sum\limits_{m=0}^{(N/2)-1} \left[ x(n) + x(n + \frac{N}{2})\right] W_{N/2}^{kn} \\ X(2k+1) &= \sum\limits_{m=0}^{(N/2)-1} \left\{ \left[ x(n) - x(n + \frac{N}{2})\right] W_N^{n} \right\} W_{N/2}^{kn} \end{aligned} \]

定义 \(\frac{N}{2}\) 点序列：
\[ \begin{aligned} g_1(n) &= x(n) + x(n + \frac{N}{2}) \\ g_2(n) &= \left[ x(n) - x(n + \frac{N}{2})\right] W_N^{n},\quad n = 0,1, \cdots, \frac{N}{2}-1 \end{aligned} \]

则有 \[ \begin{aligned} X(2k) &= \sum\limits_{m=0}^{(N/2)-1} g_1(n) W_{N/2}^{kn} \\ X(2k+1) &= \sum\limits_{m=0}^{(N/2)-1} g_2(n) W_{N/2}^{kn} \end{aligned} \]

与按时间抽取类似的，该过程可以重复抽取。
整个计算过程的抽取次数、复数乘法次数、复数加法次数均与按时间抽取相同。

基4 FFT 算法

\(N = 4^v\)

与基2 FFT 算法的推导方式类似，将隔 2 个进行一次抽取改为隔 4 个进行一次抽取。
令 \[ \begin{aligned} x_0(n) &= x(4n) \\ x_1(n) &= x(4n+1) \\ x_2(n) &= x(4n+2) \\ x_3(n) &= x(4n+3),\quad 0 \leq n \leq \frac{N}{4} - 1 \end{aligned} \]

可以得到 \[ \left[ \begin{matrix} X(k) \\ X(k + \frac{N}{4}) \\ X(k + \frac{N}{2}) \\ X(k + \frac{3N}{4}) \end{matrix} \right] = \left[ \begin{matrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \end{matrix} \right] \cdot \left[ \begin{matrix} X_0(k) \\ W_N^k X_1(k) \\ W_N^2k X_2(k) \\ W_N^3k X_3(k) \end{matrix} \right] \]

按此方法继续抽取即可。

N 为复合数的 FFT 算法

将 N 分解为两个整数的乘积：\(N = LM\)。如果 N 是质数可以通过加 0 的方式使其满足该分解。
令 \(0 \leq l \leq L-1\)，\(0 \leq m \leq M-1\)；则可以选择映射 \(n = Ml + m\)，即 \(x(n) = x(l, m)\)，相当于将原序列分为 L 组，每组的长度为 M。
同时 \(X(k)\) 被映射为 \(X(p, q)\)，即 \(k = Mp + q\ (0 \leq p \leq L-1,\ 0 \leq q \leq M-1)\)。

• 对每一组计算 M 点 DFT： \[ F(l, q) = \sum\limits_{m=0}^{M-1} x(l, m) W_M^{mq} \]

• 计算 L 点 DFT：
  \[ X(p, q) = \sum\limits_{l=0}^{L-1} W_N^{lq} F(l, q) W_L^{lp} \]

需要复数乘法 \(N(M+L+1)\) 次，复数加法 \(N(M+L-2)\) 次。

同理可以选择映射 \(n = Lm + l\)，即将原序列分为 M 组，每组长度为 L。
计算方法类似。

调频 z 变换算法

N 点数据序列 x(n) 的 DFT 可视为单位圆上均匀分布的 N 点 z 变换。

计算 \(x(n)\) 在点集 \(\{z_k\}\) 上的 z 变换
其中等高线是一个以 r 为半径的圆，且 \(z_k\) 是 N 个等间距点，即 \(z_k = r W_N^{k}\)
则有：
\[ \begin{aligned} X(z_k) &= \sum\limits_{n=0}^{N-1} x(n) z_k^{-n} \\ &= \sum\limits_{n=0}^{N-1} [x(n) r^{-n}] W_N^{-kn} \end{aligned} \]
即计算 \(x(n) r^{-n}\) 的 DFT。

离散时间系统

系统函数

用 N 阶差分方程描述：
\[ y(n) = -\sum\limits_{k=1}^{N} a_k y(n-k) + \sum\limits_{k=0}^{M} b_k x(n-k) \]

对应的系统函数：
\[ H(z) = \dfrac{Y(z)}{X(z)}= \dfrac{\sum\limits_{k=0}^{M} b_k z^{-k}}{1 + \sum\limits_{k=1}^{N} a_k z^{-k}} \]

• \(a_k \neq 0,\ b_k \neq 0\)（既有零点又有极点）：零极点系统，IIR 系统
• \(a_k = 0\)（只有零点）：全零点系统，FIR 系统
• \(b_k = 0\)（包含 N 个极点和一个 z=0 的 N 阶零点）：全极点系统，IIR 系统

FIR 系统

具有有限长冲激响应的系统。
系统函数里不包含极点。

差分方程描述：
\[ y(n) = \sum\limits_{k=0}^{M-1} b_k x(n-k) \]

系统函数描述：
\[ H(z) = \sum\limits_{k=0}^{M-1} b_k z^{-k} \]

单位冲激响应：
\[ h(n) = \left\{ \begin{aligned} &b_n,\quad 0 \leq n \leq M-1 \\ &0,\quad else \end{aligned} \right. \]

直接型结构

(P. 419)
对 \(b_k x(n-k)\)，将 \(x(n)\) 经过 k 个延时器（即 \(z^{-1}\)）后乘 \(b_k\)（即 \(h(k)\)），在将每个 k 值对应的结果通过加法器相加得到结果。

该结构需要 M-1 个存储空间来存放 M-1 个输入，每个输出需要 M 次乘法和 M-1 次加法。

当 FIR 系统具有线性相位时，单位冲激响应函数服从下面对称或不对称条件：
\[ h(n) = \pm h(M-1-n) \]

对称的部分可以在延时后先相加，再一起乘对应的 \(h(k)\)，因此乘法减为 \(\frac{M}{2}\)（M 为偶数）次或 \(\frac{M - 1}{2}\)（M 为奇数）次。

级联型结构

将 \(H(z)\) 拆分成多个 \(H_k(z)\) 相乘。
用直接型结构实现 \(H_k(z)\) 再将它们级联起来。

IIR 系统

具有无限长冲激响应的系统。
系统函数里包含极点。

直接型结构

差分方程描述和系统函数描述与系统函数中的相同。

可以视为级联的两个系统，即 \(H(z) = H_1(z) H_2(z)\)
其中 \(H_1(z)\) 包含 \(H(z)\) 的零点，\(H_2(z)\) 包含 \(H(z)\) 的极点，即 \[ H_1(z) = \sum\limits_{k=0}^{M} b_k z^{-k} \] \[ H_2(z) = \dfrac{1}{1 + \sum\limits_{k=1}^{N} a_k z^{-k}} \]

这里面 \(H_1(z)\) 是一个 FIR 系统，\(H_2(z)\) 是全极点系统。

直接 I 型

将 \(H_1(z)\) 和 \(H_2(z)\) 分别使用直接型结构实现再级联起来即可。

需要 M+N+1 次乘法，M+N 次加法和 M+N+1 个存储空间。

直接 II 型

交换直接 I 型中 \(H_1(z)\) 和 \(H_2(z)\) 对应结构的位置，即将 \(H_2(z)\) 置于 \(H_1(z)\) 之前，二者公用一组延迟器。

需要 M+N+1 次乘法，M+N 次加法和 max{M, N} 个存储空间。

级联型结构

与 FIR 的级联型结构类似，同样是将 \(H(z)\) 拆分成多个 \(H_k(z)\) 相乘，每个 \(H_k(z)\) 用直接型结构实现再级联起来。

并联型结构

假设 \(N \geq M\) 且极点是独立的，对 \(H(z)\) 使用部分分式展开法： \[ H(z) = C + \sum\limits_{k=1}^{N} \dfrac{A_k}{1 - p_k z^{-1}} = C + \sum\limits_{k=1}^{K} H_k(z) \] 其中，\(p_k\) 是极点，\(A_k\) 是部分分式展开的系数，\(C = \frac{b_N}{a_N}\)。

分别用直接型结构实现 \(H_k(z)\) 再并联起来即可。

滤波器设计

数字滤波器

技术指标

技术指标以其频率响应的幅度特性给出，有通带、阻带和过渡带。

• 通带
  • 通带截止频率：\(\omega_p\)
  • 通带波纹：\(\delta_1\)
  • 通带允许最大衰减：\(R_p \geq -20\lg (1-\delta_1)\ (dB)\)
• 阻带
  • 阻带截止频率：\(w_{st}\)
  • 阻带波纹：\(\delta_2\)
  • 阻带应达到的最小衰减：\(A_s \leq -20\lg (\delta_2)\ (dB)\)

频率响应

• 幅度平方响应：
  \[ | H(e^{j\omega}) |^2 = H(z)H(z^{-1}) |_{z=e^{j\omega}} \]

• 相位响应：
  \[ \beta (e^{j\omega}) = \dfrac{1}{2j} \ln \left[ \dfrac{H(z)}{H(z^{-1})}\right]_{z = e^{j\omega}} \]

• 群延迟响应：
  \[ \tau (e^{j\omega}) = -\dfrac{d \beta (e^{j\omega})}{d \omega} \]
  是滤波器平均延时的一个度量。
  当滤波器满足线性相位响应特性时，通带内群延迟为常数。

IIR 滤波器

模拟滤波器

巴特沃斯低通滤波器

已知的参数 - 通带截止频率：\(\Omega_p\) - 阻带截止频率：\(\Omega_{st}\) - 通带最大衰减 (dB)：\(R_p\) - 阻带最小衰减 (dB)：\(A_s\)

要求的参数 - 滤波器阶数：N
\[ N = \dfrac{\lg \left( \dfrac{10^{0.1 A_s} - 1}{10^{0.1 R_p} - 1}\right)}{2\lg \left( \dfrac{\Omega_{st}}{\Omega_p}\right)} \]

特殊情况：\(R_p = 3dB,\ \Omega_p = \Omega_c\)
\[ N \geq \dfrac{\lg (10^{0.1 A_s} - 1)}{2\lg \left( \dfrac{\Omega_{st}}{\Omega_p}\right)} \]

• 截止频率：\(\Omega_c\)
  \[ \Omega_{cp} = \dfrac{\Omega_p}{\sqrt[2N]{10^{0.1R_p} - 1}} \leq \Omega_c \leq \dfrac{\Omega_{st}}{\sqrt[2N]{10^{0.1A_s} - 1}} = \Omega_{cs} \]

得到滤波器的系统函数：
\[ H_a(s) = \dfrac{1}{\sqrt{1 + \left( \dfrac{s}{j\Omega_c}\right)^{2N}}} \]

记归一化原型低通滤波器的系统函数为 \(H_{an}(s)\)，则有：
\[ H_a(s) = H_{an} = \left( \dfrac{s}{\Omega_c}\right) \]

实际应用时可以根据所求参数查表得到 \(H_{an}(s)\)，再去归一化得到 \(H_a(s)\)

切比雪夫 I 型低通滤波器

已知的参数 - 通带截止频率：\(\Omega_p\) - 阻带截止频率：\(\Omega_{st}\) - 通带最大衰减 (dB)：\(R_p\) - 阻带最小衰减 (dB)：\(A_s\)

要求的参数 - 截止频率，指定 \(\Omega_c = \Omega_p\)

• 波纹参数：\(\varepsilon = \sqrt{10^{0.1R_p} - 1}\)

• 滤波器阶数：N
  \[ N \geq \dfrac{\lg \left[ \sqrt{\dfrac{10^{0.1 A_s} - 1}{10^{0.1 R_p} - 1}} + \sqrt{\dfrac{10^{0.1 A_s} - 1}{10^{0.1 R_p} - 1} - 1}\right]}{\lg \left[ \dfrac{\Omega_{st}}{\Omega_p} + \sqrt{\left( \dfrac{\Omega_{st}}{\Omega_p}\right)^2 - 1}\right]} \]

由 N 和 \(\varepsilon\) 查表可得 \(H_{an}(s)\)；由 \(\Omega_c\) 去归一化可得 \(H_a(s)\)

模拟频域频带变换

将归一化模拟低通滤波器转换成模拟低通、高通、带通、带阻滤波器。

记归一化 \((\bar{\Omega}_p = 1)\) 模拟低通滤波器的系统函数为 \(H_{an}(\bar{s})\)，频率响应为 \(H_{an}(j\bar{\Omega})\)。
变换后各各类模拟滤波器的系统函数为 \(H_a(s)\)，频率响应为 \(H_a(j\Omega)\)。

\(H_a(s)\)	\(\bar{s}\)	\(\bar{\Omega}\)	补充
低通 \(H_{lp}(s)\) 通带截止频率 \(\Omega_{p}\)	\(\dfrac{s}{\Omega_{p}}\)	\(\dfrac{\Omega}{\Omega_{p}}\)
高通 \(H_{hp}(s)\) 通带截止频率 \(\Omega_{p}\)	\(\dfrac{\Omega_{p}}{s}\)	\(-\dfrac{\Omega_{p}}{\Omega}\)
带通 \(H_{bp}(s)\) 通带截止频率 \(\Omega_{p1}\) \(\Omega_{p2}\)	\(\dfrac{s^2 + \Omega_{p0}^2}{B_p s}\)	\(\dfrac{\Omega^2 - \Omega_{p0}^2}{B_p \Omega}\)	\(\Omega_{p0} = \sqrt{\Omega_{p1}\Omega_{p2}}\) \(B_p = \Omega_{p2} - \Omega_{p1}\)
带阻 \(H_{bs}(s)\) 阻带截止频率 \(\Omega_{st1}\) \(\Omega_{st2}\)	\(\dfrac{\bar{\Omega}_{st} B_s s}{s^2 + \Omega_{st0}^2}\)	\(\dfrac{\bar{\Omega}_{st} B_s \Omega}{s^2 + \Omega_{st0}^2}\)	\(\Omega_{st0} = \sqrt{\Omega_{st1}\Omega_{st2}}\) \(B_s = \Omega_{st2} - \Omega_{st1}\)

• 带通滤波器计算步骤
  • 频带变换
    • 按上表中的公式求 \(\bar{\Omega}_{p0},\ B_p,\ \bar{\Omega}_{st1},\ \bar{\Omega}_{st2}\)
    • 归一化模拟低通滤波器的阻带截止频率为 \(\bar{\Omega}_{st} = \min(|\bar{\Omega}_{st1}|, |\bar{\Omega}_{st2}|)\)
  • 利用 \(\bar{\Omega}_{p} = 1,\ \bar{\Omega}_{st},\ R_p,\ A_s\) 按照相应低通滤波器的求法求出 \(H_{an}(\bar{s})\)
    • 如果是巴特沃斯滤波器，当 \(R_p \neq 3dB\) 时，\(\Omega_p \neq \Omega_c\)，需要对求出的系统函数做去归一化，得到 \(H_a(\bar{s}) = H_{an}\left(\dfrac{\bar{s}}{\bar{\Omega}_c}\right)\)
  • 再次利用上表中的公式，将 \(\bar{s}\) 变换为 \(s\) 从而得到 \(H_{bp}(s)\)
• 阻带滤波器计算步骤
  • 频带变换
    • 按上表中的公式求 \(\bar{\Omega}_{st0},\ B_s,\ \bar{\Omega}_{p1},\ \bar{\Omega}_{p2}\)
    • 归一化模拟低通滤波器的通带截止频率为 \(\bar{\Omega}_{p} = \max(|\bar{\Omega}_{p1}|, |\bar{\Omega}_{p2}|)\)（得到的是有关 \(\bar{\Omega}_{st}\) 的式子）
    • 依据 \(\bar{\Omega}_p = 1\) 计算出 \(\bar{\Omega}_{st}\)
  • 利用 \(\bar{\Omega}_{p} = 1,\ \bar{\Omega}_{st},\ R_p,\ A_s\) 按照相应低通滤波器的求法求出 \(H_{an}(\bar{s})\)
    • 如果是巴特沃斯滤波器，当 \(R_p \neq 3dB\) 时，\(\Omega_p \neq \Omega_c\)，需要对求出的系统函数做去归一化，得到 \(H_a(\bar{s}) = H_{an}\left(\dfrac{\bar{s}}{\bar{\Omega}_c}\right)\)
  • 再次利用上表中的公式，将 \(\bar{s}\) 变换为 \(s\) 从而得到 \(H_{bp}(s)\)

间接法设计数字滤波器

冲激响应不变法

将模拟滤波器的单位冲激响应 \(h_a(t)\) 映射成数字滤波器的单位冲激响应 \(h(n)\)。
记抽样周期为 \(T = \frac{1}{f_s}\)，满足 \(h(n) = h_a(nT)\)。

要求： - 要求模拟滤波器是严格带限于 \(\frac{f_s}{2}\)，不能用于设计高通滤波器和带阻滤波器。 - 只适用于并联结构的系统函数，即系统函数必须先展开成部分分式。

计算步骤： - 计算出相应模拟滤波器的系统函数 \(H_a(s)\)

• 将该系统函数展开为部分分式之和的形式，即
  \[ H_a(s) = \sum\limits_{k=1}^{N} \dfrac{A_k}{s - s_k} \]

• 依据 \(H_a(s)\) 和抽样周期 T 得到数字滤波器的系统函数
  \[ H(z) = \sum\limits_{k=1}^{N} \dfrac{A_k}{1 - e^{s_k T} z^{-1}} \]

双线性变换法

将模拟滤波器的系统函数 \(H_a(s)\) 映射成数字滤波器系统函数 \(H(z)\)。

计算步骤： - 根据模拟角频率 \(\Omega_p\)、\(\Omega_{st}\) 计算数字频率
\[ \omega_p = \Omega_p T,\ \omega_{st} = \Omega_{st} T \]

• 计算频率预畸
  \[ \Omega_p' = \dfrac{2}{T} \tan \left( \dfrac{\omega_p}{2}\right) \]
  \[ \Omega_{st}' = \dfrac{2}{T} \tan \left( \dfrac{\omega_{st}}{2}\right) \]

• 根据 \(\Omega_p',\ \Omega_{st}',\ R_p,\ A_s\) 计算出模拟滤波器的 \(H_a(s)\)

• 令 \(H_a(s)\) 中
  \[ s = \dfrac{2}{T} \cdot \dfrac{1 - z^{-1}}{1 + z^{-1}} \]
  从而根据 \(H_a(s)\) 得到 \(H(z)\)

数字频域频带变换

要求将归一化模拟低通滤波器直接转换为低通、高通、带通、带阻数字滤波器。

eg. 带通数字滤波器 - 将模拟域指标变换为数字域指标 \(\omega_{st1},\ \omega_{st2},\ \omega_{p1},\ \omega_{p2}\)

• 求中心频率
  \[ \cos(\omega_0) = \dfrac{\cos\left( \dfrac{\omega_{p2} + \omega_{p1}}{2}\right)}{\cos\left( \dfrac{\omega_{p2} - \omega_{p1}}{2}\right)} \]

• 将数字域频率指标转换，得到频域的通带截止频率和阻带截止频率

  • 选择 \(\omega_p = \max(\omega_{p1}, \omega_{p2})\)
  • 变换到频域
    \[ \Omega = \dfrac{\cos \omega_0 - \cos \omega}{\sin \omega} \]
    
  • 取 \(\Omega_{st} = \min(|\Omega_{st1}|, |\Omega_{st2}|)\)

• 根据 \(\Omega_p,\ \Omega_{st},\ R_p,\ A_s\) 按相应的模拟滤波器的设计方法得到 \(H_{a}(s)\)

  • 令 \(H_a(s)\) 中
    \[ s = \dfrac{1 - 2z^{-1} \cos \omega_0 + z^{-2}}{1 - z^{-2}} \]
    得到 \(H_{bp}(z)\)

懒得改了。。。。

FIR 滤波器

滤波器特点

• h(n) 具有对称性
  \[ h(n) = \pm h(N-1-n) \]

• 群延迟为常数
  \[ \tau = \dfrac{N-1}{2} \]

有四种情况： - \(h(n) = h(N-1-n)\)，N 为奇数：
可用于设计任一种滤波器。 - \(h(n) = h(N-1-n)\)，N 为偶数：
不能设计高通和带阻。 - \(h(n) = -h(N-1-n)\)，N 为奇数：
只能用于带通。 - \(h(n) = -h(N-1-n)\)，N 为偶数：
- 只能用于高通和带通。

窗函数设计法

给定理想频率响应 \(H_d(e^{j\omega})\)，导出 \(h_d(n)\)
对 \(h_d(n)\) 截取有限长度，形成 \(h(n)\)

计算步骤

• 求参数

  • 过渡带宽：\(\Delta \omega = \omega_{st} - \omega_p\)

  • 截止频率：\(\omega_c = \dfrac{\omega_p + \omega_{st}}{2}\)

  • 阻带最小衰减：
    \[ A_s = -20\lg \left( \dfrac{\delta}{1 + \delta}\right) \]
    因此给定 \(A_s\) 可以得到波纹 \(\delta\)：
    \[ \delta = \dfrac{10^{-\frac{A_s}{20}}}{1 - 10^{-\frac{A_S}{20}}} \]

  • 通带最大衰减：
    \[ R_p = -20\lg \left( \dfrac{1 - \delta}{1 + \delta}\right) \]

• 求 \(h_d(n)\)

  • 低通滤波器
    \[ h_d(n) = \left\{ \begin{aligned} &\dfrac{\sin [\omega_c (n-\tau)]}{\pi (n-\tau)},\ n \neq \tau \\ &\dfrac{\omega_c}{\pi},\ n = \tau \end{aligned} \right. \]

  • 带通滤波器
    \[ h_d(n) = \left\{ \begin{aligned} &\dfrac{\sin [\omega_2 (n-\tau)] - \sin [\omega_1 (n-\tau)]}{\pi (n-\tau)},\ n \neq \tau \\ &\dfrac{\omega_2 - \omega_1}{\pi},\ n = \tau \end{aligned} \right. \]

  • 带阻滤波器
    \[ h_d(n) = \left\{ \begin{aligned} &\dfrac{\sin [\pi (n-\tau)] - \sin [\omega_2 (n-\tau)] + \sin [\omega_1 (n-\tau)]}{\pi (n-\tau)},\ n \neq \tau \\ &1 - \dfrac{\omega_2}{\pi} + \dfrac{\omega_1}{\pi},\ n = \tau \end{aligned} \right. \]

  • 高通滤波器
    \[ h_d(n) = \left\{ \begin{aligned} &\dfrac{\sin [\omega_c (n-\tau)]}{\pi (n-\tau)},\ n \neq \tau \\ &1 - \dfrac{\omega_c}{\pi},\ n = \tau \end{aligned} \right. \]

• 加窗截断，得到 \(h(n) = h_d(n) w(n)\)
  记
  \[ R_N(n) = \left\{ \begin{aligned} &1,\ -\frac{N}{2} \leq n \leq \frac{N}{2} \\ & 0,\ else \end{aligned} \right. \]

  • 矩形窗
    \[ w(n) = R_N(n) \]

  • 三角形窗
    \[ w(n) = \left\{ \begin{aligned} &\dfrac{2n}{N-1},\quad 0 \leq n \leq \frac{N-1}{2} \\ &2 - \dfrac{2n}{N-1},\quad \frac{N-1}{2} < n \leq N-1 \end{aligned} \right. \]

  • 汉宁 (Hanning) 窗
    \[ w(n) = \dfrac{1}{2}\left[1 - \cos\left( \dfrac{2\pi n}{N-1}\right)\right] R_N(n) \]

  • 海明 (Hamming) 窗
    \[ w(n) = \left[0.54 - 0.46 \cos\left( \dfrac{2\pi n}{N-1}\right)\right] R_N(n) \]

  • 布莱克曼 (Blackman) 窗
    \[ w(n) = \left[0.42 - 0.5 \cos\left( \dfrac{2\pi n}{N-1}\right) + 0.08 \cos\left( \dfrac{4\pi n}{N-1}\right)\right] R_N(n) \]

频率抽样设计法

基本思路

• 给定期望逼近的频率响应 \(H_d(e^{j\omega})\)。

• 在 \(\omega\) 的一个周期 \([0, 2\pi)\) 中进行 N 点等间距抽样得到 \(H(k)\)，即
  \[ H(k) = H_d(e^{j\omega}) |_{\omega = \frac{2\pi}{N}k},\quad k = 0, 1, \cdots, N-1 \]

• 求 \(h(n)\)
  \[ h(n) = IDFT[H(k)] = \dfrac{1}{N} \sum\limits_{k=0}^{N-1} H(k) e^{j\frac{2\pi}{N} nk},\quad n = 0, 1, \cdots, N-1 \]

设计公式

PPT 上懒得打了。

计算方法：
- 在通带部分，\(|H(k)| = 1\)；在阻带部分，\(|H(k)| = 0\)；过渡带的值一般会给出。

• 通过已知的截止频率在整个周期中的占比，确定通带、阻带、过渡带各有哪几个点。

• 结合上述设计公式得到 \(H(k)\)

与 IIR 比较

• 可以实现严格的线性相位。
• h(n) 有限长，一定稳定。
• 相同的技术指标，FIR 所需阶数比 IIR 高得多。