Python中的⾳频和数字信号处理（DSP）
翻译⾃。
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)]python新手能做啥兼职
它表⽰⽣成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 = adframes(num_samples)
wav_file.close()
我们正在读取我们在上⼀个例⼦中⽣成的wave⽂件。这段代码应该⾜够清楚。adframes()函数从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

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。