使用LPC在Python中估计共振峰

Question

我刚接触信号处理（还有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])

我用这个文件作为输入，我的脚本返回了这个列表：

[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)

这是根据Warren Weckesser的建议，现在我得到了：

[631.44354635609318, 1815.8629524985781, 3421.8288991389031, 6667.5030877036006]

我觉得我还是缺少了什么，因为F3的值偏差很大。

更新 2：

我意识到传给scikits.talkbox.lpc的order参数不对，是因为采样频率不同。我把它改成了：

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]

这和Praat的估计结果更接近了！

Answer 1

这里有至少两个问题：

• 根据链接的内容，“预强调滤波器是一个高通全极点（AR(1)）滤波器”。那里的系数符号是正确的：[1, 0.63]。如果你用[1, -0.63]，你就会得到一个低通滤波器。

• 你在使用scipy.signal.lfilter时，前两个参数的位置搞反了。

所以，试着把这个：

x1 = lfilter([1., -0.63], 1, x1)

改成这个：

x1 = lfilter([1.], [1., 0.63], x1)

我还没试过运行你的代码，所以不确定这是不是唯一的问题。

Answer 2

我没有得到你期待的结果，但我注意到两件可能导致差异的事情：

1. 你的代码使用了 [1, -0.63]，而你提供的链接中的MATLAB代码是 [1 0.63]。
2. 你的处理是一次性应用到整个 x 向量上，而不是分成更小的部分来处理（看看MATLAB代码是怎么做的：x = mtlb(I0:Iend);）。

---

希望这些能帮到你。

Answer 3

这个问题和传递给lpc函数的顺序有关。根据经验法则，2 + fs / 1000是一个公式，其中fs是采样频率。

你可以在这里找到更多信息：http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect10.html