一個地球人在臺北: 使用 Python 產生聲音及繪製頻譜圖

作者：王一哲
日期：2019/8/13

前言

我在去年看到了啾啾鞋講解電話按鍵聲的影片 [1]，裡面用到的技術稱為雙音多頻訊號 (Dual-Tone Multi-Frequency, DTMF) [2]，各個按鍵對應的聲音頻率如下表所示

	1209 Hz	1336 Hz	1477 Hz	1633 Hz
697 Hz	1	2	3	A
770 Hz	4	5	6	B
852 Hz	7	8	9	C
941 Hz	*	0	#	D

另外也在 YouTube 上找到了成功大學數學系舒宇宸教授演講的影片 [3]，教授在演講時用 MATLAB 產生了上述表格中的聲音。既然用 MATLAB 可以做到，應該也可以用 Python 做到同樣的效果才對，於是我花了一天的時間研究了一下，找到以下的解決方案。

安裝所需套件

PyAudio

產生聲音會用到的套件是 PyAudio [4]，在 Windows 10 可以使用 Windows PowerShell(系統管理員) 輸入以下指令安裝

pip install pyaudio

在 Debian/Ubuntu 則是輸入以下指令安裝

sudo apt-get install python-pyaudio python3-pyaudio

如果沒有辦法抓到最新版的 PyAudio，可以改用以下指令安裝

pip3 install pyaudio

我在 Linux Mint 19.1 Cinnamon 測試過用 apt-get 安裝 pyaudio 的指令，可以成功安裝套件。

PySoundFile

為了將產生的聲音儲存成音效檔，我選用的是 PySoundFile [5]，在 Windows 10 可以使用 Windows PowerShell(系統管理員) 輸入以下指令安裝

pip install soundfile

在 Linux 則是輸入以下指令安裝

pip3 install pysoundfile

我在 Linux Mint 19.1 Cinnamon 測試以上的指令，可以成功安裝套件。

產生聲音並儲存為音效檔

以下的程式碼主要是參考子風的知識庫 [Python] Pyaudio 教學 [6] 以及 Stack Overflow 網站上的文章 [7] 改寫而成。

import numpy as np
import pyaudio
import soundfile as sf

# DTMF 頻率對照表
freq = {'1': [1209, 697], '2': [1336, 697], '3': [1477, 697], 'A': [1633, 697],
        '4': [1209, 770], '5': [1336, 770], '6': [1477, 770], 'B': [1633, 770],
        '7': [1209, 852], '8': [1336, 852], '9': [1477, 852], 'C': [1633, 852],
        '*': [1209, 941], '0': [1336, 941], '#': [1477, 941], 'D': [1633, 941]}

# 定義正弦波
def sine(frequency, duration, rate):
    n = int(duration * rate)
    interval = 2 * np.pi * frequency / rate
    return np.arange(n) / rate, np.sin(np.arange(n) * interval) 

# 產生兩個特定頻率的聲音資料
def generator(f1=697, f2=1209, duration=1, rate=44100, start=0):
    t1, d1 = sine(f1, duration, rate)
    t2, d2 = sine(f2, duration, rate)
    data = np.stack((t1+start, d1, d2), axis=1)
    return data

# 產生寫入 stream 的資料
names = ['0', '1', '2', '3', '4', '5', '6', '7',
         '8', '9', '*', '#', 'A', 'B', 'C', 'D']
data = []
duration = 0.2
for i in range(len(names)):
    tmp = generator(f1=freq[names[i]][0], f2=freq[names[i]][1], 
                    duration=duration, start=i*duration)
    data = np.append(data, tmp)

data = data.reshape(int(len(data)/3), 3)

# 產生聲音
p = pyaudio.PyAudio()
stream = p.open(format=pyaudio.paFloat32, channels=2, rate=44100, output=True)
stream.write(data[:, 1:].astype(np.float32).tostring())
stream.stop_stream()
stream.close()
p.terminate()

# 儲存資料檔
np.savetxt('DTMF_data.csv', data, delimiter=',')

# 儲存 wav 檔
sf.write('DTMF_sound.wav', data[:, 1:], 44100, subtype='PCM_24')

第 6 ~ 9 行：用字典格式儲存每個按鍵對應的聲音頻率。
第 12 ~ 15 行：定義正弦波，輸入項目為聲音頻率 frequency (Hz)、播放時間長度 duration (s)、取樣頻率 rate (Hz)，回傳項目為時刻及振幅。函式中 n 為播放數量， interval 為取樣間隔，由於 sin 裡面要放的是 $\theta = \omega t$ ，因此取樣間隔定義為 $(2 \pi f) / rate = \omega / rate$。
第 18 ~ 22 行：產生兩個特定頻率的聲音資料，輸入項目為兩個指定的頻率 f1 及 f2、播放時間長度 duration、取樣頻率 rate、開始播放時刻 start，回傳項目為時刻以及兩個頻率聲波對應的振幅。
第 25 ~ 34 行：產生寫入 stream 的資料。
第 37 ~ 42 行：產生聲音，執行程式時會聽到 0 ~ 9 以及 *、#、A、B、C、D 對應的聲音。
第 45 行：儲存資料檔，之後要畫頻譜圖。
第 48 行：儲存 wav 檔。

繪製頻譜圖

假設剛才將資料儲存成文字檔 DTMF_data.csv，我們可以用以下的程式碼繪製頻譜圖。

import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np

# 讀取資料
t, d1, d2 = np.loadtxt('DTMF_data.csv', dtype=float, delimiter=',', unpack=True)

# 產生繪圖物件
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15,8), dpi=72, sharex=True)

# 時刻 - 振幅關係圖
amp = d1 + d2
ax1.plot(t, amp, linestyle='-', color='blue')
ax1.set_xlabel('time (s)', fontsize=14)
ax1.set_ylabel('amplitude', fontsize=14)

# 頻譜圖
cmap = mpl.cm.cool
Fs = int(1.0 / (t[1]-t[0]))
data, freqs, bins, im = ax2.specgram(amp, cmap=cmap, NFFT=1024, Fs=Fs, noverlap=128)
ax2.set_ylabel('frequency (Hz)', fontsize=14)
ax2.axis(ymax=2000)
fig.colorbar(im, ax=ax2, shrink=0.6, orientation='horizontal', pad=0.2)

# 儲存及顯示圖片
fig.savefig('DTMF_specgram.png')
fig.show()

第 6 行：從 DTMF_data.csv 讀取資料，將資料按照欄位依序存入陣列 t、d1、d2。
第 9 行：產生繪圖物件，兩張小圖使用同樣的 x 軸。
第 12 ~ 15 行：繪製時刻 - 振幅關係圖。
第 18 ~ 23 行：使用 numpy.specgram 繪製頻譜圖，以下是簡單的語法，詳細的語法請參考官方說明書 [8] 及範例 [9]。

numpy.specgram([資料], 
                   cmap=[Colormaps 設定], 
                   NFFT=[計算 FFT 時每個區塊中的資料點數量], 
                   Fs=[取樣頻率],
                   noverlap=[兩個區塊重疊的資料點數量])

第 26 ~ 27 行：儲存及顯示圖片。

以下是輸出的圖片。由於資料點太密集，時刻 - 振幅關係圖的線條全部擠在一起，看不出什麼特徵。從頻譜圖可以看出資料分為16個區塊，每個區塊中看起來有兩個頻率的強度較強，如果將 y 軸最大值限制為 2000，可以更清楚地看到這些頻率，我們再跟文章最前面的 DTMF 表格對照看看，應該可以看出每一個區塊對應的按鍵。

上方為時刻 - 振幅關係圖，下方為頻譜圖（不限制 y 軸範圍）

上方為時刻 - 振幅關係圖，下方為頻譜圖（y 軸最大值為 2000）

結語

PyAudio 的功能很多，除了產生聲音之外，可以用來讀取、處理音效檔，也可配合麥克風錄音。之後應該可以用 PyAudio 和 Raspberry Pi 做一些有趣的專案。

參考資料

HackMD 版本連結：https://hackmd.io/@yizhewang/By8Jb36mS

一個地球人在臺北

置頂

我的 VPython 教學文件 (HackMD 版本)

熱門文章

2019年8月13日星期二

使用 Python 產生聲音及繪製頻譜圖

前言

安裝所需套件

PyAudio

PySoundFile

產生聲音並儲存為音效檔

繪製頻譜圖

結語

參考資料

沒有留言:

張貼留言

置頂

我的 VPython 教學文件 (HackMD 版本)

熱門文章

2019年8月13日 星期二

使用 Python 產生聲音及繪製頻譜圖

前言

安裝所需套件

PyAudio

PySoundFile

產生聲音並儲存為音效檔

繪製頻譜圖

結語

參考資料

沒有留言:

張貼留言

2019年8月13日星期二