无所不能的NumPy:我用它弹奏出了CD音质的吉他名曲“爱的罗曼史”

1 前言

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

NumPy几乎无所不能。NumPy最广为人知的能力是图像处理,而它最基础的应用是科学计算。我曾经在《C/C++/Java/Go/Rust,Python喊你来打擂:3秒钟内统计出小于1亿的素数个数》一文中应用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 π ] [0,1176\pi] [0,1176π]之间均匀采样66150点,求其正弦值,逐点连线,即可得到频率为196Hz、持续时间3秒、量化精度16位、采样频率为22.05kHz的正弦波。

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> x = np.linspace(0, 1176*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弦空弦音')

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()

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

  1. 声音幅度在振动中逐渐变小
  2. 频谱显示存在较强幅度的二倍频、三倍频、四倍频
>>> 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 12 \sqrt[12]{2} 122 。知道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')
已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 代码科技 设计师:Amelia_0503 返回首页
实付 15.20元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值