2021-12-21

音乐在频域的可视化(频谱)

这个小项目是我上学期学数字信号处理的时候想搞的，当时觉得音乐播放器里跳动的频谱很酷炫，就想着自己实现一个，于是就有了用matlab做的。但是效果不太满意，于是这几天重新用python写了一个。先给大家看看效果。

这里是效果动图

原理讲解

从正弦波到方波

正弦波大家都见过，就是sin, cos，这里统称正弦波，那么有没有想过，方波是什么样子的？如果你学过傅里叶变换，那么下意识的反应，方波是频率不同的正弦波的叠加。但是你要没有学过也没有关系，简单的来说，不同的正弦波(有不同的频率，不同的幅度，不同的相位，正弦波的三要素)可以组合出来方波。

这里是傅里叶动图

任意信号的分解

看过了上面的方波例子，信号其实可以被分解为频率不同的正弦波的叠加，然后你将频率看作一个轴(将f看作x轴)，然后不同的频率的正弦波信号所占总信号的多少是不一样，这个多少就用垂直于频率轴的一个轴可以进行表示(将他看作y轴)，然后你就有了平时在音乐播放器中所看见的频谱。而这个频谱就是我们要可视化的东西。

如何将时域的信号变成频域信号

我们在上信号与系统这门课，或者说在接触傅里叶这个人之前，我们一般是以时间为自变量理解音频，即y = f(t), t表示时间，y表示在某个时间的音量的强度，f表示对应关系，也就是哪个时间是哪个声音的强度。这就是我们所说的时域信号。那么什么是频域信号，其实差不多，就是z = g(f), f是频率，z是某个频率信号所占的多少，也可以理解为和某个频率的正弦波的相似度，越相似，说明所占的分量越大。

那么，如何将f(t)变成g(f)呢(即f->g)，这里用到的是就是傅里叶变换，其实我更愿意将傅里叶变换看作是一个函数，只不过是函数的函数。他的输入是时域信号(函数)，输出是频域信号(函数)，所以g(f) = F[ f(t) ]。

关于上面的函数的函数，其实有一门叫做泛函分析的数学课里面讲，只不过我现在还没有学。

其实上面的描述并不准确，但是对于大家理解应该有帮助，其实，z是一个复数，z的模值其实代表的是能量的大小，所以，可视化频谱需要的是z的模值。

这里是苹果切片函数

如何用计算机实现这一切

相信大家都听过这句话，计算机是二进制的，所以计算机只能处理离散的数据，所以将我们上面的东西迁移一下，就出现了数字信号处理这门课，于是我们今天所要做的东西便有了理论支撑。在离散的世界中，我们的时域音频信号变成了一串离散数据，其实还是函数，只不过输入变成了数据的序号，输出是对应下标的单个数据。就好比是c语言中学的数组，输入时数组的下标，输出是对应下标的值。又好比是字典，输入是单词，输出是单词的意思，字典本身就是一串单词的意思。

然后我们的频域信号，也就是我们所要可视化的东西，也变成了一串数据，此时下标代表的其实是频率。那么将时域变成频域的函数，就变成了fft(快速傅里叶变换)。

这里是黑客帝国效果的图

关于函数

关于函数这个话题其实能说的还挺多的，这篇博客就不写了，后面有机会专门开一个说说。

代码实现

语言选择

本来打算选matlab的，可能是我的matlab代码水平不行，好像程序执行效率不高。然后就选python了，虽然跟c语言不能比，但是方便啊，这次所用的好多工具都有对应的包，就不用重复造轮子了，好耶！比如fft这个函数就在scipy这个包中。

这里是python图

问题分解

这次要实现的东西主要实现音频播放时进行实时计算，然后将实时频域的图像绘制出来。对应此次任务就可以分解为三个部分，

1.将音频播放出来，并可以实时读取播放的音频时域数据 2.将拿到数据进行实时fft 3.将fft的数据绘制出来

任务一

这里使用sounddevice这个库，这个库中有两种播放音频的方式。一种是直接播放，直接调用play这个函数就可以，但是我们无法拿到实时的音频数据；一种是通过流(Stream)的方式进行播放，然后就可以拿到实时的音频数据。你可以将音频想象成水，然后喇叭在河的下游，然后你将水(音频流)不断倒到河道中就可以让喇叭发出声音。

代码我贴下面：

# audio stream callback function, called when chunk was palyed
def callback(indata, outdata, frames, time, status):
    frames = chunk_size # push chunk_size datas to stream every time
    global sound_chunk_idx

    # judge whether the rest of audio data is enough for fft
    if (music_data[chunk_size*sound_chunk_idx:chunk_size*\
            (sound_chunk_idx+1), 0].size < chunk_size):
        raise sd.CallbackStop()
        return
    # push to stream
    outdata[:] = music_data[chunk_size*sound_chunk_idx:chunk_size*(sound_chunk_idx+1), :]
    sound_chunk_idx += 1;
    # change buffer, for genrating animation
    buffer[:] = np.abs(fft(outdata[0:1024, 0]))[:512:8]

# open stream
with sd.Stream(callback=callback, samplerate=sr, channels=2, \
        blocksize=chunk_size, dtype=music_data.dtype):
    print("Enter q or Q to exit:")
    while True:
        response = input()
        if response in ('', 'q', 'Q'):
            break

可能有人不理解这个callback函数，其实就是通过callback函数完成上面所说的倒水这一过程(将音频推入流中)，这个流会周期性的调用callback这个函数，每次调用，我们倒入一点水(音频数据)，就可以播放音乐。

这里是数据流的图

任务二

这步比较简单，就是将时域数据通过fft变成频域数据，对应的代码在任务一的代码14行。你会看到，这里有个buffer变量，它是一个全局变量，保存频域的实时数据，这个是给后面绘图用的。然后还有个abs函数，这个函数其实是刚好对应了原理讲解中z是复数，所以要取模值，获取能量大小。然后在取完模值之后进行重采样(resample)，就是为了获取更少的数据，因为全部画出来太吃性能了。

这里是fft的图

任务三

这个任务就是绘图了，使用了matplotlib这个绘图库。基本思路就是让matplotlib调用函数不断绘图，绘图的依据就是前面提到过buffer，通过播放时不断更新buffer中的值就可以变成实时频谱图了。

下面是代码

fig, ax = plt.subplots(figsize=(18,9)) # the object of drawing
ax.axis('off') # close axies and border

# init function, first frame
def init():
    ax.set_xlim(1, 65) # x axies
    ax.set_ylim(0, 600000) # y axies
    return ax.stem(buffer, buffer)

# Update function, called at every animtion frame
def update(frame):
    xdata = np.arange(1, 65) # totally 65 data points
    ydata = buffer # change buffer via callback function
    b_c = ax.stem(xdata, ydata, linefmt='--', markerfmt=None) # draw stem plot
    return b_c # must

ani = FuncAnimation(fig, update, init_func=init, interval=20, blit=True)
# animation on
plt.show(block=False)