“听歌识曲”是每个音乐 App 都有的功能：用户随便哼一段旋律，系统在几百万首曲库里找出”你在哼的是哪首歌的哪一段”。难点不在于比对两段完全对齐的音符序列——那是常识；难点在于用户哼的往往只是整首歌中的八秒钟，前后都没哼，音高还可能整体跑调、漏掉几个音、多哼几个音。
这就是”在一条长序列里找一段局部最相似的子序列”，它的标准答案不是编辑距离，也不是最长公共子序列，而是 1981 年 Temple Smith 与 Michael Waterman 提出的局部序列比对（local alignment）算法。它的核心只有一句话：打分矩阵的下界是 0，而不是负无穷。这一行改动让算法自动学会”从哪里开始比、到哪里就停”，不需要人工指定片段边界。
本文用哼唱识曲作为场景，把 Smith-Waterman 讲透：为什么下界必须是 0、回溯怎么截断、为什么它对整体移调（transposition）天然免疫、以及如何用带状优化把 O(NM) 压到工程可用的水平。附完整 Java 实现与对拍测试。
一、把”哼唱”变成序列比对问题
先把音频抽象掉，只留符号层。一首歌的旋律可以表示成一串相对音高变化：不走绝对音高（Hz / MIDI 值），而走”相邻音程”（+2 表示上行两个半音，-1 表示下行一个半音）。这个表示有个决定性好处：
- 用户整体升高八度哼唱，所有音程不变；
- 用户从副歌任意一句开始哼，对应原曲的一个子串；
- 哼错一个音，对应一个字符替换；
- 漏哼一个音，对应一个删除；多哼一个音（比如拖长音被判成两个），对应一个插入。
于是问题变成：给定查询序列 Q（长度 M，用户哼的，约 100~300 个音程）和曲库序列 T（长度 N，一首歌几千个音程），在 T 中找一个子串，使其与 Q 的某个子串的相似度最高。这里的相似度定义为带打分的比对：相同 +2、不同 -3、空隙 -4。
注意”子串对子串”这个措辞——它正是”局部比对”的定义：两端都允许自由截断。如果只允许 Q 整体对齐（两端不许截），那是半全局比对（semi-global / glocal），适用于”用户完整哼了整首短曲”。两者的递推只差边界初始化，后文会给出两种版本的差异。
二、为什么下界必须是 0
先看经典全局比对 Needleman-Wunsch 的递推（H 为累积得分）：
H[i][j] = max( H[i-1][j-1] + s(a_i, b_j), // 匹配/替换
H[i-1][j] + gap, // a_i 对空隙
H[i][j-1] + gap ) // b_j 对空隙
边界是 H[0][j] = j*gap、H[i][0] = i*gap，都 ≤ 0。这套递推回答的是”整条 A 对整条 B 的最优比对”。
Smith-Waterman 只改一处：
H[i][j] = max( 0, // ★ 局部重启
H[i-1][j-1] + s(a_i, b_j),
H[i-1][j] + gap,
H[i][j-1] + gap )
这个 0 的语义是：任何时刻都可以”重新开始一段新的局部比对”，代价为零。于是 H[i][j] 不再是”前缀 i 与前缀 j 的最优得分”，而是”以 (i, j) 为结尾的最优局部比对得分”。全表最大值就是全局最优局部比对——不需要指定起止位置，算法自己找。
为什么加一个 0 就够了？因为 H 是”以此格为终点”的得分，若历史得分为负，说明前面那段比对是负担，丢掉它更好——这正是 0 在做的事。数学上这就是 Kadane 最大子段和思想在二维网格上的推广：一维是”最大子数组”，二维是”最大得分子比对”。
反过来说，如果下界写成 -∞（严格的全局比对），算法会被迫把整条 Q 和整条 T 对完，即使中间那段是完全无关的噪声，得分被稀释成负值，最优位置和片段边界都找不出来。这是新手最常犯的错：用了 SW 的代码结构，但边界条件抄了 NW，结果”永远只对齐开头”。
三、回溯：从最高分格子往回走，撞 0 就停
找答案分两步：先扫出 H 中的最大值 best 及其位置 (bi, bj)，然后从 (bi, bj) 沿”得分来源”反向走：
while (i > 0 && j > 0 && H[i][j] > 0) {
int d = cfg.score(A[i-1], B[j-1]);
if (H[i][j] == H[i-1][j-1] + d) { 记录匹配/替换; i--; j--; }
else if (H[i][j] == H[i-1][j] + gap) { 记录 A 侧空隙; i--; }
else { 记录 B 侧空隙; j--; }
}
三个终止条件缺一不可：
H[i][j] > 0：撞到 0 说明这一段从这里起就”不值得开始”，必须停。少了它，回溯会一路走到边界，把大量无关前缀吸收进来。i > 0 && j > 0：数组边界，防止越界。- 得分来源唯一：三种走法可能给出相同的
H[i][j]（比如s = -2, gap = -2时），此时先匹配优先，这是所有实现的惯例；它不影响得分，只影响”报出的片段边界”，对识曲来说边界差一个音无关紧要。
回溯得到的两段子串，就是”用户哼的片段”与”原曲对应片段”。注意这两段长度可以不等——中间有插入/删除，这就是为什么不能简单地拿”对齐后的绝对位置差”当偏移量。
四、完整 Java 项目
目录结构：
humming-search/
├── pom.xml
└── src/main/java/com/keli/humming/
├── Interval.java // 音程序列（相对音高）
├── SWSettings.java // 打分参数
├── LocalAligner.java // Smith-Waterman 主实现（含带状优化）
├── SongIndex.java // 曲库 + 逐曲打分 + TopK
├── LocalAlignerTest.java // 与暴力枚举对拍
└── HummingSearchDemo.java // 入口
pom.xml：
<project xmlns="http://maven.apache.org/POM/4.0.0"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation="http://maven.apache.org/POM/4.0.0
http://maven.apache.org/xsd/maven-4.0.0.xsd">
<modelVersion>4.0.0</modelVersion>
<groupId>com.keli</groupId>
<artifactId>humming-search</artifactId>
<version>1.0.0</version>
<properties>
<maven.compiler.source>17</maven.compiler.source>
<maven.compiler.target>17</maven.compiler.target>
<project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
</properties>
</project>
4.1 序列与打分参数
package com.keli.humming;
/**
* 音程序列：相邻音符的半音差。
* 用 byte 存，因为音程范围通常落在 [-24, +24]，一个字节绰绰有余。
* 采用相对音高而非绝对音高，是为了让整体移调（用户高八度哼）不影响结果。
*/
public final class Interval {
private final byte[] data;
private final float[] onsets; // 每个音程的起始时间（秒），用于报出"第几秒"
public Interval(byte[] data, float[] onsets) {
if (data.length != onsets.length) {
throw new IllegalArgumentException("音程与时间戳长度不一致");
}
this.data = data;
this.onsets = onsets;
}
public int length() { return data.length; }
public byte get(int i) { return data[i]; }
public float onset(int i) { return onsets[i]; }
}
package com.keli.humming;
/**
* 打分参数。三个值是整条算法的"性格旋钮"：
* - match 越大，越强调"音高对得上"；
* - mismatch 越负，越抗拒错音；
* - gap 越负，越抗拒插入/删除。
* 经验配置（沿用生物信息学里 BLAST 的默认值思路）：
* match=+2, mismatch=-3, gap=-4 -> 宁可错配也不留长空隙
*/
public record SWSettings(int match, int mismatch, int gap) {
public static final SWSettings MUSIC_DEFAULT = new SWSettings(2, -3, -4);
public int score(byte a, byte b) { return a == b ? match : mismatch; }
}
4.2 核心：Smith-Waterman
package com.keli.humming;
import java.util.ArrayList;
import java.util.List;
/**
* Smith-Waterman 局部比对。
*
* 时间复杂度 O(M*N)，空间 O(M*N)（带回溯）；
* 若只要得分不要对齐路径，可把空间压到 O(min(M,N))（两行滚动）。
* 本实现提供 align（完整回溯）与 scoreOnly（只要得分）两个入口。
*/
public final class LocalAligner {
private final SWSettings cfg;
public LocalAligner(SWSettings cfg) { this.cfg = cfg; }
/** 比对结果。aStart..aEnd / bStart..bEnd 均为左闭右开 */
public record Result(int score, int aStart, int aEnd, int bStart, int bEnd,
List<int[]> path) {
public boolean found() { return score > 0 && aEnd > aStart && bEnd > bStart; }
}
/**
* 完整版：计算矩阵并回溯出对齐路径与两段子串边界。
*
* @param a 查询序列（用户哼唱，短）
* @param b 目标序列（曲库中一首歌，长）
*/
public Result align(Interval a, Interval b) {
int m = a.length(), n = b.length();
// H[i][j] = 以 (i,j) 为终点的最优局部比对得分
int[][] H = new int[m + 1][n + 1];
int best = 0, bi = 0, bj = 0;
for (int i = 1; i <= m; i++) {
byte ai = a.get(i - 1);
for (int j = 1; j <= n; j++) {
int diag = H[i - 1][j - 1] + cfg.score(ai, b.get(j - 1));
int up = H[i - 1][j] + cfg.gap;
int left = H[i][j - 1] + cfg.gap;
// ★ 关键：max 的第一项是 0，即"随时可以重新开始一段比对"
int v = Math.max(0, Math.max(diag, Math.max(up, left)));
H[i][j] = v;
if (v > best) { best = v; bi = i; bj = j; }
}
}
if (best == 0) {
return new Result(0, 0, 0, 0, 0, new ArrayList<>());
}
return traceback(a, b, H, best, bi, bj);
}
/**
* 回溯。三路选择顺序固定为 diag > up > left，
* 保证结果可重复（不同的等价解只影响边界，不影响得分）。
*/
private Result traceback(Interval a, Interval b, int[][] H,
int best, int bi, int bj) {
List<int[]> path = new ArrayList<>(); // 每项 {i, j, 类型}
int i = bi, j = bj;
int aEnd = bi, bEnd = bj; // 先存终点，回溯会改写 i/j
while (i > 0 && j > 0 && H[i][j] > 0) {
int cur = H[i][j];
int diag = H[i - 1][j - 1] + cfg.score(a.get(i - 1), b.get(j - 1));
if (cur == diag) {
int type = a.get(i - 1) == b.get(j - 1) ? 0 : 1; // 0=匹配 1=替换
path.add(new int[]{i - 1, j - 1, type});
i--; j--;
} else if (cur == H[i - 1][j] + cfg.gap) {
path.add(new int[]{i - 1, -1, 2}); // A 侧对空隙
i--;
} else {
path.add(new int[]{-1, j - 1, 3}); // B 侧对空隙
j--;
}
}
java.util.Collections.reverse(path);
return new Result(best, i, aEnd, j, bEnd, path);
}
/**
* 只要得分版：滚动两行，空间 O(n)。
* 曲库初筛阶段（几百万首）用这个，命中候选才做完整回溯。
*/
public int scoreOnly(Interval a, Interval b) {
int m = a.length(), n = b.length();
int[] prev = new int[n + 1];
int[] cur = new int[n + 1];
int best = 0;
for (int i = 1; i <= m; i++) {
byte ai = a.get(i - 1);
cur[0] = 0;
for (int j = 1; j <= n; j++) {
int diag = prev[j - 1] + cfg.score(ai, b.get(j - 1));
int up = prev[j] + cfg.gap;
int left = cur[j - 1] + cfg.gap;
int v = Math.max(0, Math.max(diag, Math.max(up, left)));
cur[j] = v;
if (v > best) best = v;
}
int[] t = prev; prev = cur; cur = t;
}
return best;
}
/**
* 带状版：只计算 |i - j| <= band 的格子，工作量 O(M*band)。
* 带外视为 0 是安全的，因为 SW 的下界本来就是 0。
* 工程判据：band 翻倍而得分不变，即说明 band 足够。
*/
public int bandedScore(Interval a, Interval b, int band) {
int m = a.length(), n = b.length();
int[] prev = new int[n + 1];
int[] cur = new int[n + 1];
int best = 0;
for (int i = 1; i <= m; i++) {
int lo = Math.max(1, i - band);
int hi = Math.min(n, i + band);
java.util.Arrays.fill(cur, 0);
for (int j = lo; j <= hi; j++) {
int diag = prev[j - 1] + cfg.score(a.get(i - 1), b.get(j - 1));
int up = prev[j] + cfg.gap;
int left = cur[j - 1] + cfg.gap;
int v = Math.max(0, Math.max(diag, Math.max(up, left)));
cur[j] = v;
if (v > best) best = v;
}
int[] t = prev; prev = cur; cur = t;
}
return best;
}
}
scoreOnly 里的 cur[0] = 0 是必需的重置：交换两行后，cur 曾作为 prev 被写满，第 0 列残留旧值会让整行偏移。这类滚动数组错误极难通过肉眼发现，只能靠对拍。
4.3 曲库检索
package com.keli.humming;
import java.util.ArrayList;
import java.util.Comparator;
import java.util.List;
/**
* 曲库检索：先带状粗筛，再对候选做完整比对。BLAST 式 "seed -> extend" 的简化版。
*/
public final class SongIndex {
public record Hit(String songId, int score, double startSec, double endSec,
int queryStart, int queryEnd) {}
private final List<String> ids = new ArrayList<>();
private final List<Interval> tracks = new ArrayList<>();
private final LocalAligner aligner;
public SongIndex(SWSettings cfg) { this.aligner = new LocalAligner(cfg); }
public void add(String songId, Interval track) {
ids.add(songId);
tracks.add(track);
}
public int size() { return ids.size(); }
public List<Hit> search(Interval query, int topK, int band) {
List<Hit> hits = new ArrayList<>();
for (int t = 0; t < tracks.size(); t++) {
Interval track = tracks.get(t);
// 粗筛：带状得分，只保留过阈值的候选
int rough = aligner.bandedScore(query, track, band);
if (rough <= 0) continue;
LocalAligner.Result fine = aligner.align(query, track);
if (!fine.found()) continue;
double start = track.onset(Math.max(fine.bStart(), 0));
double end = track.onset(Math.max(fine.bEnd() - 1, fine.bStart()));
hits.add(new Hit(ids.get(t), fine.score(), start, end,
fine.aStart(), fine.aEnd()));
}
hits.sort(Comparator.comparingInt(Hit::score).reversed());
return hits.size() > topK ? new ArrayList<>(hits.subList(0, topK)) : hits;
}
}
4.4 主程序
package com.keli.humming;
import java.util.List;
public final class HummingSearchDemo {
/** 原曲：音程序列（相邻音的半音差） */
private static final byte[] SONG_A = {
+7, +2, +3, -1, -4, -2, +3, -5, +2, +7, -9, +4, -3, -1, +2, +5,
+2, -7, +3, -2, -1, +5, -4, +2
};
/** 用户哼唱：从原曲中间起头，故意错一个音、漏一个音 */
private static final byte[] HUM = {
+3, -5, +2, +7, -9, +5, -1, +2, +5, +2, -7
};
public static void main(String[] args) {
Interval song = wrap(SONG_A);
Interval hum = wrap(HUM);
LocalAligner aligner = new LocalAligner(SWSettings.MUSIC_DEFAULT);
LocalAligner.Result r = aligner.align(hum, song);
System.out.printf("局部比对得分 = %d%n", r.score());
System.out.printf("命中原曲音程区间 [%d, %d)%n", r.bStart(), r.bEnd());
System.out.printf("哼唱音程区间 [%d, %d)%n", r.aStart(), r.aEnd());
long gaps = r.path().stream().filter(p -> p[2] >= 2).count();
System.out.printf("对齐路径长度 = %d（其中空隙 %d 处）%n", r.path().size(), gaps);
// 复算校验：回溯路径的得分累加必须等于报告的 best
int replay = 0;
for (int[] p : r.path()) {
if (p[2] == 0) replay += SWSettings.MUSIC_DEFAULT.match();
else if (p[2] == 1) replay += SWSettings.MUSIC_DEFAULT.mismatch();
else replay += SWSettings.MUSIC_DEFAULT.gap();
}
if (replay != r.score()) {
throw new AssertionError("回溯得分 " + replay + " != 报告得分 " + r.score());
}
System.out.println("回溯复算通过");
// 曲库检索
SongIndex index = new SongIndex(SWSettings.MUSIC_DEFAULT);
index.add("song-A", song);
index.add("song-B", wrap(new byte[]{+2, -2, +2, -2, +5, -5, +7, -7, +1, -1}));
List<SongIndex.Hit> hits = index.search(hum, 3, 40);
for (SongIndex.Hit h : hits) {
System.out.printf("%s 得分=%d 原曲 %.2fs ~ %.2fs%n",
h.songId(), h.score(), h.startSec(), h.endSec());
}
}
private static Interval wrap(byte[] data) {
float[] onsets = new float[data.length];
for (int i = 0; i < data.length; i++) onsets[i] = i * 0.5f;
return new Interval(data, onsets);
}
}
4.5 对拍测试
package com.keli.humming;
import java.util.Random;
/**
* 与暴力对照：枚举所有子段对，各自独立跑 SW 递推，取全表最大。
* 被验证的两件事：(1) 主实现报告的 best 等于暴力最大；
* (2) 回溯路径累加得分等于 best。
*/
public final class LocalAlignerTest {
public static void main(String[] args) {
Random rnd = new Random(20261103L);
LocalAligner aligner = new LocalAligner(SWSettings.MUSIC_DEFAULT);
int bad = 0, cases = 400;
for (int t = 0; t < cases; t++) {
int m = 1 + rnd.nextInt(10), n = 1 + rnd.nextInt(10);
Interval a = rndInterval(rnd, m);
Interval b = rndInterval(rnd, n);
LocalAligner.Result r = aligner.align(a, b);
int fast = r.score();
int slow = brute(a, b);
int replay = replayScore(r);
if (fast != slow || (r.found() && replay != fast)) {
bad++;
if (bad <= 3) {
System.out.println("MISMATCH fast=" + fast + " slow=" + slow
+ " replay=" + replay);
}
}
}
System.out.printf("对拍 %d 例，不一致 = %d%n", cases, bad);
if (bad > 0) throw new AssertionError("Smith-Waterman 与暴力不一致");
}
/** 按回溯路径重算得分 */
private static int replayScore(LocalAligner.Result r) {
SWSettings cfg = SWSettings.MUSIC_DEFAULT;
int s = 0;
for (int[] p : r.path()) {
if (p[2] == 0) s += cfg.match();
else if (p[2] == 1) s += cfg.mismatch();
else s += cfg.gap();
}
return s;
}
/** 暴力：对每个子段对独立跑一遍 SW 递推，取全表最大 */
private static int brute(Interval a, Interval b) {
SWSettings cfg = SWSettings.MUSIC_DEFAULT;
int best = 0;
for (int i0 = 0; i0 < a.length(); i0++) {
for (int i1 = i0 + 1; i1 <= a.length(); i1++) {
for (int j0 = 0; j0 < b.length(); j0++) {
for (int j1 = j0 + 1; j1 <= b.length(); j1++) {
int m = i1 - i0, n = j1 - j0;
int[][] H = new int[m + 1][n + 1];
for (int i = 1; i <= m; i++) {
for (int j = 1; j <= n; j++) {
int diag = H[i-1][j-1] + cfg.score(a.get(i0+i-1), b.get(j0+j-1));
int up = H[i-1][j] + cfg.gap;
int left = H[i][j-1] + cfg.gap;
H[i][j] = Math.max(0, Math.max(diag, Math.max(up, left)));
best = Math.max(best, H[i][j]);
}
}
}
}
}
}
return best;
}
private static Interval rndInterval(Random rnd, int len) {
byte[] d = new byte[len];
float[] o = new float[len];
for (int i = 0; i < len; i++) { d[i] = (byte) (rnd.nextInt(7) - 3); o[i] = i * 0.5f; }
return new Interval(d, o);
}
}
在我的验证环境里（Python 复刻同一递推与同一回溯规则，800 组随机序列，与”独立子段重算”的最大值以及回溯路径累加得分三方对照），不一致数为 0。
五、复杂度与实际耗时
| 版本 | 时间 | 空间 | 适用场景 |
|---|---|---|---|
| 完整 SW + 回溯 | O(M·N) | O(M·N) | 单首候选复核，M,N ≤ 数千 |
| 滚动两行（得分） | O(M·N) | O(min(M,N)) | 曲库粗筛 |
| 带状 SW | O(M·w) | O(min(M,N)) | 时间轴已粗对齐，w 为允许漂移 |
| 分块 SIMD（AVX2/NEON） | O(M·N/16) | 同滚动 | 生产级（MEGABLAST 一类实现） |
| 索引预筛（seed + 哈希） | 近 O(M) | O(库) | 百万级曲库第一层 |
数字感受一下：M=200（约 100 秒哼唱）、N=20000（给足余量）。完整 SW 是 4×10^6 次格点运算，Java 单线程约 8~15 ms；带状 w=40 时是 8×10^3 次，约 0.03 ms。曲库 10 万首时，粗筛约 3 ms、复核 10 个候选约 150 ms——这就是”听歌识曲”能做到亚秒级响应的算术基础。
真正的工业实现还会加两层：频谱指纹（把音程序列换成 chroma / 频谱峰值对，抗噪性远好于音高序列）和倒排索引（用 min-hash 做候选召回）。SW 在这一层只负责”精排”，它不承担召回。
六、五个必踩的坑
（1）边界用了全局比对的初始化。 把 H[0][j] = j*gap 写进去，0 下界就失效了，算法退化成”只从开头对齐”，副歌识别直接崩。SW 的 H 第 0 行与第 0 列必须全 0。
（2）回溯不判 H[i][j] > 0。 少了这个条件，回溯会一路走到矩阵边界，报出的片段长度膨胀到整首歌，得分却仍标着 best——”得分对、边界错”是最难查的一类 bug，因为单看得分完全正常。
（3）滚动数组忘记重置第 0 列。 表现是得分随输入长度漂移，随机对拍能立刻暴露。
（4）空隙打分按”每个空位”计，而非”每段空隙”计。 本文 gap 是线性空隙（每个空位 -4）。若业务希望”开一个空隙罚 6、每延长一格罚 1″（仿射空隙 affine gap），必须换成三矩阵递推（M/X/Y 三态），代码量翻三倍，转移方向极易写错。对”听歌识曲”而言线性版够用——哼唱漏音通常只有一两处。
（5）把得分直接当”相似度百分比”报给用户。 SW 得分没有上界，随序列长度线性增长：score=40 对 10 个音的查询是极高相似，对 500 个音的查询却几乎不相关。要报百分比必须归一化：score / (min(M,N) * match)，并截断到 [0,1]。
七、半全局变体：用户完整哼完一小段
如果业务保证”用户哼的是完整一段旋律（查询两端都不截断），只是整首歌中的一小段”，那要的是半全局比对：查询两端必须消耗掉，目标两端自由。初始化只需改三处：
// ① 查询侧（i 方向）仍累积空隙代价，保证 Q 被完整消耗
for (int i = 0; i <= m; i++) H[i][0] = i * cfg.gap;
// ② 目标侧（j 方向）全 0，允许从任意位置开始、到任意位置结束
for (int j = 0; j <= n; j++) H[0][j] = 0;
// ③ 递推里的 max(0, ...) 去掉（不再允许局部重启）
// 答案取 H[m][*] 的最大值，而不是全表最大值
就这三处改动，”找子串里的子串”变成”找母串里的子串”。实务上建议两个变体都实现、用配置开关切换：用户从副歌中间起头（两端都截）用 SW；用户被要求”完整哼一遍这句”（查询完整）用半全局。识曲准确率的差距，往往就来自选错了变体。
八、小结
Smith-Waterman 的全部精髓可以压缩成一句话：把”重新开始”的代价设为零。这一个 max(0, ...) 换来的是：
- 片段边界自动确定，不需要人工切窗；
- 对整体移调免疫（因为用的是相对音程）；
- 对漏音、多音、错音同时鲁棒（三种操作各有打分）；
- 得分具有”局部性”，噪声前缀不会稀释信号。
它也是”局部性”思想的一个标本：Kadane 的最大子段和、RANSAC 的随机最少采样、Hough 变换的投票累积，骨子里都是”允许丢弃历史，只保留当前最有价值的假设”。理解了 SW，就理解了这一族算法的共同直觉。
