Kinectの信号処理
Last update: 2018/08/22
ここでは、
Kinectで収録したデータをPythonで読み込む方法
Kinectでアレイ処理を行う上で必要となるアレイ・マニフォールド・ベクトル(AMV)の生成
について、解説します。
Kinectで収録したデータの読み込み
Microsftから提供されているAudioCaptureRawでデータを収録した場合,データは,このページの最後に記述してある多チャネルに拡張されたwav形式で保存されます.この拡張形式は,scipy.io.wavfileで提供されているwav.read()関数では,読み込めません.下記のmultiWavRead()関数を用いて,データを読み込んでください.
Kinectのデータ読み出しモジュール
kinectSensor.py
#######################################
# Function: multiWavRead
#######################################
def multiWavRead(fileName,dataType):
# Opne file
try:
fo = open(fileName,"rb")
except:
print( fileName," not found" )
sys.exit()
# Chunk
chunkID = fo.read(4)
chunkSize = np.fromfile(fo,dtype=np.int32,count=1)
Format = fo.read(4)
print( "Chunk ID: ",chunkID,", Chunk Size: ",chunkSize, ", Format:",Format )
# Subchunk 1
subChunk1ID = fo.read(4)
subChunk1Size = np.fromfile(fo,dtype=np.int32,count=1)[0]
formatTag,nChannel = np.fromfile(fo,dtype=np.int16,count=2)
Fs = np.fromfile(fo,dtype=np.int32,count=1)[0]
subChunk1Data = fo.read(subChunk1Size-8)
print( "SubChunk1 ID: ",subChunk1ID,", SubChunk1 Size: ",subChunk1Size )
print( "Number of Channek:",nChannel, ", Sampling Frequency:", Fs )
# Subchunk 2
subChunk2ID = fo.read(4)
subChunk2Size = np.fromfile(fo,dtype=np.int32,count=1)
print( "SubChunk2 ID: ",subChunk2ID,", SubChunks Size: ",subChunk2Size[0] )
if dataType == 'float':
buf = np.fromfile(fo,dtype=np.float32,count=int(subChunk2Size[0]/nChannel))
elif dataType == 'int':
buf = np.fromfile(fo,dtype=np.int32,count=int(subChunk2Size[0]/nChannel))
# dividing into channels
dataLength = int(len(buf)/nChannel)
x = buf.reshape((dataLength,nChannel))
print( "Data length: ",dataLength )
# Close file
fo.close()
return(x.T,Fs)
アレイ・マニフォールド・ベクトルの生成
アレイマニフォールドベクトル(AMV)は,図1に示す音源から各マイクロホンまでの直接音の伝達関数で構成され,自由空間の場合は次式のようになります(参考文献[1]).
(1)
平面波の場合,AMVののm番目の要素 am は,次式のようにかけます.
(2)
ここで,k は波数ベクトルと呼ばれ,次式で定義されます.
(3)
u は音源の方向ベクトルで,2次元の場合は,次式のようになります.
(4)
pmは,m番目のマイクロホンの位置ベクトルです.
表1 記号の説明
図1:音源方向とベクトル
アレイ・マニフォールド・ベクトル生成モジュール
このモジュールは,Kinect用のアレイ・マニフォールド・ベクトル(AMV)を生成します.行列Aの列ベクトルがそれぞれの音源方向のAMVに対応します.音源方向などは,関数の引数で指定します.Kinectのマイクロホン位置ベクトルsensorLocationを他のマイクロホンアレイのものに変えることで,他のマイクロホンアレイの場合にも使用可能です.
kinectSensor.py
#######################################
# Kinect sensor location
#######################################
sensorLocation = np.array([[113.0, -36.0, -76.0, -113.0],
[0.0, 0.0, 0.0, 0.0]])/1000.0
#######################################
# Function: getAMatrix
#######################################
def getAMatrix(direction,sensorLocation,Fs,NFFT,SoundVelocity):
Ndim,M = sensorLocation.shape
Nfreq = int(NFFT/2) + 1
Ndirection = len(direction)
A = np.zeros((Nfreq,Ndirection,M),dtype='complex')
for n in range(Ndirection):
theta = direction[n]/180.0*np.pi
u = np.array([np.sin(theta), np.cos(theta)])
for k in range(1,Nfreq): # omitting "0"
frequency = float(Fs)*float(k)/float(NFFT)
waveLength = SoundVelocity/frequency
for m in range(M):
p = sensorLocation[:,m]
up = np.dot(u,p)
A[k,n,m] = np.exp( 1j*2.0*np.pi/waveLength*up)
a0 = A[k,n,0]
for m in range(M):
A[k,n,m] = A[k,n,m]/a0
return A
参考文献
浅野太、"音のアレイ信号処理"、コロナ社 (2011)
サンプルプログラム
このプログラムでは、図2に示す音源配置(θ=-45°、r=1m)に対して、Kinectで収録したデータを読み込み(図3)、チャネル間位相差を表示します(図4)。また、計算により求めたAMVの位相も同様に図示し(図5)、実測したチャネル間位相差と比較します。
# チャネル間位相差と理論値との比較
# 2018/8/22
# Library
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig
import scipy.fftpack as ft
import sys
sys.path.append('../lib')
import kinectSensor as kin
# Parameter
NFFT = 512 # FFTの点数
win = np.hamming(NFFT) # 窓関数
Naverage = 10 # クロススペクトルの平均回数
SoundVelocity = 340.0 # 音速
targetDirection = -45.0 # 音源方向
flagFigSave = False
# 観測値の読み込み
observationFile = '../Data/Kinect/Speech0_Noise-45.wav' # 音源数=2のデータ
dataType = 'int'
z, Fs = kin.multiWavRead(observationFile,dataType) # Kinectデータの読み込み
(M,dataLength) = z.shape
# 観測値の表示
fig=plt.figure()
tm = np.arange(dataLength)/Fs
plt.plot(tm,z[0])
plt.xlabel('Time [s]')
(ymin,ymax) = plt.ylim() # チャネル間位相差を測定するブロック
BlockStart = 0/Fs
BlockEnd = NFFT*(Naverage+1)/Fs
plt.vlines(BlockStart,ymin,ymax,colors='g',linestyles='--')
plt.vlines(BlockEnd, ymin,ymax,colors='g',linestyles='--')
if flagFigSave:
plt.savefig('Graph/kinectObservation.png',format='png')
# クロススペクトル法によるチャネル間位相差
H = np.zeros((M,NFFT),dtype=complex)
for m in range(M):
Wxy = np.zeros(NFFT,dtype=complex)
Wxx = np.zeros(NFFT,dtype=complex)
for frameNumber in range(Naverage):
t1 = NFFT*frameNumber
t2 = NFFT*frameNumber+NFFT
X = ft.fft(z[0][t1:t2]*win,NFFT)
Y = ft.fft(z[m][t1:t2]*win,NFFT)
Wxy += np.conj(X)*Y
Wxx += np.conj(X)*X
H[m] = Wxy/Wxx
Nfreq = int(NFFT/2)+1
freq = Fs/NFFT*np.arange(Nfreq)
plt.figure()
for m in range(M):
plt.subplot(2,2,m+1)
plt.plot( freq,np.angle(H[m][0:Nfreq]) )
plt.title('Channel:'+str(m))
plt.grid()
if m == 2 or m == 3:
plt.xlabel('Frequency [Hz]')
plt.tight_layout()
if flagFigSave:
plt.savefig('Graph/kinectPhase.png',format='png')
# 位相差の理論値
# アレイ・マニフォールドの計算
direction = np.array([targetDirection])
AA = kin.getAMatrix(direction,kin.sensorLocation,Fs,NFFT,SoundVelocity)
Nfreq = int(NFFT/2)+1
G = np.zeros((M,Nfreq),dtype=complex)
for m in range(M):
G[m] = AA[:,0,m]
ref = np.copy(G[0]) # 0チャネルで正規化
for m in range(M):
G[m] /= ref
# 位相差の理論値のグラフ
plt.figure()
freq = Fs/NFFT*np.arange(Nfreq)
for m in range(M):
plt.subplot(2,2,m+1), plt.plot(freq,np.angle(G[m]))
plt.title('Channel:'+str(m))
plt.ylim(-np.pi,np.pi)
plt.grid()
if m == 2 or m == 3:
plt.xlabel('Frequency [Hz]')
plt.tight_layout()
if flagFigSave:
plt.savefig('Graph/kinectManifold.png',format='png')
# センサー配置
plt.figure(figsize=(5,3))
plt.plot(kin.sensorLocation[0],kin.sensorLocation[1],'bo')
r = 1.0
theta = targetDirection
target_x = r*np.sin(theta/180.0*np.pi)
target_y = r*np.cos(theta/180.0*np.pi)
plt.plot(target_x,target_y,'dr')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.xlim(-1.0,1.0)
plt.ylim(-0.3,1.0)
plt.grid()
if flagFigSave:
plt.savefig('Graph/kinectConfiguration.png',format='png')
plt.show()
図2:音源配置(赤、θ=-45°)
図3:観測波形(z[0])
図4:観測波形のチャネル間位相差
図5:位相差の理論値(アレイ・マニフォールド)
参考 - マルチチャネルWAVファイル形式 (C++)
ここでは、マイクロホンアレイで観測したデータを保存する形式について述べます。観測データの保存形式は任意のものでかまいませんが、ここでは、マイクロソフト社のマルチチャネルWAVフォーマットを使います。このような一般に流通しているフォーマットを用いることにより、他の人とデータの交換などをする際に便利です。また、下記で述べるWAVヘッダーに、チャネル数など、データに関する属性が記述されている点も利点です。
RIFF形式
WAV形式は、Microsoft社の提唱するRIFF形式と呼ばれるデータフォーマットの1つです。RIFF形式は、下図のようなChunk(かたまり)とその下位構造のSubchunkに分かれています。WAV形式では、Subchunk1にチャネル数やサンプリング周波数を含んだWAVフォーマット(属性)を、Subchunk2にWAVデータを、それぞれ記述します。
RIFF形式
Subchunk1 - WAVフォーマット
WAVフォーマットの部分は、次に示すWAVEFORMATEX構造体とWAVEFORMATEXTENSIBLE構造体の入れ子構造になっています。これらの構造体は、Microsoft社の提供するmmreg.hに定義されています。WAVEFORMATEX構造体は、チャネル数やサンプリング周波数など、WAVファイルの基本フォーマットが記述されています。WAVEFORMATEXTENSIBLE構造体の残りの部分は、追加フォーマットが記述されており、このうちのSubFormatに、KSDATAFORMAT_SUBTYPE_IEEE_FLOATを指定することで、データを-1〜+1の浮動小数点(float型)で記述する設定となります。WAVEFORMATEX構造体のcbSizeには、この部分のサイズを指定します。すなわち、sizeof(Samples) + sizeof(dwChannelMask) + sizeof(GUID) = 2 + 4 + 16 = 22。
WAVEFORMATEX構造体
typedef struct {
WORD wFormatTag; // WAVE_FORMAT_EXTENSIBLEに設定
WORD nChannels; // チャネル数
DWORD nSamplesPerSec; // サンプリング周波数 [Hz]
DWORD nAvgBytesPerSec; // nSamplesPerSec x nBlockAlign [byte/sec]
WORD nBlockAlign; // nChannel x wBitsPerSample / 8 [byte]
WORD wBitsPerSample; // データサンプルのビット数. 32に設定
WORD cbSize; // 追加フォーマットのサイズ [byte]. 22に設定
} WAVEFORMATEX;
WAVEFORMATEXTENSIBLE構造体
typedef struct { WAVEFORMATEX Format; union { WORD wValidBitsPerSample; // WAVEFORMATEX.wBitsPerSampleを指定 WORD wSamplesPerBlock; // 未使用 WORD wReserved; // 未使用 } Samples; DWORD dwChannelMask; // 未使用 GUID SubFormat; // KSDATAFORMAT_SUBTYPE_IEEE_FLOATを指定 } WAVEFORMATEXTENSIBLE
Subchunk2 - WAVデータ
Subchunk2にはWAVデータ本体が記述されます。上述のように、本サイトでは、4byteのfloat型で記述します。下図は、データを格納する順番を示しています。z[m][t]は、m番目のチャネルの、時刻tにおけるデータを表します。通常のC++の2次元配列では、行内のデータが連続するようにデータが格納されますが、WAV形式では、時間の順序に重要な意味があり、新たな時刻のデータを追加/削除するのが容易なように、このようなデータの並びになっています。A/Dから観測データを受け取る場合も、ほとんどこのような並び方になっています。
WAVデータの並び型