在Python中使用LPC估算共振峰

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 的估计值!

- pcaisse

---

你能计算/显示UCLA信号的声谱图吗？ 在MATLAB示例中有两个元音（发声"MATLAB"），并且您可以在时频图中清楚地看到它们。 在MATLAB示例中分析介于0.1和0.25秒之间。 不知道UCLA信号的内容，很难建议要分析哪个片段。但我猜想你想先分析包含一个元音的子集。 - paisanco

Answer 1

问题与传递给lpc函数的顺序有关。根据http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect10.html，2 + fs / 1000（其中fs是采样频率）是一个经验法则。

- pcaisse

Answer 2

至少有两个问题：

• 根据链接，"预加重滤波器是一个高通全极（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)

我还没有尝试运行你的代码，所以我不知道那些是否是唯一的问题。

- Warren Weckesser

---

我将那一行代码更改为x1 = lfilter([1], [1., 0.63], x1)（第一个参数必须是类似数组的对象）。这给我带来了以下结果：[631.44354635609318, 1815.8629524985781, 3421.8288991389031, 6667.5030877036006]。但我仍然感觉缺少某些东西，因为第三个共振峰偏差很大。 - pcaisse

我修复了缺失的括号 - 谢谢。你确定 scikits.talkbox 中的 lpc 与 matlab 的 lpc 给出相同的结果吗？你有 matlab 可以比较一下结果吗？ - Warren Weckesser

Answer 3

我无法得到您期望的结果，但我注意到两件事可能会导致一些差异：

1. 您的代码使用 [1, -0.63]，而您提供的 MATLAB 代码使用 [1 0.63]。
2. 您的处理是一次性应用于整个 x 向量，而不是它的较小部分（请参见 MATLAB 代码中的此处：x = mtlb(I0:Iend);）。

---

希望这有所帮助。

- Lukeclh

---

1. 我从这里得到了预加重滤波器代码： http://pydoc.net/Python/scikits.talkbox/0.2.5/scikits.talkbox.features.mfcc/

我尝试过 [1., 0.63] 和 [1., -0.63]，但都没有给我期望的结果。

2. 我尝试将 x 设置为信号的中间三分之一，但在返回的频率上几乎没有什么区别。

- pcaisse