如何使用 python 获得声音包络
How to obtain sound envelope using python
您好，我是 python 的新手，还有声音信号分析。我想弄到一首诞生之歌（斑胸草雀）的信封。它的信号波动非常快，我尝试了不同的方法。例如，我尝试绘制信号并根据我发现的其他示例使用以下代码获取包络（我在代码上添加了注释以理解它）：
#Import the libraries
from pylab import *
import numpy
import scipy.signal.signaltools as sigtool
import scipy, pylab
from scipy.io import wavfile
import wave, struct
import scipy.signal as signal
#Open the txt file and read the wave file (also save it as txt file)
f_out = open('mike_1_44100_.txt', 'w')
w = scipy.io.wavfile.read("mike_1_44100_.wav") #here your sound file
a=w[1]
f_out.write('#time #z' + '\n')
#I print to check
print 'vector w'
print w[0],w[1]
print w
i=w[1].size
p=numpy.arange(i)*0.0000226 #to properly define the time signal with the sample rate
print 'vector p:'
print p
x=numpy.dstack([p,a])
print 'vector x:'
print x[0]
#saving file
numpy.savetxt('mike_1_44100_.txt',x[0])
f_out.close()
print 'i:'
print i
# num is the number of samples in the resampled signal.
num= np.ceil(float(i*0.0000226)/0.0015)
print num
y_resample, x_resample = scipy.signal.resample(numpy.abs(a),num, p,axis=0, window=('gaussian',150))
#y_resample, x_resample = scipy.signal.resample(numpy.abs(a), num, p,axis=-1, window=0)
#Aplaying a filter
W1=float(5000)/(float(44100)/2) #the frequency for the cut over the sample frequency
(b, a1) = signal.butter(4, W1, btype='lowpass')
aaa=a
slp =1* signal.filtfilt(b, a1, aaa)
#Taking the abs value of the signal the resample and finaly aplying the hilbert transform
y_resample2 =numpy.sqrt(numpy.abs(np.imag(sigtool.hilbert(slp, axis=-1)))**2+numpy.abs(np.real(sigtool.hilbert(slp, axis=-1)))**2)
print 'x sampled'
#print x_resample
print 'y sampled'
#print y_resample
xx=x_resample #[0]
yy=y_resample #[1]
#ploting with some style
plot(p,a,label='Time Signal') #to plot amplitud vs time
#plot(p,numpy.abs(a),label='Time signal')
plot(xx,yy,label='Resampled time signal Fourier technique Gauss window 1.5 ms ', linewidth=3)
#plot(ww,label='Window', linewidth=3)
#plot(p,y_resample2,label='Hilbert transformed sime signal', linewidth=3)
grid(True)
pylab.xlabel("time [s]")
pylab.ylabel("Amplitde")
legend()
show()
我在这里尝试了两件事，第一个是使用 scipy 的重采样函数来获取包络线，但是我对信号幅度有一些问题，我还不明白（我上传了图片用傅里叶技术获得但系统不允许我）：
第二种是用hilbert变换得到包络（现在我又上传了hilbert变换的图片系统不允许）可以运行我的代码得到这两个图片。但是我会用这个 link http://ceciliajarne.web.unq.edu.ar/?page_id=92&preview=true
现在 "envelope" 再次失败。我尝试像在某些示例中看到的那样过滤信号，但我的信号被衰减了，我无法获得包络。 任何人都可以帮助我的代码或更好的想法来获得信封吗？可以使用任何鸟鸣作为例子（我可以给你我的），但我需要看看复杂的声音而不是简单的信号会发生什么，因为它非常不同（简单的声音两种技术都可以）。
我还尝试调整我在以下位置找到的代码：http://nipy.org/nitime/examples/mtm_baseband_power.html
但是我无法为我的信号获取正确的参数，而且我不了解调制部分。我已经问过代码开发者了，等答案。
由于鸟鸣的“调制频率”可能远低于“载波频率”，即使振幅快速变化，也可以通过取信号的绝对值和然后应用长度为 20 毫秒的移动平均滤波器。
不过，您是否也对频率变化感兴趣，以充分表征歌曲？在这种情况下，对移动 window 进行傅里叶变换将为您提供更多信息，即作为时间函数的近似频率内容。这是我们人类听到的声音，可以帮助我们区分鸟类。
如果您不需要衰减，则不应应用巴特沃斯滤波器或采用移动平均，而应应用峰值检测。
移动平均值：每个输出样本是绝对值的平均值，例如。 50 个前面的输入样本。输出会衰减。
峰值检测：每个输出样本都是绝对值的最大值，例如50 个前面的输入样本。输出不会衰减。您可以在之后进行低通滤波以去除剩余的阶梯“纹波”。
您想知道为什么巴特沃斯滤波器会衰减您的信号。如果您的截止频率足够高，它几乎不会起作用，但它似乎会被强烈衰减。您的输入信号不是载波（哨声）和调制（包络）的总和，而是乘积。过滤将限制频率内容。剩下的是频率分量（项）而不是因子。您会看到衰减的调制（包络），因为该频率分量确实存在于您的信号中，比原始包络弱得多，因为它没有添加到您的载波中，而是与其相乘。由于与您的包络相乘的载波正弦波并不总是处于最大值，因此包络将被调制过程“衰减”，而不是您的滤波分析。
简而言之：如果您直接想要（乘法）包络而不是由于包络调制（乘法）而产生的（加法）频率分量，请采用峰值检测方法。
“Pythonish”伪代码中的峰值检测算法，仅供参考。
# Untested, but apart from typos this should work fine
# No attention paid to speed, just to clarify the algorithm
# Input signal and output signal are Python lists
# Listcomprehensions will be a bit faster
# Numpy will be a lot faster
def getEnvelope (inputSignal):
# Taking the absolute value
absoluteSignal = []
for sample in inputSignal:
absoluteSignal.append (abs (sample))
# Peak detection
intervalLength = 50 # Experiment with this number, it depends on your sample frequency and highest "whistle" frequency
outputSignal = []
for baseIndex in range (intervalLength, len (absoluteSignal)):
maximum = 0
for lookbackIndex in range (intervalLength):
maximum = max (absoluteSignal [baseIndex - lookbackIndex], maximum)
outputSignal.append (maximum)
return outputSignal
信号的包络可以使用相应analytic signal. Scipy implements the function scipy.signal.hilbert的绝对值来计算解析信号。
来自其文档：
我们创建一个频率从 20 Hz 增加到 100 Hz 的线性调频信号并应用调幅。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, chirp
duration = 1.0
fs = 400.0
samples = int(fs*duration)
t = np.arange(samples) / fs
signal = chirp(t, 20.0, t[-1], 100.0)
signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t))
幅度包络由分析信号的幅度给出。
analytic_signal = hilbert(signal)
amplitude_envelope = np.abs(analytic_signal)
看起来像
plt.plot(t, signal, label='signal')
plt.plot(t, amplitude_envelope, label='envelope')
plt.show()
它也可以用来计算瞬时频率（参见文档）。
您好，我是 python 的新手，还有声音信号分析。我想弄到一首诞生之歌（斑胸草雀）的信封。它的信号波动非常快，我尝试了不同的方法。例如，我尝试绘制信号并根据我发现的其他示例使用以下代码获取包络（我在代码上添加了注释以理解它）：
#Import the libraries
from pylab import *
import numpy
import scipy.signal.signaltools as sigtool
import scipy, pylab
from scipy.io import wavfile
import wave, struct
import scipy.signal as signal
#Open the txt file and read the wave file (also save it as txt file)
f_out = open('mike_1_44100_.txt', 'w')
w = scipy.io.wavfile.read("mike_1_44100_.wav") #here your sound file
a=w[1]
f_out.write('#time #z' + '\n')
#I print to check
print 'vector w'
print w[0],w[1]
print w
i=w[1].size
p=numpy.arange(i)*0.0000226 #to properly define the time signal with the sample rate
print 'vector p:'
print p
x=numpy.dstack([p,a])
print 'vector x:'
print x[0]
#saving file
numpy.savetxt('mike_1_44100_.txt',x[0])
f_out.close()
print 'i:'
print i
# num is the number of samples in the resampled signal.
num= np.ceil(float(i*0.0000226)/0.0015)
print num
y_resample, x_resample = scipy.signal.resample(numpy.abs(a),num, p,axis=0, window=('gaussian',150))
#y_resample, x_resample = scipy.signal.resample(numpy.abs(a), num, p,axis=-1, window=0)
#Aplaying a filter
W1=float(5000)/(float(44100)/2) #the frequency for the cut over the sample frequency
(b, a1) = signal.butter(4, W1, btype='lowpass')
aaa=a
slp =1* signal.filtfilt(b, a1, aaa)
#Taking the abs value of the signal the resample and finaly aplying the hilbert transform
y_resample2 =numpy.sqrt(numpy.abs(np.imag(sigtool.hilbert(slp, axis=-1)))**2+numpy.abs(np.real(sigtool.hilbert(slp, axis=-1)))**2)
print 'x sampled'
#print x_resample
print 'y sampled'
#print y_resample
xx=x_resample #[0]
yy=y_resample #[1]
#ploting with some style
plot(p,a,label='Time Signal') #to plot amplitud vs time
#plot(p,numpy.abs(a),label='Time signal')
plot(xx,yy,label='Resampled time signal Fourier technique Gauss window 1.5 ms ', linewidth=3)
#plot(ww,label='Window', linewidth=3)
#plot(p,y_resample2,label='Hilbert transformed sime signal', linewidth=3)
grid(True)
pylab.xlabel("time [s]")
pylab.ylabel("Amplitde")
legend()
show()
我在这里尝试了两件事，第一个是使用 scipy 的重采样函数来获取包络线，但是我对信号幅度有一些问题，我还不明白（我上传了图片用傅里叶技术获得但系统不允许我）：
第二种是用hilbert变换得到包络（现在我又上传了hilbert变换的图片系统不允许）可以运行我的代码得到这两个图片。但是我会用这个 link http://ceciliajarne.web.unq.edu.ar/?page_id=92&preview=true
现在 "envelope" 再次失败。我尝试像在某些示例中看到的那样过滤信号，但我的信号被衰减了，我无法获得包络。 任何人都可以帮助我的代码或更好的想法来获得信封吗？可以使用任何鸟鸣作为例子（我可以给你我的），但我需要看看复杂的声音而不是简单的信号会发生什么，因为它非常不同（简单的声音两种技术都可以）。
我还尝试调整我在以下位置找到的代码：http://nipy.org/nitime/examples/mtm_baseband_power.html
但是我无法为我的信号获取正确的参数，而且我不了解调制部分。我已经问过代码开发者了，等答案。
由于鸟鸣的“调制频率”可能远低于“载波频率”，即使振幅快速变化，也可以通过取信号的绝对值和然后应用长度为 20 毫秒的移动平均滤波器。
不过，您是否也对频率变化感兴趣，以充分表征歌曲？在这种情况下，对移动 window 进行傅里叶变换将为您提供更多信息，即作为时间函数的近似频率内容。这是我们人类听到的声音，可以帮助我们区分鸟类。
如果您不需要衰减，则不应应用巴特沃斯滤波器或采用移动平均，而应应用峰值检测。
移动平均值：每个输出样本是绝对值的平均值，例如。 50 个前面的输入样本。输出会衰减。
峰值检测：每个输出样本都是绝对值的最大值，例如50 个前面的输入样本。输出不会衰减。您可以在之后进行低通滤波以去除剩余的阶梯“纹波”。
您想知道为什么巴特沃斯滤波器会衰减您的信号。如果您的截止频率足够高，它几乎不会起作用，但它似乎会被强烈衰减。您的输入信号不是载波（哨声）和调制（包络）的总和，而是乘积。过滤将限制频率内容。剩下的是频率分量（项）而不是因子。您会看到衰减的调制（包络），因为该频率分量确实存在于您的信号中，比原始包络弱得多，因为它没有添加到您的载波中，而是与其相乘。由于与您的包络相乘的载波正弦波并不总是处于最大值，因此包络将被调制过程“衰减”，而不是您的滤波分析。
简而言之：如果您直接想要（乘法）包络而不是由于包络调制（乘法）而产生的（加法）频率分量，请采用峰值检测方法。
“Pythonish”伪代码中的峰值检测算法，仅供参考。
# Untested, but apart from typos this should work fine
# No attention paid to speed, just to clarify the algorithm
# Input signal and output signal are Python lists
# Listcomprehensions will be a bit faster
# Numpy will be a lot faster
def getEnvelope (inputSignal):
# Taking the absolute value
absoluteSignal = []
for sample in inputSignal:
absoluteSignal.append (abs (sample))
# Peak detection
intervalLength = 50 # Experiment with this number, it depends on your sample frequency and highest "whistle" frequency
outputSignal = []
for baseIndex in range (intervalLength, len (absoluteSignal)):
maximum = 0
for lookbackIndex in range (intervalLength):
maximum = max (absoluteSignal [baseIndex - lookbackIndex], maximum)
outputSignal.append (maximum)
return outputSignal
信号的包络可以使用相应analytic signal. Scipy implements the function scipy.signal.hilbert的绝对值来计算解析信号。
来自其文档：
我们创建一个频率从 20 Hz 增加到 100 Hz 的线性调频信号并应用调幅。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, chirp
duration = 1.0
fs = 400.0
samples = int(fs*duration)
t = np.arange(samples) / fs
signal = chirp(t, 20.0, t[-1], 100.0)
signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t))
幅度包络由分析信号的幅度给出。
analytic_signal = hilbert(signal)
amplitude_envelope = np.abs(analytic_signal)
看起来像
plt.plot(t, signal, label='signal')
plt.plot(t, amplitude_envelope, label='envelope')
plt.show()
它也可以用来计算瞬时频率（参见文档）。