10492 你可曾听过位运算的天籁？
.
上周末，我的好基友 Gus 给我发来两行 Matlab 代码，让我非常 excited：
t = 0:65535;
a = bitand(t, bitshift(t, -8));这段代码可以生成一个长度为 65536 的序列 a，其中第 t 个点（从第 0 个开始）的值为 t 的二进制表示中高 8 位与低 8 位的逻辑与，用 C 语言表达则是 t & (t >> 8)。
这个序列有什么特别之处呢？用 plot(t, a, '.') 命令可以画出它的散点图，如下：
一个简简单单的逻辑与运算，竟然产生了分形结构！这个图形就是大名鼎鼎的 Sierpinski triangle（谢尔宾斯基三角形）。
然而这个序列的牛逼之处不止于此 —— 它不仅能看，还能听！用 soundsc(a, 8192) 命令，可以用 8192 Hz 的采样率把这个序列当成声音信号播放出来。你听，层层叠叠的，很有太空感，仿佛一艘飞船由远及近，然后淹没在一片低频振动声中。这段音乐被称为 Sierpinski Harmony（谢尔宾斯基和谐曲）。我第一次把这个声音播放出来的时候，还以为是某个网页弹出了游戏呢！
点此欣赏 谢尔宾斯基和谐曲（原版）：sierpinski-harmony.wav
我在研究这个序列的过程中，曾经一不小心使用了 sound 命令而不是 soundsc 命令来播放。soundsc 命令会在播放前把信号缩放，使得各点的绝对值都不超过 1；而 sound 命令则会粗暴地将绝对值超过 1 的点截成 ±1。序列 a 中的值都是 0 到 255 的整数，所以 sound 命令的效果，就是把所有非零采样点全部置为 1。出乎意料的是，这个截顶操作并没有多么严重地破坏原有的信号；相反，它增强了信号中的高频部分，听起来更像电子游戏了！
点此欣赏 谢尔宾斯基和谐曲（截顶版）：sierpinski-harmony-clipped.wav
为什么几步简单的位运算，可以产生出如此丰富的视觉和听觉效果呢？这篇文章就来解答下面三个问题：
- 为什么位运算 t & (t >> 8) 可以产生分形图案？
- 为什么上述位运算可以产生富有太空感的音乐？
- 为什么截顶之后的信号听起来像电子游戏？
在文末的参考资料中，你还可以发现更多的、由位运算产生出来的奇妙的音乐。
一、位运算的视觉解释
文首提到的那个长度为 65536 的序列 a，可以分成 256 段来分析，每段 256 个点。这个序列的第 i 段（从第 0 段开始），是 0 到 255 的每个数依次与 i 取逻辑与的结果。不难看出，第 0 段的 256 个数均为 0，而第 1 段则是 0, 1, 0, 1 交替出现。但是，在画图的时候，横轴的跨度是很大的，每一段在图中只占很小的宽度。于是，第 0 段的 256 个点就挤在一起成了一个点，而第 1 段的 128 个 0 挤在一起成为一个点，128 个 1 挤在一起成为另一个点。这两段在图的左下角形成了「⣠」的形状。
再看第 2 段和第 3 段。第 2 段是 0 到 255 依次与 2 取与，即取出每个数的倒数第二个二进制位，结果是 0, 0, 2, 2 重复出现。第 3 段是 0 到 255 依次与 3 取与，即取出每个数的最后两个二进制位，相当于对 4 取余，结果是 0, 1, 2, 3 重复出现。画成图以后，第 2 段就挤成了两个点，纵坐标分别 0 和 2；第 3 段则挤成了四个点，纵坐标分别为 0, 1, 2, 3。于是，前 4 段在图的左下角形成了「⣠⣺」的形状。
我们换一种方式来理解第 2、3 段是怎么形成的。第 2、3 段各数的最低（二进制）位，与第 0、1 段各数都相同。再看次低位，第 0、1 段各数的次低位均为 0，但在第 2、3 段中，则有一半为 0，有一半为 1。于是，第 2、3 段中有一半的数跟第 0、1 段是一样的，它们在图中形成了与第 0、1 段相同的图案，即把第 0、1 段的图案向右复制了一份；而另一半则比第 0、1 段中相应的数大 2，它们把第 0、1 段的图案向右上复制了一份。需要注意的是，这两份复制品的「密度」（即每个点是由几个点挤成的）都减半了。
上面由前 2 段的图案复制得到了前 4 段的图案。用同样的思路，可以知道前 4 段到前 8 段、前 8 段到前 16 段、等等等等，都遵循同样的规律 —— 向右、向右上各复制一份，且密度减半。这种重复的过程，就是生成分形图案的关键。复制到 256 段时，我们就得到了完整的「谢尔宾斯基三角形」图案。图案越往右越稀疏，所以能看出一些毛边。到最后一段，每个点其实都只是 1 个点「挤」成的了。这一段是 0 到 255 依次与 255 取与，其实就是 0 到 255 本身，它们组成了一条斜线，从图中其实能看出右边并不是竖直的。
二、位运算的听觉解释
要解释「谢尔宾斯基和谐曲」的听觉效果，关键在于把序列中各个数的每个二进制位拆分开来。下面我们依次分析 a & 1、a & 2、a & 4、……、a & 128 这些序列的听觉效果，而序列 a 则是这八个序列的线性叠加。
序列 a 的全长为 65536 个采样点。由于播放时的采样率为 8192 Hz，整个序列的长度就是 65535 / 8192 = 8 秒。依然把序列 a 拆成 256 段来分析，于是每一段的长度就是 8 / 256 = 1/32 秒。
下面看 a & 1 这个信号。在第 0 段，这个信号一直为 0；在第 1 段，这个信号是 0, 1 交替。在第 2 段，这个信号又一直为 0；在第 3 段，这个信号又是 0, 1 交替。由此看出，这个信号是一个断续的振动：微观上看，它的振动周期是 2 个采样点，频率就是 4096 Hz；宏观上看，它的断续周期是 2 段，即 1/16 s，频率则是 16 Hz。借用通信原理中的术语，这是一个调幅波，载波频率为 4096 Hz，调制频率为 16 Hz。
再看 a & 2。它在第 0、1 段为全零；在第 2、3 段为 0, 0, 2, 2 重复。此后，在第 4、5 段为全零，在第 6、7 段为 0, 0, 2, 2 重复。于是，这是一个微观振动周期为 4 个采样点、断续周期为 4 段的调幅波，或者说 8 Hz 调制的 2048 Hz 载波。
依此类推，可以得到信号 a 的各个二进制位的载波频率和调制频率（表中其它各列将在下文讲解）：
人耳可听频率范围的下限为 20 Hz。频率高于此值的有规律振动，会被人耳感知为乐音；频率低于此值的有规律振动，则会被感知为节奏。观察刚刚求得的载波频率和调制频率，我们发现，载波频率都在 20 Hz 以上，会被感知为乐音；调制频率则都在 20 Hz 以下，会被感知为节奏。
观察各个载波频率，我们发现 a & 16 这个信号的载波频率（256 Hz）与中央 C 的频率（约 261 Hz）非常接近。事实上，在「科学调弦法」中，中央 C 的频率就是被调成 256 Hz 的。为叙述简便，下文中就把 256 Hz 称为中央 C，也记作 C4（4 表示第 4 个八度）。再看 a & 8 这个信号，它的载波频率是 512 Hz，比中央 C 高一个八度，也就是 C5。依此类推，组成「谢尔宾斯基和谐曲」的八个信号，其载波频率分别对应 C1 ~ C8 这八个音，正好横跨了整个钢琴键盘。其中，音高为 C8、C7 这两个信号的声音太尖，人耳不容易捕捉到，所以你听到第一个清晰的音是 C6。
再看各个调制频率。a & 16 的调制频率为 1 Hz，这意味着 C4 这个音符每次持续 0.5 秒。如果采用 120 拍每分钟的速度，0.5 秒恰好就是 1 拍，即四分音符。如果把组成「谢尔宾斯基和谐曲」的八个信号看成是八个声部，那么 C4 所在的声部的乐谱，就是一个四分休止符、一个四分音符这样重复下去。用同样的方法分析所有八个声部，就可以得到「谢尔宾斯基和谐曲」的总谱了（点击可查看大图）。音调越低的音，出现得越晚，每次持续和休止的时间也越长。
我们知道一个乐音有四方面的特征：时长、音高、幅度、音色。上面已经分析了乐谱中每个音的时长和音高，下面分析另外两个特征。先看比较简单的幅度。a & 1 这个信号中，去掉全零的部分，是 0, 1 交替出现；a & 2 这个信号中，去掉全零的部分，是 0, 0, 2, 2 交替出现。不难看出，后者的幅度是前者的 2 倍。依此类推，音高每降低一个八度，幅度都增加一倍。这就解释了为什么到后来所有的音符都淹没在一片低频振动中了。同时我们也发现了 C8、C7 两个声部不容易听到的另一个原因：幅度太小。幅度信息我用力度记号写在了总谱中。
最后分析音色。在时域上看，每个声部的波形都是方波。方波的傅里叶变换比较有意思，它仅含奇次谐波，不含偶次谐波；若以基波的幅度为 1，则第 k 次谐波（k 为奇数）的幅度就是 1/k，可以称为「反比例衰减」。这两个特点有什么特别之处吗？语谱图上可以见分晓。
我用 Audition 软件画出了「谢尔宾斯基和谐曲」的语谱图。窗长设为 256 个采样点，这样可以得到比较平衡的时域和频域分辨率。语谱图（下图的下半部分）的横轴为时间，范围为 0 至 8 秒；纵轴为频率，范围为 0 至 4096 Hz（采样率的一半）；图中各点的颜色表示强度，黑色表示无能量，紫色较弱，红色较强，黄色最强。每段横线的左右两端会有微弱的竖线，那是每个音符开始和结束时引起的信号突变导致的，可以忽略。
语谱图的最上方有一排密密麻麻的小细线，每秒钟有 16 段。它们的频率是 4096 Hz，是由 a & 1 产生的。语谱图的正中，是一排稍长的线段，每秒钟有 8 段。它们的频率是 2048 Hz，是由 a & 2 产生的。注意方波没有偶次谐波，所以这两个信号在频域上互不干扰。然后，在 1024 Hz 频率上，是一排每段持续 0.125 秒的线段，它们是 a & 4 这个信号的基波。在 3072 Hz 频率处，有一排同样的线段，它们是 a & 4 的三次谐波。这两排线段合在一起被感知为 C6 声部。同样地，由于方波没有偶次谐波，a & 4 这个信号在频域上正好避开了 2048 Hz 的 C7 和 4096 Hz 的 C8（这两个频率分别是 a & 4 基频的 2 倍和 4 倍）。依此类推可以发现，由于方波「没有偶次谐波」的特点，所有声部的信号在频域上都是互不干扰的。
语谱图的右半部分，叠加的信号比较多的时候，呈现出从低频到高频由黄到红的非常平滑的渐变。为什么看似独立的八个声部能够叠加出如此完美的渐变呢？我们看 C1、C2 两个声部。C1 的 1、3、5 ……次谐波的幅度之比为 1 : 1/3 : 1/5 ……。C2 的 1、3、5 ……次谐波，可以看成是 C1 的 2、6、10 ……次谐波，它们的幅度之比也是 1 : 1/3 : 1/5 ……。前面分析过，C2 的整体幅度是 C1 的一半。由于它们的音色相同（都是方波），所以对应谐波的幅度之比也是 1:2。于是，C2 的 k 次谐波（可以看成 C1 的 2k 次谐波），其幅度正是 C1 基波幅度的 1/(2k)。这样，把 C2 的谐波插到 C1 的谐波中间，得到的信号仍然符合「反比例衰减」的规律。上面的分析对于任意两个声部都成立，所以整体上就形成了平滑的渐变。这是方波谐波「反比例衰减」的特性和各个声部「幅度倍增」两个因素共同作用的结果。
从 C8 开始的连续若干个声部叠加，在时域上看，就是一群方波叠加成了锯齿波。例如，序列 a 的第 15 段（0.5 秒前的一小段），是由 C8、C7、C6、C5 四个声部叠加而成的；这一段的波形是 0 到 255 依次与 15 取与（即对 16 取余），其值从 0 上升到 15 后回到 0 并循环，确实是锯齿波。锯齿波的傅里叶变换的特点是：奇偶次谐波俱全，反比例衰减。这正好符合上一段的分析。
当多个声部叠加时，其中的高音可以看成是低音的谐波，从而整体上可以当成一个音色不是方波的低音。不过，由于高音在频繁地起起停停，与低音节奏不同，所以人耳还是会倾向于把它们感知为两个声部。
三、截顶操作对声音信号的影响
sound 函数进行的截顶操作，将序列 a 中所有的非零值全部置为 1。这是一个非线性操作，直接分析它对听觉效果的影响，恐怕十分困难。不过我们可以另辟蹊径。所有非零值置 1 这个操作，可以看成是所有二进制位进行了「逻辑或」运算。0、1 两个数之间的「逻辑或」运算用算术运算表达起来比较麻烦，但「逻辑与」则比较好表达 —— 它就是乘法。于是，截顶操作就可以看成是 a 的八个声部各自取反，然后相乘，再取反。取反操作不影响听觉效果，所以截顶操作的本质就是八个声部相乘；它与不截顶的信号的区别就在于，不截顶时八个声部是相加的。
时域中的加法，在频域中仍是加法；时域中的乘法，在频域中则变成了卷积。我们来看截顶后的信号。它的第 0 段是全零，第 1、2 段只有一个声部，这些跟未截顶的信号是一样的。而第 3 段则是 4096 Hz 和 2048 Hz 两个方波（取反）相乘（再取反），在频域上为卷积。这两个方波均只有基波，但不要忘了还有直流分量。卷积的结果，是直流、2048 Hz、4096 Hz 这三个分量均存在。
之后，第 4 段跟未截顶的信号也是一样的。第 5 段是 4096 Hz 与 1024 Hz 两个方波相乘（取反就省略了，下同）。4096 Hz 方波含有直流与 4096 Hz 分量，1024 Hz 方波含有直流、1024 Hz、3072 Hz 三个分量。可以验算，频域卷积之后的信号含有直流、1024 Hz、3072 Hz、4096 Hz 四个分量，但不含 2048 Hz 分量。也就是说，本来不存在的分量没有被卷积出来。
上面这个结论可以推广到任意多个方波在频域卷积。其中，频率最低的方波 x 含有直流和各奇次谐波；其它方波 y 的频率都是它的偶数倍，其谐波可以看成 x 的偶次谐波。x 的奇次谐波与 y 的直流分量或各谐波相卷积，得到的都是 x 的奇次谐波，这些分量本来都是存在的；x 的直流分量与 y 相卷积，得到的就是 y 本身，这也是本来存在的。
以上分析基本解答了这个问题：截顶后的信号各段含有哪些频率分量。上面这个论证其实并不完整，没有证明本来存在的分量不会在卷积过程中相消叠加而消失。不过从截顶后信号的语谱图来看，这种现象并没有发生。
要分析截顶后的信号，还要解答另一个问题：各段中各频率分量的比例如何。这个问题我也没有定量计算，不过观察语谱图之后，我可以给出一些马后炮式的解释：
- 整体上看，信号在各频段的能量分布比较均匀，不像截顶前呈反比例衰减。这是因为，在大多数（四分之三的）段落中，4096 Hz 方波或者 2048 Hz 方波参与了频域卷积。这两个信号的直流分量与基波分量的幅度是相等的，这造成卷积后能量的宏观分布比较均匀。在只有较低频的信号参与频域卷积的段落（如第 64 段，2 秒处开始），就能观察到低频较强，高频较弱了（语谱图上看起来似乎高频也比较强，那其实是第 65 段，4096 Hz 方波参与频域卷积后的结果）。
- 语谱图不像截顶前那样呈现平滑渐变，而是显示出许多水平方向的亮条或亮斑。它们是两个频率相差较大的方波在频域卷积产生的：低频方波的能量集中在直流分量附近，衰减比较快，于是它与高频方波在频域卷积的效果就是拉宽了高频方波各谐波的宽度。低频方波频率非常低的时候，拉宽效果有限，形成亮条；低频方波频率不那么低的时候，拉宽效果明显，形成亮斑。
- 在有许多方波参与频域卷积的段落，语谱图就「糊」掉了，例如最后一段。这可以这么理解：在时域上，许多方波相乘的效果是只剩个别采样点的值为 1 了；极端地，如果是从 4096 Hz 开始的连续若干个方波相乘，会得到冲激串；相乘的方波越多，冲激串越稀疏。冲激串的傅里叶变换还是冲激串，且其密度与原冲激串成反比。冲激串在语谱图上表现为竖直方向等间距、等亮度的亮条（如第 63 段，2 秒处结束），当相乘的方波多了，频域上的冲激串就会变密，呈现出「糊」掉的模样（如第 255 段）。
上面这些特点，可以解释为什么截顶后的信号听起来有电子游戏的感觉。频域能量分布均匀，避免了截顶前低频振动淹没一切的现象，丰富的高频分量使得信号听起来更清脆；语谱图中的亮条和亮斑使得一些乐音听起来更突出；较多方波在时域相乘形成的冲激串产生一种「抖动感」，与早期电子游戏中的爆炸声有相似之处。
四、扩展阅读
文首构造的序列 a，来自 2013 年 SIGINT 会议上的一场演讲「Making Music with a C Compiler」（YouTube 视频）。序列 a 位于视频的 24 分 30 秒处；从 12 分 40 秒至 30 分 35 秒，演讲者 Nils Dagsson Moskopp（网名 erlehmann）从简单到复杂地展示了许多位运算式的听觉效果，最后几个式子产生的音乐的复杂和美妙达到了匪夷所思的程度。这场演讲的主题，是强调用计算机创作音乐不一定需要乐谱，而是可以用一些简单的代码直接生成波形。从这个意义上讲，「谢尔宾斯基和谐曲」并不是一个典型的例子，因为我仍然可以把它的乐谱写出来。
视频的 30 分 35 秒处引用了 arXiv 上 2011 年 12 月的一篇论文 Discovering novel computer music techniques by exploring the space of short computer programs，文中介绍了几种利用位运算生成音乐的套路，并简要分析了它们的原理。当然，大部分的音乐是无法像本文一样拆解到骨头的，许多位运算式也不是有意构造出来，而是在胡乱尝试的过程中误打误撞发现的。作者 Ville-Matias Heikkilä（网名 viznut）在下面两篇博客中记录了他和网友们探索各种位运算式的过程，其中第二篇博客里也有一些原理分析。
Viznut 把他和网友们发现的位运算音乐汇编成了三个视频（均为 YouTube 视频，需翻墙）：
还有一些没有放进视频的位运算式，汇集在了这个文本文件里。
在 viznut 和网友们探索的过程中，他们制作了两个在线工具，可以直接输入位运算式而听到声音：
- http://wurstcaptures.untergrund.net/music/（基于 JavaScript）
- 8-bit Generative Composer（基于 ActionScript）
读者不妨以 viznut 等人发现的算式为起点，探索、创作自己的音乐。
