我对信号处理还不太熟悉(关于这一点,numpy、scipy和matlab)。我试着用Python中的LPC来估计元音共振峰,方法是修改下面的matlab代码:
http://www.mathworks.com/help/signal/ug/formant-estimation-with-lpc-coefficients.html
这是我目前的代码:
#!/usr/bin/env python
import sys
import numpy
import wave
import math
from scipy.signal import lfilter, hamming
from scikits.talkbox import lpc
"""
Estimate formants using LPC.
"""
def get_formants(file_path):
# Read from file.
spf = wave.open(file_path, 'r') # http://www.linguistics.ucla.edu/people/hayes/103/Charts/VChart/ae.wav
# Get file as numpy array.
x = spf.readframes(-1)
x = numpy.fromstring(x, 'Int16')
# Get Hamming window.
N = len(x)
w = numpy.hamming(N)
# Apply window and high pass filter.
x1 = x * w
x1 = lfilter([1., -0.63], 1, x1)
# Get LPC.
A, e, k = lpc(x1, 8)
# Get roots.
rts = numpy.roots(A)
rts = [r for r in rts if numpy.imag(r) >= 0]
# Get angles.
angz = numpy.arctan2(numpy.imag(rts), numpy.real(rts))
# Get frequencies.
Fs = spf.getframerate()
frqs = sorted(angz * (Fs / (2 * math.pi)))
return frqs
print get_formants(sys.argv[1])
使用this file作为输入,我的脚本返回以下列表:
[682.18960189917243, 1886.3054773107765, 3518.8326108511073, 6524.8112723782951]
我甚至还没到最后一步,他们根据带宽过滤频率,因为列表中的频率不对。根据Praat的说法,我应该得到这样的结果(这是元音中间的共振峰列表):
Time_s F1_Hz F2_Hz F3_Hz F4_Hz
0.164969 731.914588 1737.980346 2115.510104 3191.775838
我做错什么了?
非常感谢
我改了这个
x1 = lfilter([1., -0.63], 1, x1)
到
x1 = lfilter([1], [1., 0.63], x1)
按照沃伦·韦克瑟的建议,现在
[631.44354635609318, 1815.8629524985781, 3421.8288991389031, 6667.5030877036006]
我觉得我好像丢了什么东西,因为F3非常不好。
我意识到传递给order
的scikits.talkbox.lpc
由于采样频率的不同而关闭。更改为:
Fs = spf.getframerate()
ncoeff = 2 + Fs / 1000
A, e, k = lpc(x1, ncoeff)
现在我明白了:
[257.86573127888488, 774.59006835496086, 1769.4624576002402, 2386.7093679399809, 3282.387975973973, 4413.0428174593926, 6060.8150432549655, 6503.3090645887842, 7266.5069407315023]
更接近普拉特的估计!
至少有两个问题:
根据该链路,“预加重滤波器是高通全极点(AR(1))滤波器”。给出的系数符号是正确的:
[1, 0.63]
。如果你使用[1, -0.63]
,你会得到一个低通滤波器。前两个参数被^{} 反转。
所以,试着改变这个:
对此:
我还没试过运行你的代码,所以我不知道这些是不是唯一的问题。
这个问题与传递给lpc函数的顺序有关。
2 + fs / 1000
其中fs
是采样频率,根据经验法则是:http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect10.html
我没能得到你所期望的结果,但我确实注意到两件事可能会导致一些差异:
[1, -0.63]
,其中您提供的链接中的MATLAB代码具有[1 0.63]
。x
向量,而不是它的较小片段(请参见MATLAB代码的作用:x = mtlb(I0:Iend);
)。希望能有所帮助。
相关问题 更多 >
编程相关推荐