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返回一个复数数组。幸运的是,就像警告说的那样,虚部将被丢弃。
怕什么真理无穷,进一寸有一寸的欢喜