線形予測

Last update: 2021/7/24

線形予測は、過去のデータで現在の値を予測する手法で、株価の予測などの時系列解析に用いられます。また、線形予測モデルの伝達関数は、全極モデルとなるため、共振系のモデルとしても用いられます。携帯やスマホなどで用いられている音声符号化は、この代表的な例です。

線形予測の概要

ここでは、線形予測を簡単に解説します。線形予測では、現在の値 x(n) の予測を、過去の M 個のデータ {x(n-1), ... x(n-M)} を用いて、式(1)のように行います。過去のデータの線形結合で現在の値 x(n) を予測するため、線形予測(Linear prediction, LP)と呼ばれます。推定誤差を式(2)のように定義し、この２乗平均を最小化すると、(3)の正規方程式が導かれます。これを解くと、線形予測係数 (Linear Prediction Coefficients, LPC) a=[a₁,...,a_M]が得られます。ここで、R及びrは式(4)で定義される相関行列及び相関ベクトル、その要素r(i)は、(5)で定義される相関関数です。これは、自己回帰モデル(autoregressive model)と呼ばれます。誤差v(n)をモデルへの入力と考えて、z-変換をとり、入出力間の伝達関数を求めると式(7)のようになります。ここで、a₀=1としています。誤差v(n)は白色（スペクトルが平坦）なので、伝達関数(7)を単位円上で評価した(8)は、{x(n)}の周波数スペクトルと考えることができます。S(ω)は、ARスペクトルと呼ばれます。

参考文献

• 1. S. Haykin, "Adaptie filter theory" 4th edition, chapter 3, Prentice hall, 2002

線形予測モジュール

linearPrediction.py

コンストラクタでは、線形予測係数a及びゲインg(=σ_v²)を推定します。関数spectrum()では、ARスペクトルS(ω)を計算します。ARスペクトルをコンストラクタではなく、関数で別途計算するようにしたのは、パラメタNFで周波数分解能を変更できるようにするためです。

#

# Linear prediction

#

# Ver 1.1, 2015/7/10 F.Asano

# Ver 1.2, 2018/6/29 F.Asano

import numpy as np

import scipy.linalg as alg

import scipy.signal as sig

#############################

# Class: LPC

# x: time series data

# NLP: Order of LP

#############################

class LPC:

def __init__(self,NLP):

self.NLP = NLP

def fit(self,x): # Training of LPC

NLP = self.NLP

NSample = len(x)

R = np.zeros((NLP,NLP),dtype='float')

r = np.zeros(NLP,dtype='float')

cor = np.correlate(x,x,mode='same')

Ncor = len(cor)

rr = cor[int(Ncor/2):int(Ncor/2)+NLP+1]

for m in range(NLP):

r[m] = rr[m+1]

for n in range(NLP):

R[m,n] = rr[np.abs(m-n)]

a = alg.solve(R,-r) # Normal equation

a = np.hstack((1.0,a)) # adding a_0 (=1.0)

v = sig.lfilter(a, np.ones(1), x) # Prediction error

g = np.var(v) # Variance of Prediction error (Gain for LP Spectrum)

self.a = a

self.g = g

def spectrum(self,NF): # LP Spectrum

NLP = self.NLP

p = np.zeros(int(NF/2)+1)

for n in range(int(NF/2)+1):

s = 1.0

omega = 2*np.pi*float(n)/float(NF)

for k in range(1,NLP+1):

s += self.a[k]*np.exp(-1j*omega*float(k))

p[n] = self.g/np.abs(s)

return p

サンプルプログラム

FrequencyAnalysis.py

このサンプルでは、時系列として、音声波形を読み込み、その一部 （図1(b)に示す５１２ポイントの波形、図１(a)の点線で挟まれた部分）の時系列に対して、FFT及びLPにより周波数スペクトルを求めます。図1(c)は、LPCを同定したあとの、予測誤差 v(n) を表しています。

図２及び３は、求めた周波数スペクトルです。両者を比較すると、FFTの方は、ギザギザしていますね。これは、音声の周期性（ピッチ）により、高調波が見えているためです。一方、LP法では、周期性の情報が、誤差v(n)の方に含まれるようになるため、ARスペクトル(8)には現れません。このため、声道などにおける共振に起因したスペクトルのピークだけが見えています。このように、LP法は、共振系のスペクトル分析に向いた方法と言えます。ちなみに、図２のピークは、フォルマントと呼ばれ、人の音声知覚やコンピュータによる自動音声認識で、重要な役目を果たします。

図４は、伝達関数(7)の極の位置をプロットしたものです。極の位置は、(7)の分母の多項式の根として得られます。これと図２を比較すると、極の位置とスペクトルのピークの位置が対応しています。また、極の位置は、いずれも単位円（複素平面上の半径1の円）の内側にあります。これは、ARモデル(6)の安定性と関係しています。極が単位円の外側にあると、モデルの出力が発散してしまいます。最後に、図５と６では、スペクトログラムを求めています。これは、図２と３の周波数分析を時間（フレーム）をずらしながら行うもので、周波数スペクトルの時間的な変化を見ることができます。

#

# Example of linear prediction

# Ver 1.2, 2018/6/39, F.Asano

# Library

import numpy as np

import matplotlib.pyplot as plt

import scipy.io.wavfile as wav

import scipy.fftpack as ft

import scipy.signal as sig

import linearPrediction as lp # Class for LP analysis

# データファイル

SourceFile = '../Data/speech11.wav'

# パラメータ

NLP = 16 # LP次数

FrameLength = 512 # 分析フレーム長

FrameShift = int(FrameLength/2) # フレームシフト

NFFT = FrameLength # FFT周波数分析の点数

NF = NFFT # LP周波数分析の点数

frameToBeAnalyzed = 32 # 分析対象のフレーム番号

DataLengthInSecond = 3 # [s]

# wavファイルの読み込み

Fs,data= wav.read(SourceFile) # Fs: サンプリング周波数

DataLength = Fs*DataLengthInSecond

z = np.float64(data[0:DataLength])

# 分析フレーム

fig = plt.figure()

tm = np.float64(range(DataLength))/Fs

plt.plot(tm,z)

plt.xlabel('Time [s]')

plt.title('Observation')

n = frameToBeAnalyzed

t1 = FrameShift*n

t2 = FrameShift*n + FrameLength

ymin,ymax = plt.ylim()

plt.vlines(t1/Fs,ymin,ymax,linestyles='dashed',colors='g')

plt.vlines(t2/Fs,ymin,ymax,linestyles='dashed',colors='g')

x = z[t1:t2]

fig = plt.figure()

plt.plot(x)

plt.title('Frame data')

# LP スペクトル

lps = lp.LPC(NLP)

lps.fit(x)

s = lps.spectrum(NF)

fig = plt.figure()

freq = Fs/NF*np.arange(int(NF/2)+1)

plt.plot(freq,20.0*np.log10(s))

plt.xlabel('Frequency [Hz]')

plt.ylabel('Spectrum [dB]')

plt.title('LP spectrum')

# 残差

print('a=',lps.a)

v = sig.lfilter(lps.a, np.ones(1), x)

plt.figure()

plt.plot(v)

plt.title('Prediction error')

# 極分析

pole = np.polynomial.polynomial.polyroots(lps.a) # finding roots

pole = 1.0/pole # for z-transform z^-1

theta = np.angle(pole)

r = np.abs(pole)

fig=plt.figure()

ax = plt.subplot(111,polar=True)

ax.plot(theta,r,'x')

ax.set_rmax(1.0)

ax.grid(True)

ax.set_title('Location of poles')

# FFT

w = np.hamming(NFFT)

s = ft.fft(z[t1:t2]*w,NFFT)

fig = plt.figure()

plt.plot(freq,20.0*np.log10(abs(s[0:int(NFFT/2)+1])))

plt.xlabel('Frequency [Hz]')

plt.ylabel('Spectrum [dB]')

plt.title('FFT spectrum')

# Spectrogram # スペクトログラム

NFrame = int(np.floor((DataLength-(FrameLength-FrameShift))/FrameShift ))

S = np.zeros((int(NFFT/2)+1,NFrame),dtype='float')

#tm = np.float64(range(NFrame))*float(FrameShift)/Fs

tm = np.arange(NFrame)*FrameShift/Fs

gx,gy = np.meshgrid(tm,freq)

for sMode in ['LP','FFT']:

for n in range(NFrame):

t1 = FrameShift*n

t2 = FrameShift*n + FrameLength

if sMode == 'LP':

lps.fit(z[t1:t2])

s = lps.spectrum(NF)

elif sMode == 'FFT':

s = np.abs( ft.fft(z[t1:t2]*w,NFFT) )[0:int(NFFT/2)+1]

S[:,n] = 20.0*np.log10(s)

fig = plt.figure()

im = plt.pcolor(gx,gy,S, cmap='coolwarm')

plt.xlabel('time [s]')

plt.ylabel('Frequency [Hz]')

plt.title(sMode + ' spectrogram')

plt.show()

図１：分析する時系列（音声）。赤の点線で挟まれている部分は周波数スペクトルを推定するフレーム

図４：線形予測器の極分析。×が極の位置。