<style>
.markdown-body table{
display: unset;
}
</style>
# 使用 Python 產生聲音及繪製頻譜圖
> 作者：王一哲
> 日期：2019/8/13
## 前言
我在去年看到了啾啾鞋講解電話按鍵聲的影片[[1](https://youtu.be/48xx6H-M4ZI)]，裡面用到的技術稱為**雙音多頻訊號** (Dual-Tone Multi-Frequency, **DTMF**)[[2]( https://en.wikipedia.org/wiki/Dual-tone_multi-frequency_signaling)]，各個按鍵對應的聲音頻率如下表所示
<div style="text-align:center">
| | 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 |
</div>
<br />
另外也在 YouTube 上找到了成功大學數學系舒宇宸教授演講的影片[[3](https://youtu.be/475jWptKBtM?t=3192)]，教授在演講時用 MATLAB 產生了上述表格中的聲音。既然用 MATLAB 可以做到，應該也可以用 Python 做到同樣的效果才對，於是我花了一天的時間研究了一下，找到以下的解決方案。
<br />
## 安裝所需套件
### PyAudio
產生聲音會用到的套件是 **PyAudio** [[4]( http://people.csail.mit.edu/hubert/pyaudio/)]，在 Windows 10 可以使用 **Windows PowerShell(系統管理員)** 輸入以下指令安裝
```shell
pip install pyaudio
```
在 Debian/Ubuntu 則是輸入以下指令安裝
```shell
sudo apt-get install python-pyaudio python3-pyaudio
```
如果沒有辦法抓到最新版的 PyAudio，可以改用以下指令安裝
```shell
pip3 install pyaudio
```
我在 Linux Mint 19.1 Cinnamon 測試過用 apt-get 安裝 pyaudio 的指令，可以成功安裝套件。
<br />
### PySoundFile
為了將產生的聲音儲存成音效檔，我選用的是 **PySoundFile** [[5](https://pysoundfile.readthedocs.io/en/0.9.0/)]，在 Windows 10 可以使用 **Windows PowerShell(系統管理員)** 輸入以下指令安裝
```shell
pip install soundfile
```
在 Linux 則是輸入以下指令安裝
```shell
pip3 install pysoundfile
```
我在 Linux Mint 19.1 Cinnamon 測試以上的指令，可以成功安裝套件。
<br />
## 產生聲音並儲存為音效檔
以下的程式碼主要是參考**子風的知識庫 [Python] Pyaudio 教學** [[6]](https://zwindr.blogspot.com/2017/03/python-pyaudio.html) 以及 **Stack Overflow** 網站上的文章 [[7](https://stackoverflow.com/questions/45025892/sound-generated-not-being-saved-to-a-file-as-it-should)] 改寫而成。
<br />
```python=
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')
```
<br />
1. 第 6 ~ 9 行：用字典格式儲存每個按鍵對應的聲音頻率。
2. 第 12 ~ 15 行：定義正弦波，輸入項目為**聲音頻率** frequency (Hz)、**播放時間長度** duration (s)、**取樣頻率** rate (Hz)，回傳項目為**時刻**及**振幅**。函式中 n 為**播放數量**， interval 為**取樣間隔**，由於 sin 裡面要放的是 $\theta = \omega t$ ，因此取樣間隔定義為 $(2 \pi f) / rate = \omega / rate$。
3. 第 18 ~ 22 行：產生兩個特定頻率的聲音資料，輸入項目為**兩個指定的頻率** f1 及 f2、**播放時間長度** duration、**取樣頻率** rate、**開始播放時刻** start，回傳項目為**時刻以及兩個頻率聲波對應的振幅**。
4. 第 25 ~ 34 行：產生寫入 stream 的資料。
5. 第 37 ~ 42 行：產生聲音，執行程式時會聽到 0 ~ 9 以及 \*、\#、A、B、C、D 對應的聲音。
6. 第 45 行：儲存資料檔，之後要畫頻譜圖。
7. 第 48 行：儲存 wav 檔。
<br />
## 繪製頻譜圖
假設剛才將資料儲存成文字檔 DTMF_data.csv，我們可以用以下的程式碼繪製頻譜圖。
```python=
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()
```
<br />
1. 第 6 行：從 DTMF_data.csv 讀取資料，將資料按照欄位依序存入陣列 t、d1、d2。
2. 第 9 行：產生繪圖物件，兩張小圖使用同樣的 x 軸。
3. 第 12 ~ 15 行：繪製時刻 - 振幅關係圖。
4. 第 18 ~ 23 行：使用 **numpy.specgram** 繪製頻譜圖，以下是簡單的語法，詳細的語法請參考官方說明書 [[8](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.specgram.html)] 及範例 [[9](https://matplotlib.org/3.1.1/gallery/images_contours_and_fields/specgram_demo.html#sphx-glr-gallery-images-contours-and-fields-specgram-demo-py)]。
```python
numpy.specgram([資料],
cmap=[Colormaps 設定],
NFFT=[計算 FFT 時每個區塊中的資料點數量],
Fs=[取樣頻率],
noverlap=[兩個區塊重疊的資料點數量])
```
5. 第 26 ~ 27 行：儲存及顯示圖片。
<br />
以下是輸出的圖片。由於資料點太密集，時刻 - 振幅關係圖的線條全部擠在一起，看不出什麼特徵。從頻譜圖可以看出資料分為16個區塊，每個區塊中看起來有兩個頻率的強度較強，如果將 y 軸最大值限制為 2000，可以更清楚地看到這些頻率，我們再跟文章最前面的 DTMF 表格對照看看，應該可以看出每一個區塊對應的按鍵。
<br />
<img height="100%" width="100%" src="https://imgur.com/rIne5Oi.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">上方為時刻 - 振幅關係圖，下方為頻譜圖（不限制 y 軸範圍）</div><br />
<img height="100%" width="100%" src="https://imgur.com/RbcBQM5.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">上方為時刻 - 振幅關係圖，下方為頻譜圖（y 軸最大值為 2000）</div><br />
## 結語
PyAudio 的功能很多，除了產生聲音之外，可以用來讀取、處理音效檔，也可配合麥克風錄音。之後應該可以用 PyAudio 和 Raspberry Pi 做一些有趣的專案。
<br />
## 參考資料
1. [電話按鍵聲的秘密，原來不只是響好聽的? | 超邊緣冷知識 第35集 | 啾啾鞋](https://youtu.be/48xx6H-M4ZI)
2. [WikiPedia DTMF]( https://en.wikipedia.org/wiki/Dual-tone_multi-frequency_signaling)
3. [探索10-8講座：醫學影像中的數學 / 舒宇宸助理教授](https://youtu.be/475jWptKBtM?t=3192)
4. [PyAudio 官方說明書]( http://people.csail.mit.edu/hubert/pyaudio/)
5. [PySoundFile 官方說明書](https://pysoundfile.readthedocs.io/en/0.9.0/)
6. [子風的知識庫 [Python] Pyaudio 教學](https://zwindr.blogspot.com/2017/03/python-pyaudio.html)
7. [Stack Overflow](https://stackoverflow.com/questions/45025892/sound-generated-not-being-saved-to-a-file-as-it-should)
8. [matplotlib.pyplot.specgram 官方說明書](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.specgram.html)
9. [matplotlib.pyplot.specgram 範例](https://matplotlib.org/3.1.1/gallery/images_contours_and_fields/specgram_demo.html#sphx-glr-gallery-images-contours-and-fields-specgram-demo-py)
10. [Matplotlib Colormaps 官方教學文件](https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html)
---
###### tags:`Python`
Sign in with Wallet
Connect another wallet