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)

平面波の場合，AMVののm番目の要素 a_m は，次式のようにかけます．

(2)

ここで，k は波数ベクトルと呼ばれ，次式で定義されます．

(3)

u は音源の方向ベクトルで，２次元の場合は，次式のようになります．

(4)

p_mは，m番目のマイクロホンの位置ベクトルです．

表１記号の説明

図１：音源方向とベクトル

アレイ・マニフォールド・ベクトル生成モジュール

このモジュールは，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)

サンプルプログラム

このプログラムでは、図２に示す音源配置（θ=-45°、r=1m）に対して、Kinectで収録したデータを読み込み（図３）、チャネル間位相差を表示します（図４）。また、計算により求めたAMVの位相も同様に図示し（図５）、実測したチャネル間位相差と比較します。

# チャネル間位相差と理論値との比較

# 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()

図２：音源配置（赤、θ=-45°)

図３：観測波形（z[0]）

図４：観測波形のチャネル間位相差

図５：位相差の理論値(アレイ・マニフォールド）

参考 - マルチチャネル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データ本体が記述されます。上述のように、本サイトでは、４byteのfloat型で記述します。下図は、データを格納する順番を示しています。z[m][t]は、m番目のチャネルの、時刻tにおけるデータを表します。通常のC++の２次元配列では、行内のデータが連続するようにデータが格納されますが、WAV形式では、時間の順序に重要な意味があり、新たな時刻のデータを追加／削除するのが容易なように、このようなデータの並びになっています。A/Dから観測データを受け取る場合も、ほとんどこのような並び方になっています。

WAVデータの並び型