誤差関数and最急降下法を用いて音声の周波数のスペクトルをフィッティングで予測＆行列計算でフィッティングを効率化
前回はこちら。
引き続き音声の周波数のスペクトルを調べます。FFT(高速フーリエ変換)を使わずに音声ファイルから元の声の周波数の分布を調べることに拘っていこうと思います。
今回はこちらのページを参考にしました。
こちらのページでは、離散データ(波形データなど)を$${=a cos(x_i+b)+c}$$の形にして、コサイン関数に当てはめて元のデータの波形を類推しています。類推というのは、xを時間、yを振幅データとするとき、できるだけ元のデータとの距離が近づくように変数a,b,cを動かしていくことで、一番元のデータに近い別のx,yに関する関数を作りあげて、当てはめます。
この記事で元の関数を求めるプロセスを応用して、$${x_i}$$の係数を$${2π×f}$$(fは求める波形型関数の周波数)複数設定して、それらの総和を取ることで、さまざまな周波数を合成した波と仮定した時の、波形データを推測する関数を求めることができます。
前回の記事で作った音波の周波数のスペクトルを求めるプログラムと同じものを求めるためのプロセスを、違う角度からアプローチします。また、同時に、求めたスペクトルから音声ファイルを作るためのプログラムも組みました。後でフィッティングした関数から作った音声データを上げるのですが、実際には人間の発声とはほど遠い音声ができることになります。いろいろな要因があって、実際の波形に近いデータを求めることはできなかったのであろうと思います。しかし、音声の周波数を考えるためにいい素材がたくさんできたので、記録のために書き残しておこうと思います。
コードはこちらです。
import numpy as np
import wave
import struct
import csv
#パラメータ
#--------------------------------------------------------------------------------------------------------------------
volume=0.5#max->1
myu=0.00030#study rate
A=4.0e-6#output shrink
reverse=1#-1->progress reverse slope
repeat=0#repeat set
part=1100
fs = 44100
f0 = 1050
sec = 2 #秒
#--------------------------------------------------------------------------------------------------------------------
def create_wave(A,f0,fs,t):
heltz=[1100,1000,900,800,706,654,628,592,560,528,498,471,444,419,396,373,353,333,314,296,280,264,249,235,222,209,198,186,176.5,166,157,148,140,132,124,117,111,104,99,93,88.25,83,78.62,74.20,70.04,66.11,62.40,58.89,55.59,52.47,49.52,46.74,44.125,41.648,39.310,37.104,35.022,33.056,31.201,29.449,27.797,26.236,24.764,23.374,22.0625,20.824,19.655,18.552,17.511,16.528,15.600,14.724,13.898,13.118,12.382,11.687,11.03125,10.412,9.8277,9.2761,8.7555,8.2640,7.8002,7.3624,6.9492,6.5592,6.1910,5.8436,5.515625,5.2060,4.9138,4.6380,4.3777,4.1320,3.9001,3.6812,3.4746,3.2796,3.0955,2.9218,2.7578,2.5,2.3,2.1,1.9,1.7,1.5,1.3]#96
arrsin=[4.648648150089610,7.125893363745060,8.12128929144061,8.995253173061090,9.019686893945860,8.013992837798230,5.828934870532150,2.618399670427710,-1.8577641395080900,-7.001092508602710,-12.190048402042500,-17.06966125591930,-20.13872741675840,-21.01265690282240,-18.806120309797500,-13.488244182197700,-6.5484309181973500,0.8215731842666290,5.325985475256930,4.279713626403720,-3.261143352564350,-16.419255529029000,-7.96542823648640,-3.217201672707100,-6.86869968606550,-10.209905427906100,1.7524334883575500,-3.8430858377678600,-27.868658210467400,-1.676403346626200,-41.60480013863580,-5.66919429368504,0.41467946639648400,-4.01141141814760,-8.75111181462200,-18.44174160224190,4.260547320863400,-270.6468542574310,-9.5310379467390,-0.22522412538200,8.64234755414400,0.29209236074570,-4.1247171719036,7.3770978240737,-4.06868664106200,0.39174151878330,-4.99612583759700,4.142885693936300,-2.61697324162500,20.207060462070500,-2.9181576379800,907.145709977173,19.470109890193500,72.07795849386550,3.22122381156500,111.2348714529820,2.2440159986023,4.753062154594100,0.15326190156680,3.111509360685500,8.21491838945830,1.36905543556010,16.71791383119230,9.003846712536600,32.90936744659320,9.605809780918080,28.706476283615800,21.28875323166190,13.53411924510910,11.13394553957920,9.344813019712000,9.701360191948100,20.959425256237200,10.128894551045500,9.220292181701900,8.715136777400750,19.741149898801600,-3.816222611862580,1.9954380740322500,7.09980203958077,23.453057947548900,6.129996545803530,5.241976452131390,4.741589855078980,3.4927810790820800,2.71312342772176,2.0772169840177300,2.546182086754480,2.950857834497970,0.32592955288792800,1.5029298799370100,1.1887813049544200,0.5818171944845080,3.656763129258720,2.92469631728601,0.4682249973694840,0,0,0,0,0,0,0,0,0,0,0]
arrcos=[4.648648150089610,7.125893363745060,8.12128929144061,8.995253173061090,9.019686893945860,8.013992837798230,5.828934870532150,2.618399670427710,-1.8577641395080900,-7.001092508602710,-12.190048402042500,-17.06966125591930,-20.13872741675840,-21.01265690282240,-18.806120309797500,-13.488244182197700,-6.5484309181973500,0.8215731842666290,5.325985475256930,4.279713626403720,-3.261143352564350,-16.419255529029000,-7.96542823648640,-3.217201672707100,-6.86869968606550,-10.209905427906100,1.7524334883575500,-3.8430858377678600,-27.868658210467400,-1.676403346626200,-41.60480013863580,-5.66919429368504,0.41467946639648400,-4.01141141814760,-8.75111181462200,-18.44174160224190,4.260547320863400,-270.6468542574310,-9.5310379467390,-0.22522412538200,8.64234755414400,0.29209236074570,-4.1247171719036,7.3770978240737,-4.06868664106200,0.39174151878330,-4.99612583759700,4.142885693936300,-2.61697324162500,20.207060462070500,-2.9181576379800,907.145709977173,19.470109890193500,72.07795849386550,3.22122381156500,111.2348714529820,2.2440159986023,4.753062154594100,0.15326190156680,3.111509360685500,8.21491838945830,1.36905543556010,16.71791383119230,9.003846712536600,32.90936744659320,9.605809780918080,28.706476283615800,21.28875323166190,13.53411924510910,11.13394553957920,9.344813019712000,9.701360191948100,20.959425256237200,10.128894551045500,9.220292181701900,8.715136777400750,19.741149898801600,-3.816222611862580,1.9954380740322500,7.09980203958077,23.453057947548900,6.129996545803530,5.241976452131390,4.741589855078980,3.4927810790820800,2.71312342772176,2.0772169840177300,2.546182086754480,2.950857834497970,0.32592955288792800,1.5029298799370100,1.1887813049544200,0.5818171944845080,3.656763129258720,2.92469631728601,0.4682249973694840,0,0,0,0,0,0,0,0,0,0,0]
a_dush=[0]*107
b_dush=[0]*107
#print("Hello")
with open("output2.csv") as f:
reader=list(csv.reader(f))
#print(reader)
a_dush=[float(i)for i in reader[0]]
b_dush=[float(i)for i in reader[1]]
print(len(a_dush),len(b_dush))
for i in range(107):
float(a_dush[i])
float(b_dush[i])
#--------------------------------------------------------------------------------------------------------------------
rad=[0]*107
theta=[0]*107
#for i in range(107):
rad=np.sqrt(np.power(a_dush,2)+np.power(b_dush,2))
theta=np.arctan2(b_dush,a_dush)
for i in range(107,107):
if rad[i]<5.0e+18/np.power(2,i/12):
rad[i]=0
for i in range(107,107):
rad[i]=0
for i in range(0):
rad[i]=0
#print(rad)
point = np.arange(0,fs*t)
sin_wave=[0]*fs*t
# print(sin_wave,rad,theta)
#for i in range(fs*t):
for n in range(107):
sin_wave =sin_wave +A*rad[n]*np.cos(2*np.pi*point/heltz[n]+theta[n])
max=0
for i in range(fs*t):
if abs(sin_wave[i])>max:
max=abs(sin_wave[i])
#print(sin_wave)
sin_wave = [int(x * 32767.0/max/2) for x in sin_wave]
print(sin_wave[:10])
binwave = struct.pack("h" * len(sin_wave), *sin_wave)
#print("Hello")
w = wave.Wave_write("output.wav")
p = (1, 2, fs, len(binwave), 'NONE', 'not compressed')
w.setparams(p)
w.writeframes(binwave)
w.close()
#create_wave(A, f0, fs, sec)
#print("Hello ")
from scipy.io import wavfile
import math
import array
import csv
def jijou():
wavfil = 'output.wav'
fs , data = wavfile.read(wavfil,mmap=False)
print (fs)
length=len(data)/1102
length=math.floor(length)-1
print (length)
print( 'HelloWorld')
print(data[:10])
heltz=[1100,1000,900,800,706,654,628,592,560,528,498,471,444,419,396,373,353,333,314,296,280,264,249,235,222,209,198,186,176.5,166,157,148,140,132,124,117,111,104,99,93,88.25,83,78.62,74.20,70.04,66.11,62.40,58.89,55.59,52.47,49.52,46.74,44.125,41.648,39.310,37.104,35.022,33.056,31.201,29.449,27.797,26.236,24.764,23.374,22.0625,20.824,19.655,18.552,17.511,16.528,15.600,14.724,13.898,13.118,12.382,11.687,11.03125,10.412,9.8277,9.2761,8.7555,8.2640,7.8002,7.3624,6.9492,6.5592,6.1910,5.8436,5.515625,5.2060,4.9138,4.6380,4.3777,4.1320,3.9001,3.6812,3.4746,3.2796,3.0955,2.9218,2.7578,2.5,2.3,2.1,1.9,1.7,1.5,1.3]#96
wavfil = 'aaaonkai.wav'
fs ,flow= wavfile.read(wavfil,mmap=False)
flow2=[]
basho=51
for i in range(part):
sa=flow[1102*basho+i+1]-flow[1102*basho+i]
zero=0
flow2.append(flow[1102*basho+i])
for j in range(0):#データの間を埋めたい時に使用。
zero=zero+sa/4
flow2.append(flow[1102*basho+i]+zero)
if i==10:
print(flow2)
gosasum=[]
gosaa=[0]*107
gosab=[0]*107
a_dush_=[0]*107
b_dush_=[0]*107
with open("output2.csv") as f:
reader=list(csv.reader(f))
#print(reader)
a_dush_=[float(i)for i in reader[0]]
b_dush_=[float(i)for i in reader[1]]
for i in range(107):
float(a_dush_[i])
float(b_dush_[i])
det1=np.zeros((214,214))
det2=np.zeros((214,1))
print("Wait a little.")
for i in range(part):#*4
for j in range(107):
sinij=np.sin(2*np.pi*i/heltz[j])
cosij=np.cos(2*np.pi*i/heltz[j])
#最急降下法による計算。誤差関数の計算。
#gosaa[j]=gosaa[j]+a_dush_[j]*cosij*cosij-b_dush_[j]*sinij*cosij-flow2[i]*cosij*volume
#gosab[j]=gosab[j]-(a_dush_[j]*cosij*sinij-b_dush_[j]*sinij*sinij-flow2[i]*sinij*volume)
#行列解法
for l in range(107):
det1[l,j]+=cosij*cosij
det1[l+107,j]+=cosij*sinij
det1[l,j+107]+=-sinij*cosij
det1[l+107,j+107]+=-sinij*sinij
det2[j,0]+=flow2[i]*cosij*volume
det2[j+107,0]+=-flow2[i]*sinij*volume
#print(det1[1,1])
#行列解法。逆行列の計算。
a=det1**-1 * det2
#print("Hello")
for i in range(107):
a_dush_[i]=a[i,0]
b_dush_[i]=a[i+107,0]
gosasum.append(gosaa)
gosasum.append(gosab)
alpha=myu
absum=0
for i in range(2):
for j in range(107):
absum=absum+abs(gosasum[i][j])
print(absum)
#最急降下法による計算。残差によるパラメータの更新。
#for h in range(107):
#a_dush_[h]=a_dush_[h]-alpha*gosaa[h]*reverse
#b_dush_[h]=b_dush_[h]-alpha*gosab[h]*reverse
print("Calculate ended.")
#print(a_dush,b_dush)
with open("output2.csv","w") as f2:
writer=csv.writer(f2)
writer.writerow(a_dush_)
writer.writerow(b_dush_)
#jijou()
create_wave(A, f0, fs, sec)
for i in range(repeat):
#create_wave0(A, f0, fs, sec)
create_wave(A, f0, fs, sec)
jijou()
print(i)fall.py
$${a_k}$$と$${b_k}$$によって予測式$${/sigma_{k=1}^n a_k cos(2πfx_i+b_k)}$$ 、$${a’_k=a_k cos b_k , b’_k=a_k sin b_k}$$と置くことで、$${a’_k ,b’_k}$$による予測式と、実際のデータとの誤差関数を求めて、パラメータ$${a',b'}$$を更新します。プログラムを反復して実行し、ある程度収束して値が変わらなくなったら、元のa,bを求めます。
ただし、今回は記事にある通り、行列を解くことでa'b'を求めることができることに気付いたので、途中からプログラムに書き足して、行列計算による求め方を追加しました。
コードには両方の解き方を載せてあります。両者で若干最終解の数値が違うような気がしますが、ほとんどは同じです。最急降下法によるフィッティングは50回以上繰り返してプログラムを実行しますが、行列解は一回の計算で終わります。200行程度変数を入力しても1分くらいで終わるので、圧倒的に早いです。解の再現性も高いと思われます。
上のwav音源、イヤホン注意⚠、不快音です。行列解の上のグラフを、音源としてそのまま出力したデータです。僕の声をそのまま分析用に使ったのですが、そのまま出力される訳ではありません。もう少し分析が必要なようです。
以上。
次の記事はこちら。