Appearance
数字信号处理（DSP）
概念
前面四章把工具备齐了：卷积描述系统、傅里叶看频率成分、拉普拉斯给连续系统定性、Z 变换把这套语言搬到采样序列上。数字信号处理（Digital Signal Processing，DSP）就是把它们变成计算机能跑的算法。
一句话定义：DSP 是"对一串采样值做算术，来改变或看清它包含的频率成分"。
它只有三个动作：
| 动作 | 输入 | 输出 | 用什么 |
|---|---|---|---|
| 变换 | 时域序列 x[n] | 频域 X[k] | DFT / FFT |
| 滤波 | 时域序列 x[n] | 时域序列 y[n] | 差分方程（IIR）或卷积（FIR） |
| 分析 | 一段长序列 | 谱估计 | 加窗 + FFT + 平均 |
这三个动作里的"变换"是整门课的心脏：只要能把一段采样变成频谱，滤波、检测、压缩、识别就都有了统一的出发点。而变换之所以在工程上可用，靠的是一次算法层面的革命——FFT。
DSP 还回答了一个上一章留下的问题：H(z) 写出来到底是什么？答案是一个循环里的一行（上一章结尾已经见过 y = 0.5*y + x）。DSP 做的就是系统地把这只"一行"组织成滤波器、频谱仪和通信机的算法。
原理
一、DFT：把 N 个采样变成 N 条谱线
离散傅里叶变换（DFT）的定义就是上一章 Z 变换在单位圆上均匀取 N 点：
X[k] = Σ_{n=0}^{N−1} x[n] · e^(−j2πnk/N) k = 0, 1, …, N−1它的物理含义可以直接读出来：
谱线 k | 对应的物理频率 | 说明 |
|---|---|---|
k = 0 | 0 Hz（直流） | 序列的平均值乘 N |
k = 1 | fs/N | 最低的一个"非零频" |
k = N/2 | fs/2（奈奎斯特） | 最高可表示的频率 |
k > N/2 | 与 N−k 共轭对称 | 实序列的下半段只是镜像 |
三条必须记住的缩放关系：
- 频率分辨率
Δf = fs / N（等价地Δf = 1/T，T 是采集时长）； - 单频幅度：解析幅度为
A的余弦，其|X[k]| = A·N/2——所以要从谱线读幅度，得乘2/N； - 实序列只算一半：
X[k]与X[N−k]共轭，有效信息只有N/2 + 1条。
二、FFT：同一个 DFT，快 N/log₂N 倍
按定义直算 DFT 要 N² 次复数乘法。FFT 不是近似，它靠"分治"把同样的结果算得更快：把 N 点序列按奇偶拆成两个 N/2 点 DFT，再合并——这个合并动作叫蝶形：
A' = A + W·B
B' = A − W·B W = e^(−j2πk/N)一共 log₂N 级，每级 N/2 个蝶形：
| N | DFT 复数乘 N² | FFT 复数乘 (N/2)log₂N | 加速比 |
|---|---|---|---|
| 8 | 64 | 12 | 5.33 |
| 64 | 4096 | 192 | 21.33 |
| 1024 | 1048576 | 5120 | 204.80 |
| 8192 | 67108864 | 53248 | 1260.31 |
结果一模一样，代价差三个数量级——这就是为什么 1965 年 FFT 一出来，频谱分析从"做不起"变成了"随手做"。
读这张表的两个要点：
- 加速比随
N增长而增长（N/log₂N），所以 FFT 对大N才划算； - radix-2 要求
N是 2 的幂，先补零到最近的 2 的幂即可（补零的代价见下面第四节）。
三、加窗：泄漏与主瓣宽度是一笔交易
DFT 只看到有限的一段，等于把无限长的信号乘了一个矩形窗。时域相乘 = 频域卷积，于是原来一条干净的谱线被"抹开"成 sinc 形状——这条 sinc 的裙摆就是频谱泄漏，旁瓣会把旁边的小信号淹掉。
换一个窗函数，就是换一段"更平滑的截断"。四类常用窗的实测指标（本页示例 2 用窗的离散时间傅里叶变换算出）：
| 窗函数 | 主瓣总宽（谱线数） | 最大旁瓣 | 相干增益 | 最大扇形损失 |
|---|---|---|---|---|
| 矩形 | 2.00 | −13.25 dB | 1.0000 | −3.92 dB |
| 汉宁 | 4.00 | −31.47 dB | 0.5000 | −1.42 dB |
| 汉明 | 4.00 | −42.45 dB | 0.5400 | −1.75 dB |
| 布莱克曼 | 6.00 | −58.11 dB | 0.4200 | −1.10 dB |
这张表要横着读：旁瓣压得越低，主瓣就越宽。矩形窗主瓣最窄但旁瓣只有 −13 dB；布莱克曼把旁瓣压到 −58 dB，代价是主瓣宽到 3 倍。
两个必须会用窗的修正量：
- 相干增益：
Σw[n]/N。加窗后幅度整体变小，读出来的峰值要除以相干增益才是真实幅度（汉宁窗习惯上乘 2）。 - 扇形损失（scalloping loss）：真实频率落在两条谱线正中间时，读到的峰值比真值低。矩形窗最多低 3.92 dB，加窗后缩到 1~1.8 dB——这是加窗的第二份好处。
选窗的口诀：只要幅度精度、不管泄漏 → 矩形；要把强信号旁边的小信号挖出来 → 汉明或布莱克曼。这正是 数字音频 与 数字图像处理 里做谱分析时的第一道选择题。
四、频率分辨率：只有采集时长说了算
一个极常见的误会是"补零能提高分辨率"。不是。
Δf = fs / N = 1 / (N·T_s) = 1 / T T 是采集时长决定分辨率的是 T，不是 N，更不是补了多少零。
| 采样率 | 点数 | 采集时长 T | 分辨率 Δf |
|---|---|---|---|
| 8 kHz | 256 | 0.0320 s | 31.2500 Hz |
| 8 kHz | 1024 | 0.1280 s | 7.8125 Hz |
| 44.1 kHz | 1024 | 0.0232 s | 43.0664 Hz |
| 44.1 kHz | 4096 | 0.0929 s | 10.7666 Hz |
| 48 kHz | 48000 | 1.0000 s | 1.0000 Hz |
注意后两行：同样 1024 点，采样率越高分辨率越差——因为 1024 点在 44.1 kHz 下只录了 23 ms。想分清两个靠得很近的频率，唯一的办法是录得更久。
补零的作用是把谱线画得更密（曲线更光滑），让人眼看得清峰值落在哪，但它不会创造出原来不存在的信息（示例 2 第 5 节会用一个"凹陷比"把这件事量化出来）。
五、FIR 与 IIR：两种滤波器的分工
| FIR（有限冲激响应） | IIR（无限冲激响应） | |
|---|---|---|
| 结构 | y[n] = Σ b_k x[n−k] | y[n] = Σ b_k x[n−k] + Σ a_k y[n−k] |
| 有没有反馈 | 没有 | 有（输出回到输入） |
| 相位 | 可做到严格线性 | 一般非线性 |
| 达到同样陡度 | 需要更多阶 | 少数几阶就够 |
| 稳定性 | 恒稳定 | 要看极点是否在单位圆内 |
| 实现代价 | 每点 N 次乘加 | 每点几次乘加 |
线性相位的条件与代价：h[n] = ±h[N−1−n]，此时群延迟恒为 (N−1)/2 拍。它意味着所有频率被延迟相同的时间——波形的形状不会失真，这对通信的鉴相、音频的相位一致性都是硬要求。
选择的分界：
- 要求相位线性、或要求绝对稳定（比如航空电子）→ FIR；
- 只要陡峭的幅度特性、算力紧张（比如语音编码）→ IIR。
六、定点实现：Q 格式与溢出
真实 DSP 芯片大多没有浮点单元，数用定点表示。最常见的 Q15：
Q15：1 位符号 + 15 位小数，可表示范围 [−1, +1)
0.5 → 16384 = 0x4000
0.75 → 24576 = 0x6000乘法的位宽会变宽：16 位 × 16 位 = 32 位，所以累加器必须是 32 位甚至 40 位，只在最后一步截回 16 位。这一步截断有两种策略：
- 环绕（wrap）：直接丢掉高位 → 1.25 变成 −0.75，符号翻转，听感是爆音；
- 饱和（saturate）：超出就卡在最大值 → 1.25 变成 0.999969，只是削顶。
工程上一律选饱和。这也解释了为什么写滤波器时系数都要先归一化、为什么求和顺序会影响结果——定点运算里每一步都可能出界。
示例
例 1：手写 radix-2 FFT，与按定义直算逐项对照（C）
下面这段代码同时跑两条路：按定义直算 DFT（N² 次复数乘）和 radix-2 时域抽取 FFT（倒序重排 + 三级蝶形），把两边的结果逐项打印出来对比，最后按公式给出不同 N 的代价对照。
/* 同一个 DFT，两种算法：radix-2 FFT 与按定义直算，结果逐项对照、代价逐档对照 */
#include <stdio.h>
#include <math.h>
#define N 8
#define PI 3.14159265358979323846
static const double x[N] = {1.0, 2.0, 3.0, 4.0, 4.0, 3.0, 2.0, 1.0};
/* 路径一：按定义直算 X[k] = sum x[n] * e^(-j2nk/N) */
static void dft_direct(double *re, double *im)
{
for (int k = 0; k < N; k++) {
double sr = 0.0, si = 0.0;
for (int n = 0; n < N; n++) {
double th = -2.0 * PI * k * n / N;
sr += x[n] * cos(th);
si += x[n] * sin(th);
}
re[k] = sr;
im[k] = si;
}
}
/* 路径二：radix-2 时域抽取 FFT（倒序重排 + log2(N) 级蝶形） */
static void fft_radix2(double *re, double *im)
{
for (int i = 0; i < N; i++) {
re[i] = x[i];
im[i] = 0.0;
}
int j = 0; /* 位倒序重排 */
for (int i = 0; i < N; i++) {
if (i < j) {
double t = re[i];
re[i] = re[j];
re[j] = t;
}
int m = N >> 1;
while (m >= 1 && j >= m) {
j -= m;
m >>= 1;
}
j += m;
}
for (int len = 2; len <= N; len <<= 1) { /* 蝶形：每级 N/2 个 */
int half = len >> 1;
for (int s = 0; s < N; s += len) {
for (int t = 0; t < half; t++) {
double th = -2.0 * PI * t / len;
double wr = cos(th), wi = sin(th);
double br = re[s + half + t], bi = im[s + half + t];
double tr = wr * br - wi * bi;
double ti = wr * bi + wi * br;
re[s + half + t] = re[s + t] - tr;
im[s + half + t] = im[s + t] - ti;
re[s + t] += tr;
im[s + t] += ti;
}
}
}
}
int main(void)
{
double ar[N], ai[N], br[N], bi[N];
dft_direct(ar, ai);
fft_radix2(br, bi);
printf("%-6s %-12s %-12s %-12s\n", "k", "Re(FFT)", "Im(FFT)", "abs(FFT)");
for (int k = 0; k < N; k++) {
printf("%-6d %-12.6f %-12.6f %-12.6f\n", k, br[k], bi[k],
sqrt(br[k] * br[k] + bi[k] * bi[k]));
}
double worst = 0.0;
for (int k = 0; k < N; k++) {
double dr = fabs(br[k] - ar[k]);
double di = fabs(bi[k] - ai[k]);
double d = sqrt(dr * dr + di * di);
if (d > worst) {
worst = d;
}
}
printf("max difference between FFT and direct DFT = %.3e\n", worst);
printf("\n%-8s %-18s %-24s %-10s\n", "N", "DFT complex mul", "FFT complex mul", "speedup");
int sizes[4] = {8, 64, 1024, 8192};
for (int i = 0; i < 4; i++) {
int n = sizes[i];
double d = (double) n * (double) n;
int lg = 0;
while ((1 << (lg + 1)) <= n) {
lg++;
}
double f = (double) n / 2.0 * (double) lg;
printf("%-8d %-18.0f %-24.0f %-10.2f\n", n, d, f, d / f);
}
printf("\nSame spectrum, cost differs by three orders of magnitude.\n");
return 0;
}
c 本站为静态站，不提供在线运行；可复制到本地用 gcc / python 执行
预期输出：
k Re(FFT) Im(FFT) abs(FFT)
0 20.000000 0.000000 20.000000
1 -5.828427 -2.414214 6.308644
2 0.000000 0.000000 0.000000
3 -0.171573 -0.414214 0.448342
4 0.000000 0.000000 0.000000
5 -0.171573 0.414214 0.448342
6 0.000000 0.000000 0.000000
7 -5.828427 2.414214 6.308644
max difference between FFT and direct DFT = 7.700e-15
N DFT complex mul FFT complex mul speedup
8 64 12 5.33
64 4096 192 21.33
1024 1048576 5120 204.80
8192 67108864 53248 1260.31
Same spectrum, cost differs by three orders of magnitude.读数：两条路径的最大偏差只有 7.7e-15——这就是浮点舍入的极限，可以认为完全相同。所以 FFT 换来的 1260 倍加速，一个比特的信息都没丢。
再看这条谱本身：输入 1,2,3,4,4,3,2,1 是关于中点对称的实序列，所以 X[k] 是实数加一点共轭对称的虚部；X[0] = 20 正好是 8 个数的和；X[2] = X[4] = X[6] = 0 是全零谱线。谱线的对称性和零点都能用手推出来——这是校验自己 FFT 写得对不对最快的方法。
例 2：分辨率、加窗、补零、定点化一次验完（Python）
C 段证明了"两种算法同一个结果"。这一段把 DSP 里其余五件事补齐：FFT 与 DFT 的等价性、分辨率由谁决定、四类窗的实测指标、扇形损失与相干增益、补零为什么无效、Q15 的溢出长什么样。
import math
def pad(s, w):
dw = sum(2 if ord(c) > 0x2000 else 1 for c in str(s))
return str(s) + " " * max(0, w - dw)
def table(head, rows, gap=2):
data = [[str(c) for c in r] for r in rows]
w = [max([sum(2 if ord(c) > 0x2000 else 1 for c in str(head[i]))]
+ [sum(2 if ord(c) > 0x2000 else 1 for c in r[i]) for r in data]) + gap
for i in range(len(head))]
print(" " + "".join(pad(head[i], w[i]) for i in range(len(head))))
for r in data:
print(" " + "".join(pad(r[i], w[i]) for i in range(len(head))))
def fft(x):
n = len(x)
re = [float(v) for v in x]
im = [0.0] * n
j = 0
for i in range(n):
if i < j:
re[i], re[j] = re[j], re[i]
m = n >> 1
while m >= 1 and j >= m:
j -= m
m >>= 1
j += m
length = 2
while length <= n:
half = length >> 1
for s in range(0, n, length):
for t in range(half):
th = -2.0 * math.pi * t / length
wr, wi = math.cos(th), math.sin(th)
br, bi = re[s + half + t], im[s + half + t]
tr = wr * br - wi * bi
ti = wr * bi + wi * br
re[s + half + t] = re[s + t] - tr
im[s + half + t] = im[s + t] - ti
re[s + t] += tr
im[s + t] += ti
length <<= 1
return re, im
def dft(x):
n = len(x)
out = []
for k in range(n):
sr = si = 0.0
for i in range(n):
th = -2.0 * math.pi * k * i / n
sr += x[i] * math.cos(th)
si += x[i] * math.sin(th)
out.append(complex(sr, si))
return out
def win_make(kind, n):
w = []
for i in range(n):
if kind == "rect":
v = 1.0
elif kind == "hann":
v = 0.5 - 0.5 * math.cos(2 * math.pi * i / n)
elif kind == "hamming":
v = 0.54 - 0.46 * math.cos(2 * math.pi * i / n)
else:
v = 0.42 - 0.5 * math.cos(2 * math.pi * i / n) \
+ 0.08 * math.cos(4 * math.pi * i / n)
w.append(v)
return w
def win_spectrum(kind, n, pts):
w = win_make(kind, n)
out = []
for p in range(pts + 1):
om = math.pi * p / pts
sr = si = 0.0
for i in range(n):
sr += w[i] * math.cos(om * i)
si -= w[i] * math.sin(om * i)
out.append(abs(complex(sr, si)))
return out
print("=== 1. FFT 与直接 DFT：同一个答案，两种算法 ===")
x = [1, 2, 3, 4, 4, 3, 2, 1]
re, im = fft(x)
direct = dft(x)
rows = []
worst = 0.0
for k in range(len(x)):
d = abs(complex(re[k], im[k]) - direct[k])
worst = max(worst, d)
rows.append(["k = %d" % k, "%.6f" % re[k], "%.6f" % im[k], "%.6f" % abs(complex(re[k], im[k]))])
table(["谱线", "实部 Re", "虚部 Im", "幅度"], rows)
print(" 两种算法最大偏差 = %.3e（浮点舍入级别，可视为完全相同）" % worst)
rows = [["N", "DFT 复数乘 N^2", "FFT 复数乘 (N/2)log2N", "加速比"]]
for n in (8, 64, 1024, 8192):
lg = 0
while (1 << (lg + 1)) <= n:
lg += 1
d = float(n * n)
f = n / 2.0 * lg
rows.append(["%d" % n, "%d" % d, "%.0f" % f, "%.2f" % (d / f)])
table(rows[0], rows[1:])
print(" 结果一模一样，代价差三个数量级 —— FFT 不是近似，而是同一个 DFT 的快速算法。")
print()
print("=== 2. 频率分辨率：谁决定能不能分开两个频率 ===")
rows = []
for fs_v, n in ((8000, 256), (8000, 1024), (44100, 1024), (44100, 4096), (48000, 48000)):
rows.append(["%d Hz" % fs_v, "N = %d" % n, "%.5f s" % (n / fs_v), "%.4f Hz" % (fs_v / n)])
table(["采样率", "点数", "采集时长 T = N/fs", "分辨率 df = fs/N"], rows)
print(" df 只跟采集时长有关：N = 1024 在 8 kHz 下只有 128 ms（df = 7.8125 Hz），")
print(" 而在 44.1 kHz 下更短（23.2 ms，df = 43.07 Hz）—— 想分得更细只能录得更久。")
print()
print("=== 3. 加窗：旁瓣压得越低，主瓣就越宽 ===")
pts = 16384
rows = []
wn = 64
for kind, cn in (("rect", "矩形"), ("hann", "汉宁"), ("hamming", "汉明"), ("blackman", "布莱克曼")):
sp = win_spectrum(kind, wn, pts)
pk = sp[0]
step = math.pi / pts
null = None
for p in range(1, pts):
if sp[p] > sp[p - 1]:
null = p - 1
break
width_bins = 2 * null * step / (2 * math.pi / wn)
sl = max(sp[null:]) / pk
half = sp[pts // wn]
rows.append([cn, "%.2f" % width_bins, "%.2f" % (20 * math.log10(sl)),
"%.4f" % (sum(win_make(kind, wn)) / wn),
"%.2f" % (20 * math.log10(half / pk))])
table(["窗函数", "主瓣总宽（谱线数）", "最大旁瓣 (dB)", "相干增益", "最大扇形损失 (dB)"], rows)
print(" 主瓣宽与旁瓣电平是一对反向指标：矩形窗主瓣最窄（2 条谱线）但旁瓣只压到 -13.25 dB；")
print(" 布莱克曼主瓣宽到 6 条谱线，换来 -58.11 dB 的旁瓣。工程上按'弱信号要被看见'来选。")
print()
print("=== 4. 一个正弦测出来的幅度：扇形损失与相干增益 ===")
fs = 8000.0
n = 1024
rows = []
for kind, cn in (("rect", "矩形"), ("hann", "汉宁"), ("hamming", "汉明"), ("blackman", "布莱克曼")):
w = win_make(kind, n)
cg = sum(w) / n
line = []
for f0 in (1000.0, 1050.0):
sig = [math.cos(2 * math.pi * f0 * i / fs) * w[i] for i in range(n)]
re2, im2 = fft(sig)
mag = [abs(complex(re2[k], im2[k])) for k in range(n // 2 + 1)]
line.append(max(mag) * 2.0 / n)
rows.append([cn, "%.6f" % line[0], "%.6f" % line[1], "%.6f" % (line[1] / cg)])
table(["窗函数", "1000 Hz（整谱线）测幅", "1050 Hz（半谱线）测幅", "1050 Hz 修正后"], rows)
print(" 1000 Hz 正好落在第 128 条谱线上，矩形窗测得 1.000000，一点泄漏都没有；")
print(" 1050 Hz 落在 134.4 号谱线（两条谱线之间），矩形窗只测得 0.756679 —— 这就是扇形损失。")
print(" 加窗后用相干增益修正（汉宁 0.5、汉明 0.54、布莱克曼 0.42），幅度才回到 1 附近。")
print()
print("=== 5. 补零不提高分辨率：只是把谱线画密 ===")
fs = 8000.0
f1, f2 = 1000.0, 1020.0
rows = []
for n, pad_to, tag in ((256, 256, "N = 256 采集"), (256, 4096, "N = 256 采集 + 补零到 4096"),
(1024, 1024, "N = 1024 采集")):
s = [math.cos(2 * math.pi * f1 * i / fs) + math.cos(2 * math.pi * f2 * i / fs) for i in range(n)]
s = s + [0.0] * (pad_to - n)
re3, im3 = fft(s)
mag = [abs(complex(re3[k], im3[k])) for k in range(pad_to // 2 + 1)]
lo, hi = int(900 * pad_to / fs), int(1120 * pad_to / fs)
peaks = [k for k in range(lo + 1, hi) if mag[k] > mag[k - 1] and mag[k] >= mag[k + 1]]
peaks.sort(key=lambda k: -mag[k])
mid = (int(f1 * pad_to / fs) + int(f2 * pad_to / fs)) // 2
if len(peaks) >= 2:
top = sorted(peaks[:2])
dip = mag[mid] / min(mag[top[0]], mag[top[1]])
verdict = "分开" if dip <= 0.5 else ("勉强" if dip <= 0.9 else "没分开")
desc = "%.1f/%.1f Hz，峰间凹陷比 %.4f → %s" % (top[0] * fs / pad_to, top[1] * fs / pad_to, dip, verdict)
else:
desc = "只有一个峰（%.1f Hz）→ 没分开" % (peaks[0] * fs / pad_to)
rows.append([tag, "%.4f Hz" % (fs / n), "%.4f Hz" % (fs / pad_to), desc])
table(["采集方式", "真实分辨率 1/T", "补零后谱线间隔", "1000 与 1020 Hz 的结果"], rows)
print(" 补零把谱线间隔从 31.25 Hz 加密到 1.95 Hz，凹陷比却只从'看不到'升到 0.9992 ——")
print(" 分辨两个频率靠的是采集时长，补零只能把同一条曲线画得更光滑。")
print()
print("=== 6. 定点化：Q15 的溢出怎么发生的 ===")
def wrap16(q):
q &= 0xFFFF
return q - 65536 if q >= 32768 else q
rows = []
for expr, a, b in (("0.5 + 0.25", 0.5, 0.25), ("0.5 + 0.5", 0.5, 0.5), ("0.5 + 0.75", 0.5, 0.75),
("-0.75 - 0.5", -0.75, -0.5)):
qa = int(round(a * 32768))
qb = int(round(b * 32768))
s = qa + qb
sat = 32767 if s > 32767 else (-32768 if s < -32768 else s)
rows.append([expr, "%d + %d = %d" % (qa, qb, s),
"%.6f" % (wrap16(s) / 32768.0), "%.6f" % (sat / 32768.0)])
table(["运算", "累加器（32 位中间值）", "截回 Q15（环绕）", "截回 Q15（饱和）"], rows)
print(" Q15 的可表示范围是 [-1, 1)，真值 1.25 已经出界：")
print(" 环绕得到 -0.750000（符号翻转，听感是爆音），饱和得到 +0.999969（只是削顶）。")
print(" 真实 DSP 的累加器有 32 位甚至 40 位，只在最后一步截回 16 位，就是为了不在这发生溢出。")
python 本站为静态站，不提供在线运行；可复制到本地用 gcc / python 执行
预期输出：
=== 1. FFT 与直接 DFT：同一个答案，两种算法 ===
谱线 实部 Re 虚部 Im 幅度
k = 0 20.000000 0.000000 20.000000
k = 1 -5.828427 -2.414214 6.308644
k = 2 0.000000 0.000000 0.000000
k = 3 -0.171573 -0.414214 0.448342
k = 4 0.000000 0.000000 0.000000
k = 5 -0.171573 0.414214 0.448342
k = 6 0.000000 0.000000 0.000000
k = 7 -5.828427 2.414214 6.308644
两种算法最大偏差 = 7.700e-15（浮点舍入级别，可视为完全相同）
N DFT 复数乘 N^2 FFT 复数乘 (N/2)log2N 加速比
8 64 12 5.33
64 4096 192 21.33
1024 1048576 5120 204.80
8192 67108864 53248 1260.31
结果一模一样，代价差三个数量级 —— FFT 不是近似，而是同一个 DFT 的快速算法。
=== 2. 频率分辨率：谁决定能不能分开两个频率 ===
采样率 点数 采集时长 T = N/fs 分辨率 df = fs/N
8000 Hz N = 256 0.03200 s 31.2500 Hz
8000 Hz N = 1024 0.12800 s 7.8125 Hz
44100 Hz N = 1024 0.02322 s 43.0664 Hz
44100 Hz N = 4096 0.09288 s 10.7666 Hz
48000 Hz N = 48000 1.00000 s 1.0000 Hz
df 只跟采集时长有关：N = 1024 在 8 kHz 下只有 128 ms（df = 7.8125 Hz），
而在 44.1 kHz 下更短（23.2 ms，df = 43.07 Hz）—— 想分得更细只能录得更久。
=== 3. 加窗：旁瓣压得越低，主瓣就越宽 ===
窗函数 主瓣总宽（谱线数） 最大旁瓣 (dB) 相干增益 最大扇形损失 (dB)
矩形 2.00 -13.25 1.0000 -3.92
汉宁 4.00 -31.47 0.5000 -1.42
汉明 4.00 -42.45 0.5400 -1.75
布莱克曼 6.00 -58.11 0.4200 -1.10
主瓣宽与旁瓣电平是一对反向指标：矩形窗主瓣最窄（2 条谱线）但旁瓣只压到 -13.25 dB；
布莱克曼主瓣宽到 6 条谱线，换来 -58.11 dB 的旁瓣。工程上按'弱信号要被看见'来选。
=== 4. 一个正弦测出来的幅度：扇形损失与相干增益 ===
窗函数 1000 Hz（整谱线）测幅 1050 Hz（半谱线）测幅 1050 Hz 修正后
矩形 1.000000 0.756679 0.756679
汉宁 0.500000 0.450492 0.900984
汉明 0.540000 0.474987 0.879605
布莱克曼 0.420000 0.387423 0.922436
1000 Hz 正好落在第 128 条谱线上，矩形窗测得 1.000000，一点泄漏都没有；
1050 Hz 落在 134.4 号谱线（两条谱线之间），矩形窗只测得 0.756679 —— 这就是扇形损失。
加窗后用相干增益修正（汉宁 0.5、汉明 0.54、布莱克曼 0.42），幅度才回到 1 附近。
=== 5. 补零不提高分辨率：只是把谱线画密 ===
采集方式 真实分辨率 1/T 补零后谱线间隔 1000 与 1020 Hz 的结果
N = 256 采集 31.2500 Hz 31.2500 Hz 只有一个峰（1000.0 Hz）→ 没分开
N = 256 采集 + 补零到 4096 31.2500 Hz 1.9531 Hz 1003.9/1017.6 Hz，峰间凹陷比 0.9992 → 没分开
N = 1024 采集 7.8125 Hz 7.8125 Hz 1000.0/1023.4 Hz，峰间凹陷比 0.2799 → 分开
补零把谱线间隔从 31.25 Hz 加密到 1.95 Hz，凹陷比却只从'看不到'升到 0.9992 ——
分辨两个频率靠的是采集时长，补零只能把同一条曲线画得更光滑。
=== 6. 定点化：Q15 的溢出怎么发生的 ===
运算 累加器（32 位中间值） 截回 Q15（环绕） 截回 Q15（饱和）
0.5 + 0.25 16384 + 8192 = 24576 0.750000 0.750000
0.5 + 0.5 16384 + 16384 = 32768 -1.000000 0.999969
0.5 + 0.75 16384 + 24576 = 40960 -0.750000 0.999969
-0.75 - 0.5 -24576 + -16384 = -40960 0.750000 -1.000000
Q15 的可表示范围是 [-1, 1)，真值 1.25 已经出界：
环绕得到 -0.750000（符号翻转，听感是爆音），饱和得到 +0.999969（只是削顶）。
真实 DSP 的累加器有 32 位甚至 40 位，只在最后一步截回 16 位，就是为了不在这发生溢出。五条结论：
- FFT 与直算 DFT 的结果完全一致（偏差
7.7e-15），代价差 1260 倍——加速不是靠近似换来的。 - 分辨率只由采集时长决定：8 kHz 下 1024 点只有 128 ms（
Δf = 7.8125 Hz），44.1 kHz 下同样 1024 点只有 23.2 ms（Δf = 43.0664 Hz），采样率越高、同样的点数分辨率越差。 - 窗函数是"旁瓣 vs 主瓣"的交易：矩形窗 2 条谱线主瓣 / −13.25 dB 旁瓣，布莱克曼 6 条 / −58.11 dB。幅度要除以相干增益才准（汉宁 0.5、汉明 0.54、布莱克曼 0.42）。
- 扇形损失是"半谱线"的代价：1000 Hz 落在整谱线上时矩形窗测得
1.000000，1050 Hz 落在两条谱线之间只测得0.756679（−2.42 dB）；加窗后能收回到0.90左右。 - 补零不会提高分辨率：把谱线间隔从 31.25 Hz 加密到 1.95 Hz，两个相隔 20 Hz 的频率峰间凹陷比仍然只有
0.9992——曲线更光滑，峰还是分不开；真正分开它们的还是更长的采集时间（N = 1024 时凹陷比降到0.2799）。
考点
考点
1. DFT 的定义与三条缩放关系
X[k] = Σ_{n=0}^{N−1} x[n]e^(−j2πnk/N)；Δf = fs/N。- 解析幅度
A的余弦 →|X[k]| = A·N/2；读幅度要乘2/N。 - 实序列
X[k] = X*[N−k]，只有N/2 + 1条独立谱线。
2. FFT 的代价与前提（高频）
- 复数乘：DFT
N²，FFT(N/2)log₂N；复数加：N(N−1)对N log₂N。 - radix-2 要求
N = 2^m，否则先补零到最近的 2 的幂（补零不影响分辨率）。 - 蝶形：
A' = A + W·B、B' = A − W·B，W = e^(−j2πk/N)；共log₂N级、每级N/2个。 - FFT 与 DFT 结果完全相同，不是近似算法。
3. 加窗
- 泄漏的来源：有限长度截断 = 乘矩形窗 = 频域卷积 sinc。
- 主瓣宽 ↔ 旁瓣电平是反向指标：矩形 2 条谱线 / −13 dB，汉宁 4 条 / −31 dB，汉明 4 条 / −42 dB，布莱克曼 6 条 / −58 dB（以谱线数为单位）。
- 相干增益：加窗后峰值要除以它（汉宁 0.5、汉明 0.54、布莱克曼 0.42）。
- 扇形损失：频率落在两条谱线之间时峰值降低；矩形窗最多 3.92 dB，加窗后 ≤ 1.8 dB。
4. 分辨率与补零（易错重灾区）
分辨率
Δf = 1/T，只由采集时长决定。补零只加密谱线（叫"谱线间隔"），不提高分辨率。
- 采样率固定时，点数越多 = 采得越久 = 分辨率越好；但"同样点数、采样率更高"反而更差。
- 要分开间隔
Δ的两个频率，至少需要T ≥ 1/Δ。
5. FIR 与 IIR
| FIR | IIR | |
|---|---|---|
| 结构 | 只有 Σ b_k x[n−k] | 含 Σ a_k y[n−k] 反馈 |
| 线性相位 | h[n] = ±h[N−1−n] 时成立 | 一般做不到 |
| 群延迟（线性相位 FIR） | (N−1)/2 拍 | — |
| 稳定性 | 恒稳定 | 极点全在单位圆内才稳定 |
6. 定点运算
Q15：1 位符号 + 15 位小数，范围[−1, +1)；0.5 = 16384。- 累加器必须比乘数宽（16×16 → 32 位），只在最后一步截回。
- 溢出：环绕会翻转符号（1.25 → −0.75），饱和只削顶（1.25 → 0.999969），工程上必须饱和。
- 定标不当会让"看起来正常"的滤波器在某个输入下突然爆音——这是定点 DSP 最难查的一类 bug。
7. 易错点清单
- 把
Δf = fs/N和"谱线间隔"混为一谈：补零改变后者，不改变前者。 - 忘记
2/N或相干增益，读出来的幅度偏一倍。 - 用偶数长度的窗定义（分母
N−1）与 DFT 型的周期窗（分母N）混用：本页统一用分母N（周期型），因为后者在谱分析里与 DFT 天然对齐。 - 认为 FFT 是"近似算法"：它是精确的分治算法，只是要求长度为 2 的幂。
- 认为"点数越多分辨率越高"而不看采样率：关键量是时长
T。 - 拿加窗后的幅度直接当真实幅度报出去，不做相干增益修正。
小结
- DFT 是"把 N 个采样变成 N 条谱线"，分辨率
Δf = fs/N = 1/T，实序列只有一半独立谱线。 - FFT 是同一个 DFT 的快速算法：
N² → (N/2)log₂N，N = 8192时快 1260 倍，结果逐位一致（本页实测偏差7.7e-15）。 - 泄漏来自截断，加窗是唯一对策，代价是主瓣变宽、幅度要除以相干增益。
- 分辨率由采集时长决定，补零只能把曲线画光滑——这是本页最容易考错的一条。
- FIR 保证线性相位与绝对稳定，IIR 用少数几阶换陡度；线性相位 FIR 的群延迟是
(N−1)/2拍。 - 定点实现里溢出是头号敌人：
Q15的 1.25 环绕后变成 −0.75，这是"能跑但会爆"的典型故障。
回到分支：上一章把 s 换成 z，让整套变换语言在采样系统里成立，并且给出了 H(z)——这一章把 H(z) 变成能在机器上跑的算法：DFT/FFT 负责"看"，差分方程负责"改"，窗函数与定点化负责"在真实的有限精度和有限时长里还能用"。到这里，信号与系统 那四章埋下的工具全部落地了。
往下走有两条路：要设计滤波器（给定指标反求系数、巴特沃斯与窗函数法）去 滤波器设计基础；要对信号做变换与传输（调制、编码、信道）去 通信原理：调制、编码、信道。再往下，采样定理 会回过头把"为什么必须先抗混叠"这件事算清楚——那正是这一章所有结论成立的前提。前置数学仍是 复变函数与积分变换。
这一章回答的是"有一段采样，怎么把它变成频谱、再变成滤波后的序列"；下一章问的是"给定指标，怎么把滤波器造出来"。
下一篇：滤波器设计基础
评论（0）
当前浏览器不允许本地存储，评论无法保存。
还没有评论，来说两句。