音波をフーリエ変換と少しだけ違う方法で、好きな周波数だけで近似して、音波の構成を求める(試作)
前回はこちら。
前回に引き続き、音波解析用のプログラムを自作してみます。
import numpy as np
import wave
import struct
import csv
from scipy.linalg import lu_factor,lu_solve,solve
import scipy.linalg as linalg
#パラメータ
#--------------------------------------------------------------------------------------------------------------------
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#1秒辺りのサンプル数
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.366,4.3777,4.1320,3.9001,3.6812,3.4746,3.2796,3.0955,2.9218,2.7578,2.6,2.5,2.4,2.3,2.1,2.05,2]#107種類のデータを求めたい周波数のセット。44100/heltz[i]が周波数であり、周波数の逆数をとっている。
#heltz=[1102.5,551.25,367.5,275.625,220.5,183.75,157.5,137.8125,122.5,110.25,100.227272727273,91.875,84.8076923076923,78.75,73.5,68.90625,64.8529411764706,61.25,58.0263157894737,55.125,52.5,50.1136363636364,47.9347826086957,45.9375,44.1,42.4038461538462,40.8333333333333,39.375,38.0172413793103,36.75,35.5645161290323,34.453125,33.4090909090909,32.4264705882353,31.5,30.625,29.7972972972973,29.0131578947368,28.2692307692308,27.5625,26.890243902439,26.25,25.6395348837209,25.0568181818182,24.5,23.9673913043478,23.4574468085106,22.96875,22.5,22.05,21.6176470588235,21.2019230769231,20.8018867924528,20.4166666666667,20.0454545454545,19.6875,19.3421052631579,19.0086206896552,18.6864406779661,18.375,18.0737704918033,17.7822580645161,17.5,17.2265625,16.9615384615385,16.7045454545455,16.455223880597,16.2132352941176,15.9782608695652,15.75,15.5281690140845,15.3125,15.1027397260274,14.8986486486486,14.7,14.5065789473684,14.3181818181818,14.1346153846154,13.9556962025316,13.78125,13.6111111111111,13.4451219512195,13.2831325301205,13.125,12.9705882352941,12.8197674418605,12.6724137931034,12.5284090909091,12.3876404494382,12.25,12.1153846153846,11.9836956521739,11.8548387096774,11.7287234042553,11.6052631578947,11.484375,11.3659793814433,11.25,11.1363636363636,11.025,10.9158415841584,10.8088235294118,10.7038834951456,10.6009615384615,10.5,10.4009433962264,10.303738317757]
a_dush=[0]*107
b_dush=[0]*107
with open("output2.csv") as f:#読み取り元のcsvデータファイルを読み込む。1行目と2行目の数値を使う。
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])
#--------------------------------------------------------------------------------------------------------------------
point = np.arange(0,fs*t)
sin_wave=[0]*fs*t
for n in range(107):#波形を作る。a cos(nt)+b cos(nt)の全ての周波数の和を取る。
sin_wave =sin_wave +A*a_dush[n]*np.cos(2*np.pi*point/heltz[n])-A*a_dush[n]*np.sin(2*np.pi*point/heltz[n])
max=0
for i in range(fs*t):
if abs(sin_wave[i])>max:
max=abs(sin_wave[i])
sin_wave = [int(x * 32767.0/max/4) for x in sin_wave]
print(sin_wave[:10])
binwave = struct.pack("h" * len(sin_wave), *sin_wave)
w = wave.Wave_write("output.wav")#波形をデータファイルに書き込み。
p = (1, 2, fs, len(binwave), 'NONE', 'not compressed')
w.setparams(p)
w.writeframes(binwave)
w.close()
from scipy.io import wavfile
import math
import array
import csv
def jijou():#全周波数を合わせた合成波と、一つの周波数のみに焦点を当てた単調な正弦波について、元の波形に一番近似する波形の振幅を求めてcsvファイルに出力する。
wavfil = 'aaaonkai.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.366,4.3777,4.1320,3.9001,3.6812,3.4746,3.2796,3.0955,2.9218,2.7578,2.6,2.5,2.4,2.3,2.1,2.05,2]#96
#heltz=[1102.5,551.25,367.5,275.625,220.5,183.75,157.5,137.8125,122.5,110.25,100.227272727273,91.875,84.8076923076923,78.75,73.5,68.90625,64.8529411764706,61.25,58.0263157894737,55.125,52.5,50.1136363636364,47.9347826086957,45.9375,44.1,42.4038461538462,40.8333333333333,39.375,38.0172413793103,36.75,35.5645161290323,34.453125,33.4090909090909,32.4264705882353,31.5,30.625,29.7972972972973,29.0131578947368,28.2692307692308,27.5625,26.890243902439,26.25,25.6395348837209,25.0568181818182,24.5,23.9673913043478,23.4574468085106,22.96875,22.5,22.05,21.6176470588235,21.2019230769231,20.8018867924528,20.4166666666667,20.0454545454545,19.6875,19.3421052631579,19.0086206896552,18.6864406779661,18.375,18.0737704918033,17.7822580645161,17.5,17.2265625,16.9615384615385,16.7045454545455,16.455223880597,16.2132352941176,15.9782608695652,15.75,15.5281690140845,15.3125,15.1027397260274,14.8986486486486,14.7,14.5065789473684,14.3181818181818,14.1346153846154,13.9556962025316,13.78125,13.6111111111111,13.4451219512195,13.2831325301205,13.125,12.9705882352941,12.8197674418605,12.6724137931034,12.5284090909091,12.3876404494382,12.25,12.1153846153846,11.9836956521739,11.8548387096774,11.7287234042553,11.6052631578947,11.484375,11.3659793814433,11.25,11.1363636363636,11.025,10.9158415841584,10.8088235294118,10.7038834951456,10.6009615384615,10.5,10.4009433962264,10.303738317757]
wavfil ='440hz.wav'#音源のwavファイルの読み込み。
fs ,flow= wavfile.read(wavfil,mmap=False)
flow2=[]
basho=20
for i in range(part):
sa=flow[1102*basho+i+1]-flow[1102*basho+i]
zero=0
flow2.append(flow[1102*basho+i])
if i==10:
print(flow2)
gosasum=[]
gosaa=[0]*107
gosab=[0]*107
a_dush_=[0]*107
b_dush_=[0]*107
print("Wait a little.")
sinlist=[[np.sin(2*np.pi*i/heltz[j]) for i in range(part)]for j in range(107)]
coslist=[[np.cos(2*np.pi*i/heltz[j]) for i in range(part)]for j in range(107)]
l=0
for i in range(part):#①求めたい周波数の波を全て含む合成波で元の音波に最も近いものを作った場合のそれぞれの周波数毎の振幅を全て求める。(使いづらいので以後の分析では使わないダミーのプログラムになりました。)
l+=1
if l==30:
print(i)
l=0
coscos=[[coslist[j][i]*coslist[k][i]for j in range(107)]for k in range(107)]
cossin=[[sinlist[j][i]*coslist[k][i]for j in range(107)]for k in range(107)]
sinsin=[[sinlist[j][i]*sinlist[k][i]for j in range(107)]for k in range(107)]
for j in range(107):
for k in range(107):
det1[k,j]+=coscos[j][k]
det1[k+107,j]+=cossin[j][k]
det1[k,j+107]-=cossin[k][j]
det1[k+107,j+107]-=sinsin[j][k]
det1[214,j]+=coslist[j][i]
det1[214,j+107]-=sinlist[j][i]
det1[j,214]+=coslist[j][i]
det1[107+j,214]+=sinlist[j][i]
det1[214,214]+=1
det2[j,0]+=flow2[i]*coslist[j][i]*volume
det2[j+107,0]+=flow2[i]*sinlist[j][i]*volume
det2[214,0]+=flow2[i]*volume
print(det1[212,205:])
a=lu_solve(lu_factor(det1),det2)
print(a[107,0])
det1=np.zeros((215,215))
det2=np.zeros((215,1))
for i in range(107):
a_dush_[i]=a[i,0]
b_dush_[i]=a[i+107,0]
c_dush_=[0]*107
d_dush_=[0]*107
e_dush_=[0]*107
det3=np.zeros((3,3))
det4=np.zeros((3,1))
print("Wait a little.")
for j in range(107):#②求めたい周波数を一つだけ含む音波で元の音波に最も近いものを作った時の振幅を全ての求めたい周波数で計算する。
print(j)
coscos=[[coslist[j][i]*coslist[k][i]for i in range(part)]for k in range(107)]
cossin=[[sinlist[j][i]*coslist[k][i]for i in range(part)]for k in range(107)]
sinsin=[[sinlist[j][i]*sinlist[k][i]for i in range(part)]for k in range(107)]
sincos=[[sinlist[k][i]*coslist[j][i]for i in range(part)]for k in range(107)]
for i in range(part):#*4
det3[0,0]+=coscos[j][i]
det3[1,0]+=cossin[j][i]
det3[0,1]-=sincos[j][i]
det3[1,1]-=sinsin[j][i]
det3[2,0]+=coslist[j][i]
det3[2,1]-=sinlist[j][i]
det3[0,2]+=coslist[j][i]
det3[1,2]+=sinlist[j][i]
det3[2,2]+=1
det4[0,0]+=flow2[i]*coslist[j][i]*volume
det4[1,0]+=flow2[i]*sinlist[j][i]*volume
det4[2,0]+=flow2[i]*volume
b=lu_solve(lu_factor(det3),det4)
c_dush_[j]=b[0,0]
d_dush_[j]=b[1,0]
e_dush_[j]=b[2,0]
det3=np.zeros((3,3))
det4=np.zeros((3,1))
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)
print("Calculate ended.")
with open("output2.csv","w") as f2:#①sin波の振幅を1行目、cos波の振幅を2行目。②sin波の振幅を3行目、cos波の振幅を4行目、切片cを5行目に書き出す。
writer=csv.writer(f2)
writer.writerow(a_dush_)
writer.writerow(b_dush_)
writer.writerow(c_dush_)
writer.writerow(d_dush_)
writer.writerow(e_dush_)
#jijou()
#create_wave(A, f0, fs, sec)
#jijou()
def jijou2():#jijouの②のデータから元の周波数に近い波を最小二乗法(ムーア・ペンローズの逆行列、一般化逆行列)で推定する。(…③)
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.366,4.3777,4.1320,3.9001,3.6812,3.4746,3.2796,3.0955,2.9218,2.7578,2.6,2.5,2.4,2.3,2.1,2.05,2]#96
#heltz=[1102,551,367.333333333333,275.5,220.4,183.666666666667,157.428571428571,137.75,122.444444444444,110.2,100.181818181818,91.8333333333333,84.7692307692308,78.7142857142857,73.4666666666667,68.875,64.8235294117647,61.2222222222222,58,55.1,52.4761904761905,50.0909090909091,47.9130434782609,45.9166666666667,44.08,42.3846153846154,40.8148148148148,39.3571428571429,38,36.7333333333333,35.5483870967742,34.4375,33.3939393939394,32.4117647058824,31.4857142857143,30.6111111111111,29.7837837837838,29,28.2564102564103,27.55,26.8780487804878,26.2380952380952,25.6279069767442,25.0454545454545,24.4888888888889,23.9565217391304,23.4468085106383,22.9583333333333,22.4897959183673,22.04,21.6078431372549,21.1923076923077,20.7924528301887,20.4074074074074,20.0363636363636,19.6785714285714,19.3333333333333,19,18.6779661016949,18.3666666666667,18.0655737704918,17.7741935483871,17.4920634920635,17.21875,16.9538461538462,16.6969696969697,16.4477611940299,16.2058823529412,15.9710144927536,15.7428571428571,15.5211267605634,15.3055555555556,15.0958904109589,14.8918918918919,14.6933333333333,14.5,14.3116883116883,14.1282051282051,13.9493670886076,13.775,13.6049382716049,13.4390243902439,13.2771084337349,13.1190476190476,12.9647058823529,12.8139534883721,12.6666666666667,12.5227272727273,12.3820224719101,12.2444444444444,12.1098901098901,11.9782608695652,11.8494623655914,11.7234042553191,11.6,11.4791666666667,11.360824742268,11.2448979591837,11.1313131313131,11.02,10.9108910891089,10.8039215686275,10.6990291262136,10.5961538461538,10.4952380952381,10.3962264150943,10.2990654205607]
gosasum=[]
gosaa=[0]*107
gosab=[0]*107
a_dush_=[0]*107
b_dush_=[0]*107
c_dush_=[0]*107
det3=np.zeros((1,3))
det4=np.zeros((3,1))
with open("output2.csv") as f:#csvファイルの読み込み(3行目〜5行目を使う)
reader=list(csv.reader(f))
a_dush_=[float(i)for i in reader[2]]
b_dush_=[float(i)for i in reader[3]]
c_dush_=[float(i)for i in reader[4]]
for i in range(107):
float(a_dush_[i])
float(b_dush_[i])
float(c_dush_[i])
det1=np.zeros((214,214))
det2=np.zeros((214,107))
print("Wait a little.")
sinlist=[[np.sin(2*np.pi*i/heltz[j]) for i in range(part)]for j in range(107)]
coslist=[[np.cos(2*np.pi*i/heltz[j]) for i in range(part)]for j in range(107)]
msinlist=[[np.sin(np.pi*i/heltz[j]) for i in range(part)]for j in range(107)]
mcoslist=[[np.cos(np.pi*i/heltz[j]) for i in range(part)]for j in range(107)]
for j in range(107):
print(j)
for i in range(part):#②のデータから作る推定用の音波データと、推定したい振幅を指定する変数を含むデータを行列に入力する。(ここの作業のスピードはかなり重いですが、複数の行列の積の形にすることで、もう少し軽量化を図れると思います。たぶん。。)
for k in range(107):
det1[k,j]+=coslist[j][i]*mcoslist[k][i]
det1[k,j+107]-=sinlist[j][i]*mcoslist[k][i]
det1[k+107,j]+=coslist[j][i]*msinlist[k][i]
det1[k+107,j+107]-=coslist[j][i]*mcoslist[k][i]
det2[k,j]+=a_dush_[j]*coslist[j][i]*mcoslist[k][i]-b_dush_[j]*sinlist[j][i]*mcoslist[k][i]+c_dush_[j]*mcoslist[k][i]
det2[k+107,j]+=a_dush_[j]*coslist[j][i]*msinlist[k][i]-b_dush_[j]*sinlist[j][i]*msinlist[k][i]+c_dush_[j]*msinlist[k][i]
detX=np.zeros((214,107))
for i in range(107):
detX[i,i]=1
detX[107+i,i]=1
print(det1[1,1])
detY=np.ones((1,107))
try:
a=linalg.pinv(det1)@det2*detX@linalg.pinv(detY)#ムーア・ペンローズの逆行列、アダマール積、普通の行列積を使い、求めたい周波数の振幅を求める変数を計算する。
except:
print('pinv errored inverse used.')
Q,R=linalg.qr(det1.T@det1)
#a=linalg.solve(R,np.dot(Q.T,det1.T@det2*detX@linalg.pinv(detY)))
# a=linalg.inv(det1)@det2*detX@linalg.pinv(detY)
#a=solve(det1.T@det1,det1.T@det2*detX@linalg.pinv(detY))
giji=(det1.T @ det1)**-1@det1.T
a=giji@det2*detX@linalg.pinv(detY)
#a=det1**-1 @ det2
#k=det1@a
print("Wait a little.")
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
print("Calculate ended.")
with open("output2.csv","a") as f2:#w or a
writer=csv.writer(f2)#結果をcsvファイルに書き込む。1行目にsin波の振幅、2行目にcos波の振幅。wを指定すれば新規、aを指定すれば追記で書き込みができる。書き込んだファイルはcreate_waveでそのまま音波を合成するのに使うことができる。
writer.writerow(a_dush_)
writer.writerow(b_dush_)
jijou()
jijou2()
print('calculate ended')
create_wave(A, f0, fs, sec)
fall.py
このプログラムは、主に3つの要素で構成されます
① 求めたい周波数の波を全て含む合成波で元の音波に最も近いものを近似する。
② 求めたい周波数を一つだけ含む音波で元の音波に最も近いものを近似する。
③②のデータと同じような聞こえ方になる(と思われる)合成波を逆算して作る。
③のプログラムは分析用というよりは、音源の合成用です。好きな聞こえ方になる音を
(無理矢理に)作るためのものです。
②を参考にして、③のプログラムの出力が決定されます。②は前回の記事のおまけに、最後に少しだけ書き込んでいたグラフとほぼ同じになります。何故②を参考にしたのかと言うと、②が一番内耳の音の聞こえる仕組みに近いグラフになると考えたからです。
耳では、外から入ってきた音を鼓膜で受け取り、耳小骨で増幅した上で、内耳のカタツムリ状の構造で感知します。内耳の中では、内耳の管の太さに応じた、固有の周波数の共鳴により音を感知して、内耳の神経への入力へと繋げています。フーリエ変換により見えるグラフは、おそらくは音を構成する波を、さまざまな周波数の波の合成と捉えた時の、振幅のスペクトルをそれぞれ分解して表示するものと思います。
僕はこれを考えた場合に、分解した波が必ずしも耳で感知する音の大きさと一致するとは限らない。むしろ、もっと直感的に、耳の中の共鳴の仕方に近いように分解すれば、分解した後の波を再合成した時に、好きな周波数で似たような聞こえ方のする音が作れるのではないかと思いました。また、②を参考にしてもう一度いろんな周波数での波の再合成することを考えれば、例えば声の周波数の細かい構成を、フーリエ変換よりも柔軟に分析したり、再構成することができるのではと考えました。
上のプログラムが試行錯誤の結果となります。
実際に同じ音を分析してみた結果のグラフは以上のようになります。1枚目の3つのグラフが今回のプログラムでできた振幅の分布データとなります。それぞれの形に若干の違いがあるものの、ほとんどの周波数で同じ程度の大きさであり、明確な違いがあるようにはあまり見えません。
2枚のスクリーンショットと同じグラフを、縦軸のスケールを線形軸から対数軸に切り替えてもう一度見てみます。
ほとんど同じですが、細かい音の大きさはFFT(高速フーリエ変換：pythonに標準で備え付けられている関数)とそれ以外とで少し違いが見られます。僕のプログラムが間違っている可能性とFFTの見方をある程度間違えている可能性が考えられます。どっちがどう違うのかがなかなか指摘しにくいため、どちらをより信じればいいかはわかりません。
音波のエネルギーに関する運動方程式から、振幅とエネルギーとの関係は高周波になればなるほど、振幅の割に多くのエネルギーを抱えていることになります。また、FFTではある数の自然数倍の周波数での振幅をそれぞれ調べています。今回作ったデータと比べると、高周波になればなるほど、FFTの方がたくさんの周波数を調べて、スペクトルを算出しています。
それがFFTと今回作ったデータの違いに影響があるのかもしれないと考えています。しかし、もう少し調べてみないとわかりません。
後は位相の話。細かくチェックしないとわからない微妙な変化の話ではありますが、今回のプログラムで近似させたデータについて、位相を表すデータはarctan(3行目/4行目)で計算することができますので、計算した結果を見ると、ちょっと上のスクリーンショットではわかりにくくはありますが、周波数が大きくなるごとに、正の値と負の値を行ったり来たりを繰り返しているように見えます。これは、同じ波形を見ている場合でも、周波数を漸増させていけばいった分だけ、近似させた波の位相は同方向に連続的にずれていっていくものであることを示していると思われます。
この性質は、フーリエ変換によって表した波の位相についても同様のものと思われます。読んでる方でFFTのプログラムを使うことができる人がいれば、absで振幅を取る前のデータを取り出して、確認してみてください。たぶん同じ特徴が目視できると思います。
音源です。ご参考までに。
以上です。
この記事が気に入ったらサポートをしてみませんか？