Python中的音頻和數字信號處理（DSP）

翻譯自Python For Engineers。

1. 建立一個正弦波

在這個項目中，咱們將建立一個正弦波，並將其保存爲wav文件。

但在此以前，你應該知道一些理論。

頻率：頻率是正弦波重複一秒的次數。我將使用1KHz的頻率。

採樣率：大多數現實世界的信號是模擬的，而計算機是數字的。所以，咱們須要一個模數轉換器將模擬信號轉換爲數字信號。有關轉換器如何工做的詳細信息超出了本書的範圍。關鍵是採樣率，即轉換器每秒採樣模擬信號的次數。

如今，採樣率對咱們來講並不重要，由於咱們正在以數字方式完成全部工做，但咱們的正弦波公式須要它。我將使用48000的採樣率值，這是專業音頻設備中使用的值。
正弦波公式：公式以下。

y(t)= A * sin(2 * pi * f *t)

y(t) 是對於某個 x 軸 時間 t , y 軸的值

A 是幅值

pi 是咱們的老朋友 3.14159.

f 是頻率.

t 是樣本. 因爲咱們須要將其轉換爲數字，咱們將按採樣率進行劃分。

關於幅值：

上述的幅值A。在大多數書中，他們只選擇A的隨機值，一般爲1.但咱們這裏不能這麼取。咱們生成的正弦波將處於浮點狀態，雖然這對於繪製圖形來講已經足夠了，可是當咱們寫入文件時這將沒法工做。緣由是咱們須要處理整數。若是查看wave文件，它們將被寫爲16位短整數。若是咱們寫一個浮點數，它將沒法正確被表示。

爲了解決這個問題，咱們必須將浮點數轉換爲固定值。其中一種方法是將其乘以固定常數。咱們如何計算這個常數？帶符號的16位數的最大值是32767（2 ^ 15-1）。 （由於最左邊的位是爲符號保留的，留下15位。咱們將2增長到15的冪，而後減去1，由於計算機從0開始計數）。

因此咱們想要全音階音頻，咱們將它與32767相乘。但我想要的音頻信號是滿量程的一半，因此我將使用16000的幅度。

代碼以下：

import numpy as np
 
import wave
 
import struct
 
import matplotlib.pyplot as plt
 
# frequency is the number of times a wave repeats a second
 
frequency = 1000
 
num_samples = 48000
 
# The sampling rate of the analog to digital convert
 
sampling_rate = 48000.0
 
amplitude = 16000
 
file = "test.wav"

設置咱們剛纔聲明過的波形變量：

sine_wave = [np.sin(2 * np.pi * frequency * x/sampling_rate) for x in range(num_samples)]

它表示生成0到num_samples範圍內的 x，而且對於每一個x值，生成一個值爲該值的值。你能夠將此值視爲y軸值。而後將全部這些值放入列表中。十分簡單。

nframes=num_samples
 
comptype="NONE"
 
compname="not compressed"
 
nchannels=1
 
sampwidth=2

好的，如今是時候將正弦波寫入文件了。咱們將使用Python的內置wave庫。在這裏咱們設置參數。 nframes是幀數或樣本數。

comptype和compname都表示一樣的事情：數據未壓縮。 nchannels是通道數，即1. sampwidth是以字節爲單位的樣本寬度。正如我前面提到的，波形文件一般是每一個樣本16位或2個字節。

wav_file=wave.open(file, 'w')
 
wav_file.setparams((nchannels, sampwidth, int(sampling_rate), nframes, comptype, compname))

打開波形文件而且設置參數:

for s in sine_wave:
   wav_file.writeframes(struct.pack('h', int(s*amplitude)))

這裏可能須要解釋一下。咱們正在按樣本寫sine_wave文件。writeframes是寫正弦波的函數。一切都很簡單。可能讓你感到困惑的是下面這一行：

struct.pack('h', int(s*amplitude))

分解上述句子：

int(s*amplitude)

s 是咱們正在寫的sine_wave的單個樣本。我將它與此處的幅度相乘（轉換爲固定點）。放在這裏處理的緣由是寫入文件時須要轉化爲整數處理。

如今，咱們擁有的數據只是一個數字列表。若是咱們將其寫入文件，音頻播放器將沒法讀取。

Struct是一個Python庫，它接收咱們的數據並將其打包爲二進制數據。代碼中的h表示16位數。

要了解這個包的做用，讓咱們看看IPython中的一個例子。

In [1]: import numpy as np
 
In [2]: np.sin(0.5)
 
Out[2]: 0.47942553860420301
 
In [5]: 0.479*16000
 
Out[5]: 7664.0

上面是以0.5爲例。

所以咱們取0.5的sin，並將其乘以16000將其轉換爲固定點數。如今若是咱們將其寫入文件，它只會將7664寫爲字符串，這是錯誤的。讓咱們看一下struct作了什麼：

In [6]: struct.pack('h', 7664)
 
Out[6]: 'xf0x1d'

x 表示數字是十六進制。如你所見，struct已將咱們的數字7664轉換爲2個十六進制值：0xf0和0x1d。

爲何是0xf0 0x1d？好吧，若是你將7664轉換爲十六進制，你將得到0xf01d。
爲何兩個值？由於咱們使用的是16位值，並且咱們的數字不能合二爲一。所以，struct將其分爲兩個數字。

因爲數字如今是十六進制，所以其餘程序能夠讀取它們，包括咱們的音頻播放器。

回到咱們的代碼：

for s in sine_wave:
   wav_file.writeframes(struct.pack('h', int(s*amplitude)))

這將採用咱們的正弦波樣本並將其寫入咱們的文件test.wav，打包爲16位音頻。

在你擁有的任何音頻播放器中播放文件 - Windows Media Player，VLC等。你應該能聽到很是短的蜂鳴。

如今，咱們須要檢查音調的頻率是否正確。我將使用Audacity，一個具備大量功能的開源音頻播放器。其中之一就是咱們能夠找到音頻文件的頻率。讓咱們打開Audacity。

咱們獲得了一個正弦波。請注意，波形高達0.5，而1.0是最大值。還記得咱們乘以16000，這是36767的一半。

如今找頻率。轉到編輯 - >全選（或按Ctrl A），而後選擇分析 - >頻譜分析。

能夠看到峯值大約爲1000 Hz，這就是咱們建立波形文件的方式。

2. 計算正弦波的頻率

我在個人學位課程中參加了一門信號處理課程，而且不瞭解任何事情。咱們被要求解一百個方程式，沒有任何意義或邏輯。我發現這個主題很無聊和迂腐。

當我不得再也不爲個人碩士學習時，我不高興。但我很幸運。

此次，老師是一名執業工程師。他經營本身的公司併兼職教學。與大學教師不一樣，他實際上知道方程式的用途。

可是這位老師（我忘了他的名字，他是一個丹麥人）向咱們展現了一個嘈雜的信號，而且使用了DFT。而後他在圖形窗口中顯示結果。咱們清楚地看到了原始的正弦波和噪聲頻率，我第一次理解了DFT的做用。

使用 DFT 獲取頻率

要得到正弦波的頻率，你須要得到其離散傅立葉變換（DFT）。你不須要了解如何推導變換。你只須要知道如何使用它。

最簡單的說，DFT接收信號並計算其中存在哪些頻率。在更多技術術語中，DFT將時域信號轉換爲頻域。那是什麼意思？讓咱們來看看咱們的正弦波。

波形隨着時間而變化。若是這是一個音頻文件，你能夠想象播放器在文件播放時向右移動。

在頻域中，你能夠看到信號的頻率部分。此圖片將從本章後面的內容中獲取，以顯示頻域的外觀：

若是添加或刪除頻率，信號將發生變化，但不會及時更改。例如，若是你採用1000 Hz的音頻並採用其頻率，不管你看多長時間，頻率都將保持不變。可是若是你在時域中看它，你會看到信號在移動。

DFT在計算機上運行得很慢（早在70年代），所以發明了快速傅里葉變換（FFT）。 FFT現在被普遍使用。

它的工做方式是，接收信號並在其上運行FFT，而後你會獲得信號的頻率。

若是你從未使用過（甚至沒有聽過）FFT，請不要擔憂。我將教你如何開始使用它，若是你願意，你能夠在線閱讀更多內容。不管如何，大多數教程或書籍都不會教你太多。他們一般會用公式來轟炸你，而不會告訴你如何處理它們。

frame_rate = 48000.0
 
infile = "test.wav"
 
num_samples = 48000
 
wav_file = wave.open(infile, 'r')
 
data = wav_file.readframes(num_samples)
 
wav_file.close()

咱們正在讀取咱們在上一個例子中生成的wave文件。這段代碼應該足夠清楚。wave.readframes()函數從wave文件中讀取全部音頻幀。

data = struct.unpack('{n}h'.format(n=num_samples), data)

還記得咱們必須打包數據以使其以二進制格式可讀嗎？好吧，咱們如今正好相反。該函數的第一個參數是格式字符。告訴解包器按照num_samples16位字來解壓縮文件（記住h表示16位）。

data = np.array(data)

而後咱們將數據轉換爲numpy數組。

data_fft = np.fft.fft(data)

咱們接受數據的fft。這將建立一個包含信號中全部頻率的陣列。

如今，這裏就是問題所在。 fft返回一個複數的數組，不告訴咱們任何東西。若是我打印出fft的前8個值，會獲得：

In [3]: data_fft[:8]
 
Out[3]:
 
array([ 13.00000000 +0.j        ,   8.44107682 -4.55121351j,
 
6.24696630-11.98027552j,   4.09513760 -2.63009999j,
 
-0.87934285 +9.52378503j,   2.62608334 +3.58733642j,
 
4.89671762 -3.36196984j,  -1.26176048 +3.0234555j ])

Numpy 能夠將複數轉換爲實數值。

# This will give us the frequency we want
 
frequencies = np.abs(data_fft)

numpy abs()函數將獲取咱們的複數信號並生成它的實數值。

旁註

稍微解釋一下FFT如何返回其結果。

FFT返回信號中全部可能的頻率。它返回的方式是每一個索引包含一個頻率元素。假設你將FFT結果存儲在名爲data_fft的數組中。而後：

data_fft [1]將包含1 Hz的頻率部分。

data_fft [2]將包含2 Hz的頻率部分。

...

data_fft [8]將包含8 Hz的頻率部分。

...

data_fft [1000]將包含1000 Hz的頻率部分。

如今若是你的信號中沒有1Hz頻率怎麼辦？你仍然能夠在data_fft [1]得到一個值，但它將是微不足道的。舉個例子，咱們來看看1000 Hz波的實際fft：

data_fft = np.fft.fft(sine_wave)
 
abs(data_fft[0])
 
Out[7]: 8.1289678326462086e-13
 
abs(data_fft[1])
 
Out[8]: 9.9475299243014428e-12
 
abs(data_fft[1000])
 
Out[11]: 24000.0

若是你查看data_fft [0]或data_fft [1]的絕對值，你會發現它們很小。最後的e-12意味着它們被提高到-12的冪，因此對於data_fft [0]來講就像0.00000000000812。可是若是你看一下data_fft [1000]，那麼這個值就是24000.這能夠很容易地繪製出來。

若是咱們想要找到具備最高值的數組元素，咱們能夠經過如下方式找到它：

print("The frequency is {} Hz".format(np.argmax(frequencies)))

np.argmax將返回信號中的最高頻率，而後打印出來。正如咱們手動看到的，這是1000Hz（或存儲在data_fft [1000]的值）。如今咱們也能夠繪製數據。

plt.subplot(2,1,1)
 
plt.plot(data[:300])
 
plt.title("Original audio wave")
 
plt.subplot(2,1,2)
 
plt.plot(frequencies)
 
plt.title("Frequencies found")
 
plt.xlim(0,1200)
 
plt.show()

subplot（2,1,1）意味着咱們正在繪製一個2×1網格。第三個數字是圖號，是惟一一個會改變的號碼。

就是這樣。咱們獲取了音頻文件並計算了它的頻率。接下來，咱們將爲咱們的情景添加噪音，而後嘗試清理它。

3. 簡單分離噪聲

frequency  = 1000
noisy_freq = 50
num_samples = 48000

sampling_rate = 48000.0

非噪聲的頻率是1000Hz，咱們加入50Hz的噪聲。

# Create the sine wave and noise

sine_wave = [np.sin(2*np.pi*frequency*x1/sampling_rate) for x1 in range(num_samples)]

sine_noise = [np.sin(2*np.pi*noisy_freq*x1/sampling_rate) for x1 in range(num_samples)]

# Convert them to numpy arrays

sine_wave = np.array(sine_wave)
sine_noise = np.array(sine_noise)

上面的構建代碼和以前的相似。咱們生成了兩個正弦波，其中是一個信號另外一個是噪聲，而且咱們將他們都轉化爲了一個numpy的數組

# Add them to create a noisy signal

combined_signal = sine_wave+sine_noise

我把噪聲加到了信號裏。正如以前提到的，numpy纔有這麼方便的累加方式。若是是普通的python列表，則須要寫一個for循環。總之numpy很方便啦。

把到目前爲止獲得的信號圖像化：

plt.subplot(3,1,1)

plt.title('Original sine wave')

# Need to add empty space, else everything looks scrunched up!

plt.subplots_adjust(hspace=.5)

plt.plot(sine_wave[:500])

plt.subplot(3,1,2)

plt.title('Noise wave')

plt.plot(sine_noise[:4000])

plt.subplot(3,1,3)

plt.title('Original + Noise')

plt.plot(combined_signal[:3000])

plt.show()

data_fft = np.fft.fft(combined_signal)

freq = (np.abs(data_fft[:len(data_fft)]))

data_fft 包含了噪聲+信號波的 fft. freq 包含了其中發現的頻率的絕對值。

plt.plot(freq)

plt.title('Before filtering: Will have main signal(1000Hz) + noise frequency(50Hz)')

plt.xlim(0,1200)

如今來過濾信號。我不會詳細介紹過濾的每一個細節(由於那樣可能要講一整本書)。我將使用 fft 來建立一個簡單的過濾器。

下面是完整的代碼：

filtered_freq = []
index = 0
for f in freq:
    # Filter between lower and upper limits
    # Chooing 950, as closest to 1000. In real world, won't get exact numbers like these
    if index > 950 and index<1050:
        # Has a real value. I'm choosing >1 , as many values are like 0.00000001 etc
        if f>1:
            filtered_freq.append(f)
        else:
            filtered_freq.append(0)
    else:
        filtered_freq.append(0)
    index+=1

上述代碼建立一個空列表 filtered_freq。若是你還記得的話， freq 存儲了 fft 的絕對值。

index 是數組 freq 當前的索引。正如我說的，fft 返回了信號中全部的頻率。
它們基於索引存儲在數組中，所以 freq[1] 的頻率爲1Hz， freq[2] 的頻率爲2Hz，依此類推。

由於我知道個人目標信號頻率是1000Hz，因此我會在這個值附近搜索它。在現實世界中，咱們永遠不會獲得確切的頻率，由於噪聲意味着一些數據將會丟失。這裏使用950的下限和1050的上限，代碼檢查咱們循環的頻率是否在這個範圍內。

至於嵌套的 if else，一樣在前文中有所說起：雖然全部頻率都會有值來表示，但它們的絕對值將是微小的，一般小於1.所以，若是咱們找到大於1的值，咱們將其保存到filtered_freq數組中。

若是咱們的頻率不在咱們正在尋找的範圍內，或者若是該值過低，直接0。這是爲了刪除咱們不想要的全部頻率。而後咱們移動索引到下一這個值。

接下來畫出咱們濾波後的圖像：

plt.plot(filtered_freq)

plt.title('After filtering: Main signal only(1000Hz)')

plt.xlim(0,1200)

plt.show()

plt.close()

recovered_signal = np.fft.ifft(filtered_freq)

如今咱們使用 ifft，即 FFT 的逆過程。這將接收咱們的信號並將其轉換回時域。咱們如今能夠將它與原始噪聲信號進行比較。

plt.subplot(3,1,1)
 
plt.title("Original sine wave")
 
# Need to add empty space, else everything looks scrunched up!
 
plt.subplots_adjust(hspace=.5)
 
plt.plot(sine_wave[:500])
 
plt.subplot(3,1,2)
 
plt.title("Noisy wave")
 
plt.plot(combined_signal[:4000])
 
plt.subplot(3,1,3)
 
plt.title("Sine wave after clean up")
 
plt.plot((recovered_signal[:500]))
 
plt.show()

注意咱們會看到一個警告：

ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)

這是由於fft返回一個複數數組。幸運的是，就像警告說的那樣，虛部將被丟棄。