1 前言

大家好，欢迎来到 Crossin的编程教室！

在Python的世界里，没有一个模块能够像NumPy那样支撑并影响着整个生态系统：从科学计算到数据处理，从视觉识别到机器学习，从神经网络到虚拟现实，处处都有它的身影。无论是OpenCV、OpenGL，还是Pandas、Matplotlib，抑或是Scikti-learn、TensorFlow、Keras、Theano、PyTorch，无不依赖于NumPy，尤其是依赖它所创造的数组对象（numpy.ndarray）。NumPy几乎无所不能。NumPy最广为人知的能力是图像处理，而它最基础的应用是科学计算。

本文则是独辟蹊径，讨论如何使用NumPy发出声音，以及如何模拟吉他音色，最终生成CD级别的音乐文件。除了生成wave文件时用到标准模块wave，全部代码仅依赖NumPy单个模块。此外，采集吉他音频数据（供频谱分析用）和播放wave文件时，还用到了PyAudio模块。如果有同学想重复频谱分析、使用代码播放wave文件，请使用如下的命令安装PyAudio模块。

pip install PyAudio

2 发出单一频率的声音

想用NumPy弹奏吉他，先得让NumPy发出声音，比如，生成频率为196Hz（吉他3弦空弦音的频率）且持续3秒钟的声音。这需要3x196=588个周期的波形。假定采样频率为22.05kHz（相当于调频收音机的音质），3秒钟会采集到66150个数据。如果每个数据用两个字节的整数表示，则声音幅度的动态范围从-32768到32767。

以下代码在 [0, 1176π] 之间均匀采样66150点，求其正弦值，逐点连线，即可得到频率为196Hz、持续时间3秒、量化精度16位、采样频率为22.05kHz的正弦波。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(0, 1158*np.pi, 66150, endpoint=False)
>>> y = np.sin(x) * 32767
>>> y.dtype
dtype('float64')
>>> y = y.astype(np.int16)
>>> y.dtype
dtype('int16')
>>> plt.plot(x, y, color='green')
>>> plt.show()

下图左侧画出了全部66150个点，图像变成绿色块了；右侧只画出了前1000个点，可以看出清晰的正弦波形。

现在，声波有了，怎样才能听到声音呢？用PyAudio模块可以直接在声卡上播放这个数据（后面的代码有实现），也可以借助于Python的标准模块wave将声波数据保存为wave文件。wave模块是一个读写后缀名为.wav文件的工具模块，用起来非常简单。

>>> import wave
>>> with wave.open(r'd:\sound_196Hz_3s.wav', 'wb') as fp:
        fp.setparams((1, 2, 22050, 0, 'NONE', 'NONE'))
        fp.writeframes(y.tobytes())

去D盘下看看，多了一个名为sound_196Hz_3s.wav的文件。点击并播放它，就会听到类似轮船开动时的汽笛声。

3 模拟吉他音色

刚才生成的声音，虽然和吉他3弦的空弦音频率相同，但听起来却一点不像吉他。那么吉他发出的声音是什么样的呢？下图是我用声卡采集的吉他1弦和2弦的空弦音的波形图，时长大约4秒钟左右。

下面是采集声音的代码。有兴趣的同学可以拿出吉他或者其他乐器录制一段，然后做一下频谱分析。

import pyaudio
import numpy as np
import time

def capture(rate, chunk):
    pa = pyaudio.PyAudio()
    stream = pa.open(
        format = pyaudio.paInt16,       # 设置量化精度（每个采样数据占用的位数）
        channels = 1,                   # 设置单声道模式           
        rate = 22050,                   # 设置采样频率
        frames_per_buffer = 2205,       # 设置声卡读写缓冲区     
        input = True                    # 设置声卡输出模式
    )

    data = list()
    while len(data) < 50:
        data.append(np.frombuffer(stream.read(2205), dtype=np.int16))
    stream.close()
    pa.terminate()

    return np.hstack(data)

if __name__ == '__main__':
    for i in range(5):
        print(5-i)
        time.sleep(1)
    print('Start')

    data = capture(22050, 1024)
    np.save('吉他10.npy', data)

NumPy提供了快速傅里叶分析工具，可以分析刚才采集到的吉他声音的频率成分及其强度。

import numpy as np
import matplotlib.pyplot as plt
import time
import wave

plt.rcParams['font.sans-serif'] = ['FangSong']
plt.rcParams['axes.unicode_minus'] = False

d1 = np.load('吉他10.npy')
d2 = np.load('吉他20.npy')

plt.subplot(221)
plt.plot(d1, c='g')
plt.title('吉他1弦空弦音')
plt.subplot(222)
plt.plot(d2, c='m')
plt.title('吉他2弦空弦音')
#plt.show()

fd1 = np.fft.fft(d1[22050:44100]) # 截取吉他1弦空弦音1秒钟的数据进行傅里叶分析
fd2 = np.fft.fft(d2[22050:44100]) # 截取吉他2弦空弦音1秒钟的数据进行傅里叶分析

plt.subplot(223)
plt.plot(np.abs(fd1[:11025]/22050), c='g')
plt.title('吉他1弦单边频谱图')
plt.subplot(224)
plt.plot(np.abs(fd2[:11025]/22050), c='m')
plt.title('吉他2弦单边频谱图')
plt.show()

不难看出，吉他发出的声音有以下两个特点：

声音幅度在振动中逐渐变小
频谱显示存在较强幅度的二倍频、三倍频、四倍频

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(0, 3*np.pi, 22050, endpoint=False)
>>> y1 = 1-x/(10*np.pi)+(1-x/(6*np.pi))*np.sin(x)*0.5
>>> x = np.arange(3*22050)/22050
>>> y2 = 0.7*np.exp(-x)
>>> GUITAR_EFFECT_ARRAY = np.hstack((y1, y2))
>>> y_guitar = y*GUITAR_EFFECT_ARRAY[:y.shape[0]]
>>> plt.subplot(121)
>>> plt.plot(GUITAR_EFFECT_ARRAY)
>>> plt.subplot(122)
>>> plt.plot(y_guitar)
>>> plt.show()

上面的代码，生成了一个模拟吉他波形的包络线，将等幅的正弦波模拟成吉他波形。效果如下图所示。

解决了声波的振幅衰减，还需要叠加二倍频、三倍频、四倍频。对于NumPy来说，这都是小菜一碟。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> plt.rcParams['font.sans-serif']=['FangSong']
>>> plt.rcParams['axes.unicode_minus']=False
>>> x1 = np.linspace(0, 4*np.pi, 300)
>>> x2 = np.linspace(0, 8*np.pi, 300)
>>> x3 = np.linspace(0, 12*np.pi, 300)
>>> y1 = np.sin(x1)
>>> y2 = np.sin(x2)
>>> y3 = np.sin(x3)
>>> y = np.sum(np.dstack((y1,y2,y3))[0], axis=1)
>>> plt.subplot(221)
>>> plt.title('基波')
>>> plt.plot(y1, c='be')
>>> plt.subplot(222)
>>> plt.title('二倍频')
>>> plt.plot(y2, c='cn')
>>> plt.subplot(223)
>>> plt.title('三倍频')
>>> plt.plot(y3, c='y')
>>> plt.subplot(224)
>>> plt.title('复合波形')
>>> plt.plot(y, c='m')
>>> plt.show()

不同频率的波形叠加在一起，效果如下图所示。

4 弹出吉他上所有的音符

写到这里，必须要感谢中国古代的一位贵族——明朝皇室子弟朱载堉。他是一位被严重低估了的天才，他也是我认为最类似以解决数学问题为乐趣的欧洲贵族的中国人之一。他的主要学术成果是十二平均律，经过传教士利玛窦传到西方后，直接导致了西方现代音乐理论的兴起。朱载堉同时也是手工开方的先驱，他自制81档算盘开12次方，堪称前无古人后无来者。

扯远了。如果理解十二平均律的话，很容易计算出吉他上所有品格的音频频率。简单说，同一根弦上，从低到高所有品格的频率是一个等比数列，比例为2的1/12次方。知道A=440Hz的话，也很容易根据十二平均律计算出从1弦到6弦的频率分别为329.6Hz、246.9Hz、196.0Hz、146.8Hz、110.0Hz、82.4Hz。

这里定义1弦空弦用字符串'10'表示，1弦3品用字符串'13'表示，字符串'212'则表示2弦12品，字符串'0'则表示休止。计算吉他上指定弦品音频频率的代码如下。

def get_frequency(pos):
    """返回指定弦品pos的频率"""

    fs = (329.6, 246.9, 196.0, 146.8, 110.0, 82.4)
    if pos[0] == '0':
        return 0
    else:
        return fs[int(pos[0])-1] * pow(2, int(pos[1:])/12)

5 吉他谱的格式约定

仅有吉他弦品的表示法还不够，还需要定义每个音符的时值表示。这里约定使用节拍数表示。如果以4分音符为1拍，则4表示全音符，2表示2分音符，1表示4分音符，0.5表示8分音符，0.25表示16分音符，依次类推。

有了吉他弦品和时值的表示，如何描述吉他谱呢？这里约定，一个吉他谱是若干乐谱小节的有序集合，每个乐谱小节又是多个声部的有序集合，每个声部又是一根琴弦上品格和时值的元组的有序集合。这里说的有序集合，就是Pytho的列表。下面是根据这个约定翻译过来的“天空之城”的序曲部分。

castle_in_the_sky= [
    [ # 第1节
        [('10',0.5),('12',0.5)] # 1弦
    ],
    [ # 第2节
        [('13',1.5),('12',0.5),('13',1),('13',0.5),('17',0.5)], # 1弦
        [('0',1),('20',3)], # 2弦
        [('0',0.5),('30',2),('30',1.5)], # 3弦
        [('60',2),('60',2)] # 6弦
    ],
    [ # 第3节
        [('12',4)], # 1弦
        [('0',1),('23',2),('20',1)], # 2弦
        [('0',0.5),('32',1),('32',2.5)], # 3弦
        [('40',2),('40',2)] # 4弦
    ],
    [ # 第4节
        [('10',2),('10',1),('13',1)], # 1弦
        [('21',1.5),('23',2.5)], # 2弦
        [('0',1),('30',1.5),('30',1),('30',0.5)], # 3弦
        [('0',0.5),('42',3.5)], # 4弦
        [('53',4)] # 5弦
    ],
    [ # 第5节
        [('23',3),('20',1)], # 2弦
        [('0',1),('30',1),('30',2)], # 3弦
        [('0',0.5),('40',1),('40',2.5)], # 4弦
        [('52',4)] # 5弦
    ],
    [ # 第6节
        [('0',2.5),('13',1.5)], # 1弦
        [('21',1.5),('20',0.5),('21',2)], # 2弦
        [('0',1),('30',2.5),('30',0.5)], # 3弦
        [('0',0.5),('42',2.5),('42',1)], # 4弦
        [('50',2),('50',2)] # 5弦
    ],
    [ # 第7节
        [('0',3),('13',0.5),('13',0.5)], # 1弦
        [('20',2),('20',2)], # 2弦
        [('0',1),('30',3)], # 3弦
        [('0',0.5),('42',1),('42',2.5)], # 4弦
        [('60',4)] # 6弦
    ],
    [ # 第8节
        [('12',3),('12',0.5),('10',0.5)], # 1弦
        [('0',1.5),('22',0.5),('22',2)], # 2弦
        [('0',1),('33',1.5),('33',1.5)], # 3弦
        [('0',0.5),('44',3.5)], # 4弦
        [('62',2),('62',2)] # 6弦
    ],
    [ # 第9节
        [('12',2),('12',2)], # 1弦
        [('0',1.5),('24',0.5),('24',2)], # 2弦
        [('0',1),('34',1),('34',2)], # 3弦
        [('0',0.5),('44',1.5),('44',2)], # 4弦
        [('51',4)] # 5弦
    ]
]

6 弹奏吉他谱

有了以上这些储备，写一段弹奏吉他谱的代码就是水到渠成了。全部代码不足100行，重点位置已有注释。

import numpy as np
import wave
import pyaudio

SPEED = 80 # 用每分钟节拍数表示弹奏速度
FRAME_RATE = 44100 # 采样速率（44100为CD音质，22050为调频广播音质）
STEREO = True # 立体声（双声道）

# 生成吉他音色包络线
x = np.linspace(0, 3*np.pi, 2*int(FRAME_RATE*60/SPEED), endpoint=False)
y1 = 1 - x/(10*np.pi) + (1-x/(6*np.pi))*np.sin(x)*0.5
x = np.arange(6*int(FRAME_RATE*60/SPEED))/int(FRAME_RATE*60/SPEED)
y2 = 0.7*np.exp(-x)
GUITAR_EFFECT_ARRAY = np.hstack((y1, y2))

def get_frequency(pos):
    """返回指定弦品pos的频率"""

    fs = (329.6, 246.9, 196.0, 146.8, 110.0, 82.4)
    if pos[0] == '0':
        return 0
    else:
        return fs[int(pos[0])-1] * pow(2, int(pos[1:])/12)

def get_wave(f, beat):
    """返回指定频率和节拍数的波形数据"""

    data = list()
    duration = beat*60/SPEED
    sample_num = int(duration*FRAME_RATE)

    for k, p in [(1,0.4), (2,0.3), (3,0.2), (4,0.1)]:
        x = np.linspace(0, 2*duration*f*k*np.pi, sample_num, endpoint=False)
        y = np.sin(x)*p
        data.append(y)

    return guitar_effect(np.sum(np.dstack(data)[0], axis=1))

def guitar_effect(data):
    """将等幅声波变成吉他音色的声波数据"""

    return data*GUITAR_EFFECT_ARRAY[:data.shape[0]]

def play(melody, wave_file=None):
    """弹奏吉他谱，若wave_file存在，同时生成.wav文件"""

    data = list()
    for section in melody:
        data_section = list()
        for cord in section:
            data_cord = list()
            for pos, beat in cord:
                f = get_frequency(pos)
                dw = get_wave(f, beat)
                data_cord.append(dw)
            data_cord = np.hstack(data_cord)
            data_section.append(data_cord)

        d = data_section[0]
        for i in range(1, len(data_section)):
            if d.shape[0] > data_section[i].shape[0]:
                d[:data_section[i].shape[0]] += data_section[i]
            else:
                data_section[i][:d.shape[0]] += d
                d = data_section[i]
        data.append(d)

    data = np.hstack(data)
    data = data*20000/data.max()
    data = data.astype(np.int16)

    if STEREO:
        blank = np.zeros(int(0.006*FRAME_RATE), dtype=np.int16)
        d_left = np.hstack((data, blank))
        d_right = np.hstack((blank, data))
        data = np.dstack((d_left, d_right))[0].ravel()

    if wave_file:
        with wave.open(wave_file, 'wb') as fp:
            fp.setparams((int(STEREO)+1, 2, FRAME_RATE, 0, 'NONE', 'NONE'))
            fp.writeframes(data.tobytes())

    pa = pyaudio.PyAudio()
    stream = pa.open(
        format = pyaudio.paInt16,   # 设置量化精度：每个采样数据占用的位数
        channels = int(STEREO)+1,   # 设置通道数量           
        rate = FRAME_RATE,          # 设置采样频率
        frames_per_buffer = 1024,   # 设置声卡读写缓冲区     
        output = True               # 设置声卡输出模式
    )

    for i in range(0, data.shape[0], 1024):
        stream.write(data[i:i+1024].tobytes())

    stream.stop_stream()
    stream.close()
    pa.terminate()

用哪一首曲子来体验一下代码的效果呢？爱的罗曼史！这首曲子第二段的指法实在变态，对于手指僵硬的我来说，这辈子是不可能弹出第二段的。现在终于有机会弹出我所理解的这首世界名曲了。下面是根据约定翻译过来的吉他名曲“爱的罗曼史”的乐谱。

romance= [
    [
        [('17',1),('17',1),('17',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('17',1),('15',1),('13',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('13',1),('12',1),('10',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('10',1),('13',1),('17',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('112',1),('112',1),('112',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('112',1),('110',1),('18',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('18',1),('17',1),('15',1)],
        [('0',0.33),('25',1),('25',1),('25',0.66)],
        [('0',0.66),('35',1),('35',1),('35',0.33)],
        [('50',3)]
    ],
    [
        [('15',1),('17',1),('18',1)],
        [('0',0.33),('25',1),('25',1),('25',0.66)],
        [('0',0.66),('35',1),('35',1),('35',0.33)],
        [('50',3)]
    ],
    [
        [('17',1),('18',1),('17',1)],
        [('0',0.33),('27',1),('27',1),('27',0.66)],
        [('0',0.66),('38',1),('38',1),('38',0.33)],
        [('67',3)]
    ],
    [
        [('111',1),('18',1),('17',1)],
        [('0',0.33),('27',1),('27',1),('27',0.66)],
        [('0',0.66),('38',1),('38',1),('38',0.33)],
        [('67',3)]
    ],
    [
        [('17',1),('15',1),('13',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('13',1),('12',1),('10',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('60',3)]
    ],
    [
        [('12',1),('12',1),('12',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('32',1),('32',1),('32',0.33)],
        [('52',3)]
    ],
    [
        [('12',1),('13',1),('12',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('32',1),('32',1),('32',0.33)],
        [('52',3)]
    ],
    [
        [('10',1),('10',1),('10',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('30',1),('30',1),('30',0.33)],
        [('42',1),('52',1),('63',1)]
    ],
    [
        [('10',3)],
        [('20',3)],
        [('30',3)],
        [('60',3)]
    ],
    [
        [('14',1),('14',1),('14',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('31',1),('31',1),('31',0.33)],
        [('60',3)]
    ],
    [
        [('14',1),('12',1),('10',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('31',1),('31',1),('31',0.33)],
        [('60',3)]
    ],
    [
        [('25',1),('24',1),('24',1)],
        [('0',0.33),('32',1),('32',1),('32',0.66)],
        [('0',0.66),('44',1),('44',1),('44',0.33)],
        [('50',3)]
    ],
    [
        [('24',1),('23',1),('24',1)],
        [('0',0.33),('32',1),('32',1),('32',0.66)],
        [('0',0.66),('44',1),('44',1),('44',0.33)],
        [('50',3)]
    ],
    [
        [('19',1),('19',1),('19',1)],
        [('0',0.33),('27',1),('27',1),('27',0.66)],
        [('0',0.66),('38',1),('38',1),('38',0.33)],
        [('67',3)]
    ],
    [
        [('19',1),('111',1),('19',1)],
        [('0',0.33),('27',1),('27',1),('27',0.66)],
        [('0',0.66),('38',1),('38',1),('38',0.33)],
        [('67',3)]
    ],
    [
        [('19',1),('17',1),('17',1)],
        [('0',0.33),('29',1),('29',1),('29',0.66)],
        [('0',0.66),('39',1),('39',1),('39',0.33)],
        [('60',3)]
    ],
    [
        [('17',1),('19',1),('111',1)],
        [('0',0.33),('29',1),('29',1),('29',0.66)],
        [('0',0.66),('39',1),('39',1),('39',0.33)],
        [('60',3)]
    ],
    [
        [('112',1),('112',1),('112',1)],
        [('0',0.33),('29',1),('29',1),('29',0.66)],
        [('0',0.66),('39',1),('39',1),('39',0.33)],
        [('60',3)]
    ],
    [
        [('112',1),('111',1),('110',1)],
        [('0',0.33),('29',1),('29',1),('29',0.66)],
        [('0',0.66),('39',1),('39',1),('39',0.33)],
        [('60',3)]
    ],
    [
        [('19',1),('19',1),('19',1)],
        [('0',0.33),('25',1),('25',1),('25',0.66)],
        [('0',0.66),('36',1),('36',1),('36',0.33)],
        [('50',3)]
    ],
    [
        [('19',1),('17',1),('15',1)],
        [('0',0.33),('25',1),('25',1),('25',0.66)],
        [('0',0.66),('36',1),('36',1),('36',0.33)],
        [('50',3)]
    ],
    [
        [('14',1),('14',1),('14',1)],
        [('0',0.33),('24',1),('24',1),('24',0.66)],
        [('0',0.66),('32',1),('32',1),('32',0.33)],
        [('52',3)]
    ],
    [
        [('14',1),('15',1),('12',1)],
        [('0',0.33),('24',1),('24',1),('24',0.66)],
        [('0',0.66),('32',1),('32',1),('32',0.33)],
        [('52',3)]
    ],
    [
        [('10',1),('10',1),('10',1)],
        [('0',0.33),('20',1),('20',1),('20',0.66)],
        [('0',0.66),('31',1),('31',1),('31',0.33)],
        [('42',1),('52',1),('64',1)]
    ],
    [
        [('10',3)],
        [('20',3)],
        [('30',3)],
        [('60',3)]
    ]
]

激动人心的时刻来了，运行下面的代码，美妙的旋律就会从你的声卡中源源不断地飞出来。

play(romance)

若想同时保存成文件，请使用下面的命令。

play(romance, '爱的罗曼史.wav')

大功告成！来听下最终效果：

还不错吧？

完整代码都在文中了，感兴趣的同学可以自己复制到代码中试一下，也可以改成你喜欢的曲目。

如果文章对你有帮助，欢迎转发/点赞/收藏~

作者：天元浪子

来源：Python作业辅导员

_往期文章推荐_

谁在用python弹奏一曲《菊花台》

如需了解付费精品课程及教学答疑服务

请在Crossin的编程教室内回复: 666

用NumPy“弹奏”一首CD级音质的《爱的罗曼史》