matlab lpc求共振峰频率,在Python中使用LPC估计共振峰

我是信号处理的新手(而且是numpy,scipy和matlab).我正在尝试通过调整这个matlab代码来估计

Python中使用LPC的元音共振峰：

到目前为止,这是我的代码：

#!/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)

按照Warren Weckesser的建议,我现在正在接受

[631.44354635609318,1815.8629524985781,3421.8288991389031,6667.5030877036006]

我觉得我错过了一些东西,因为F3非常关闭.

更新2：

我意识到由于采样频率的不同,传递给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]

更接近Praat的估计！