線形予測
Last update: 2021/7/24
線形予測は、過去のデータで現在の値を予測する手法で、株価の予測などの時系列解析に用いられます。また、線形予測モデルの伝達関数は、全極モデルとなるため、共振系のモデルとしても用いられます。携帯やスマホなどで用いられている音声符号化は、この代表的な例です。
線形予測の概要
ここでは、線形予測を簡単に解説します。線形予測では、現在の値 x(n) の予測を、過去の M 個のデータ {x(n-1), ... x(n-M)} を用いて、式(1)のように行います。過去のデータの線形結合で現在の値 x(n) を予測するため、線形予測(Linear prediction, LP)と呼ばれます。推定誤差を式(2)のように定義し、この2乗平均を最小化すると、(3)の正規方程式が導かれます。これを解くと、線形予測係数 (Linear Prediction Coefficients, LPC) a=[a1,...,aM]が得られます。ここで、R及びrは式(4)で定義される相関行列及び相関ベクトル、その要素r(i)は、(5)で定義される相関関数です。これは、自己回帰モデル(autoregressive model)と呼ばれます。誤差v(n)をモデルへの入力と考えて、z-変換をとり、入出力間の伝達関数を求めると式(7)のようになります。ここで、a0=1としています。誤差v(n)は白色(スペクトルが平坦)なので、伝達関数(7)を単位円上で評価した(8)は、{x(n)}の周波数スペクトルと考えることができます。S(ω)は、ARスペクトルと呼ばれます。
参考文献
S. Haykin, "Adaptie filter theory" 4th edition, chapter 3, Prentice hall, 2002
線形予測モジュール
linearPrediction.py
コンストラクタでは、線形予測係数a及びゲインg(=σv2)を推定します。関数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)に示す512ポイントの波形、図1(a)の点線で挟まれた部分)の時系列に対して、FFT及びLPにより周波数スペクトルを求めます。図1(c)は、LPCを同定したあとの、予測誤差 v(n) を表しています。
図2及び3は、求めた周波数スペクトルです。両者を比較すると、FFTの方は、ギザギザしていますね。これは、音声の周期性(ピッチ)により、高調波が見えているためです。一方、LP法では、周期性の情報が、誤差v(n)の方に含まれるようになるため、ARスペクトル(8)には現れません。このため、声道などにおける共振に起因したスペクトルのピークだけが見えています。このように、LP法は、共振系のスペクトル分析に向いた方法と言えます。ちなみに、図2のピークは、フォルマントと呼ばれ、人の音声知覚やコンピュータによる自動音声認識で、重要な役目を果たします。
図4は、伝達関数(7)の極の位置をプロットしたものです。極の位置は、(7)の分母の多項式の根として得られます。これと図2を比較すると、極の位置とスペクトルのピークの位置が対応しています。また、極の位置は、いずれも単位円(複素平面上の半径1の円)の内側にあります。これは、ARモデル(6)の安定性と関係しています。極が単位円の外側にあると、モデルの出力が発散してしまいます。最後に、図5と6では、スペクトログラムを求めています。これは、図2と3の周波数分析を時間(フレーム)をずらしながら行うもので、周波数スペクトルの時間的な変化を見ることができます。
#
# 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()
(a)全体の波形
(b)分析フレームの波形
(c) 予測誤差
図1:分析する時系列(音声)。赤の点線で挟まれている部分は周波数スペクトルを推定するフレーム
図2:線形予測による周波数スペクトル
図3:FFTによる周波数スペクトル
図4:線形予測器の極分析。×が極の位置。
図5:線形予測によるスペクトログラム
図6:FFTによるスペクトログラム