[PYTHON] Linear Predictive Analysis (LPC) (= formant analysis)

Hello. "Linear Predictive Analysis (LPC) (20th audio signal processing in Python (2011/05/14): Finding LPC spectrum entrainment) I found an article called "(= formant analysis)", so I tried to run the code there almost as it is. This is an example of the vowel "A".

The code I ran is the same as it is there, but I'll write it (`` `levinson_durbin.py``` is also taken from there):

• The article "Performant analysis with Python (NumPy, SciPy)" But it seems to be using this.
• I also got the vowel audio file a.wav from somewhere.

# coding:utf-8
from __future__ import division
import wave
import numpy as np
import scipy.io.wavfile
import scipy.signal
import pylab as P
from levinson_durbin import autocorr, LevinsonDurbin

"""Find the LPC spectral envelope"""

def wavread(filename):
    wf = wave.open(filename, "r")
    fs = wf.getframerate()
    x = wf.readframes(wf.getnframes())
    x = np.frombuffer(x, dtype="int16") / 32768.0  # (-1, 1)Normalized to
    wf.close()
    return x, float(fs)

def preEmphasis(signal, p):
    """Emphasis filter"""
    #coefficient(1.0, -p)Create FIR filter for
    return scipy.signal.lfilter([1.0, -p], 1, signal)

def plot_signal(s, a, e, fs, lpcOrder, file):
    t = np.arange(0.0, len(s) / fs, 1/fs)
    #Find the positively predicted signal with LPC
    predicted = np.copy(s)
    #Since it is predicted from the past lpcOrder, the start index is from lpcOrder.
    #Prior to that, I was copying the original signal unpredictably
    for i in range(lpcOrder, len(predicted)):
        predicted[i] = 0.0
        for j in range(1, lpcOrder):
            predicted[i] -= a[j] * s[i - j]
    #Plot the original signal
    P.plot(t, s)
    #Plot positively predicted signals with LPC
    P.plot(t, predicted,"r",alpha=0.4)
    P.xlabel("Time (s)")
    P.xlim((-0.001, t[-1]+0.001))
    P.title(file)
    P.grid()
    P.show()
    return 0

def plot_spectrum(s, a, e, fs, file):
    #Find the amplitude spectrum of the LPC coefficient
    nfft = 2048   #Number of FFT samples
    fscale = np.fft.fftfreq(nfft, d = 1.0 / fs)[:nfft/2]
    #Logarithmic spectrum of the original signal
    spec = np.abs(np.fft.fft(s, nfft))
    logspec = 20 * np.log10(spec)
    P.plot(fscale, logspec[:nfft/2])
    #LPC logarithmic spectrum
    w, h = scipy.signal.freqz(np.sqrt(e), a, nfft, "whole")
    lpcspec = np.abs(h)
    loglpcspec = 20 * np.log10(lpcspec)
    P.plot(fscale, loglpcspec[:nfft/2], "r", linewidth=2)
    P.xlabel("Frequency (Hz)")
    P.xlim((-100, 8100))
    P.title(file)
    P.grid()
    P.show()
    return 0

def lpc_spectral_envelope(file):
    #Load audio
    wav, fs = wavread(file)
    t = np.arange(0.0, len(wav) / fs, 1/fs)
    #Cut out the central part of the audio waveform
    center = len(wav) / 2  #Center sample number
    cuttime = 0.04         #Length to cut out[s]
    s = wav[center - cuttime/2*fs : center + cuttime/2*fs]
    #Apply pre-emphasis filter
    p = 0.97         #Emphasis coefficient
    s = preEmphasis(s, p)
    #Humming window
    hammingWindow = np.hamming(len(s))
    s = s * hammingWindow
    #Find the LPC coefficient
#    lpcOrder = 12
    lpcOrder = 32
    r = autocorr(s, lpcOrder + 1)
    a, e  = LevinsonDurbin(r, lpcOrder)
    plot_signal(s, a, e, fs, lpcOrder, file)
    plot_spectrum(s, a, e, fs, file)
    return 0

if __name__ == "__main__":
    file = "a.wav"
    lpc_spectral_envelope(file)
    exit(0)