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_samples 16位字來解壓縮文件（記住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返回一個復數數組。幸運的是，就像警告說的那樣，虛部將被丟棄。