Wavelet変換とデータ圧縮(1)

Masahiro MIZUTANI
11 Nov  2011 2012 Nov 29

データ数が2べき乗 L=(2^n)のモノラル音声信号をwaveletを使ってデータ圧縮する方法を考えてみる。
Wavelet変換は局所的平均を次々と求めていく操作と考えることができる（実際、Haar wavelet変換はscaleling係数として隣り合う2つのデータ（またはデータブロック）の平均を、wavelet係数としてその平均からのズレを次々と求める操作だ）。

単純な非可逆データ圧縮として考えられる方法として、ここでは元データをレベル p のwavelet変換で得られる平均値として近似する方法を考える。
Haar Wavelet変換を使うレベル p 近似のデータ数は(1/(2^p))倍に減っている。レベルpであることを知って、そのレベル p 近似のデータで構成されるデータに「納得」できるなら圧縮に成功と見なせるはずだ。

高速Haar wavelet変換

[高速Wavelet変換アルゴリズム] Mathematicaの配列指定は 1 から始まることを考慮している：signal=

L=(2^J)の標本データ signal=((Subscript[c, 0]^J), (Subscript[c, 1]^J), (Subscript[c, 2]^J),...., (Subscript[c, L-1]^J)) を与えて、J回wavelet変換を繰り返して得られるリスト ((Subscript[c, 0]^((J))), (Subscript[d, 0]^((J-1))) , (Subscript[d, (Subscript[k, 2])]^((2))),...., (Subscript[d, (Subscript[k, i])]^((n)))) を返す
 
L=(2^J)
data = signal 
for j = J to 1
	for k = 0 to (2^(j-1))-1
 		tmp = data[(2^(J+1-j))k + 1]
 		data[(2^(J+1-j))k + 1]= (1/2) (tmp + data[(2^(J+1-j))k + (2^(J-j)) + 1])
  		data[(2^(J+1-j))k + (2^(J-j)) + 1] = (1/2) (tmp - data[(2^(J+1-j))k + (2^(J-j)) + 1])
dataを返り値とする

FirstHaarWaveletTransform[signal_List]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[signal]];
data=signal;
For[j=scale,j≥1,j--,
For[k=0,k≤2^(j-1) -1,k++,
tmp=data[[2^{scale+1-j} k+1]];
data[[2^{scale+1-j} k+1]]=(tmp+data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
data[[2^{scale+1-j} k+2^(scale-j)+1]]=(tmp-data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
]
];
data
]

wavelet変換レベル 0≤level≤scale 毎のスケール係数

関数 f のスケール J の標本デー(Subscript[タf, J])（(2^J)(Subscript[個の関数値c, k]^J)）に対して、Haarl近似（隣り合うデータを足して2で割る平均操作）をJ-level回（J≥level≥0）繰り返して得られる 近似関数 (Subscript[f, level])の（平均値）係数リスト（デー((タ数は2)^J)から (2^level)個に減少）

HaarScalingCoefficient[signal_List,level_Integer]:=Module[{scale,tmp,data,j,k,rec={}},
scale=Log[2,Length[signal]];
data=signal;
For[j=scale,j>level,j--,
For[k=0,k≤2^(j-1) -1,k++,
tmp=data[[2^{scale+1-j} k+1]];
data[[2^{scale+1-j} k+1]]=(tmp+data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
data[[2^{scale+1-j} k+2^(scale-j)+1]]=(tmp-data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
]
];
For[k=0, k≤2^level-1,k++,
AppendTo[rec,data[[2^(scale-level)  k+1]]];
];
rec
]

u={3,1,0,4,8,6,9,9};

変換Level pのHaar近似値

HaarScalingCoefficient[signal_List,level_Integer]

HaarScalingCoefficient[u,3]

HaarScalingCoefficient[u,2]

HaarScalingCoefficient[u,1]

HaarScalingCoefficient[u,0]

Haar wavelet逆変換

[高速逆Wavelet変換]
L=(2^J)の標本データ signal=((Subscript[c, 0]^J), (Subscript[c, 1]^J), (Subscript[c, 2]^J),...., (Subscript[c, L-1]^J)) をJ回wavelet変換を繰り返して得られるリスト ((Subscript[c, 0]^((J))), (Subscript[d, 0]^((J-1))) , (Subscript[d, (Subscript[k, 2])]^((2))),...., (Subscript[d, (Subscript[k, i])]^((n)))) が与えられたとして、逆変換を繰り返して「元データ」を返す。 Mathematicaの配列指定は 1 から始まることを考慮している：

L=(2^J)
data = waveletdata
for j = 0 to J-1
	 for k = 0 to  (2^j)-1
 		tmp = data[(2^(J-j))k + 1]
 		data[(2^(J-j))k + 1] =tmp + data[(2^(J-j))k + (2^(J-j-1))+1]
  		data[(2^(J-j))k + (2^(J-j-1))+1] = tmp - data[(2^(J-j))k+ (2^(J-j-1))+1]
dataを返り値とする

FirstHaarInverseWaveletTransform[waveletdata_List]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[waveletdata]];
data=waveletdata;
For[j=0,j<scale,j++,
For[k=0,k<=2^j -1,k++,
tmp=data[[2^(scale-j) k+1]];
data[[2^(scale-j) k+1]]=tmp + data[[2^(scale-j) k+2^(scale-j-1)+1]];
data[[2^(scale-j) k+2^(scale-j -1)+1]]=tmp-data[[2^(scale-j) k+2^(scale-j-1)+1]];
]
];
data
]

レベル指定のwavelet変換

レベルpまでのwavelet変換

((長さを2)^J)のスケールJの信号 (Subscript[f, J])（信号係数{(Subscript[c, k]^J)}）をlevel Jとする。レベル J≥p≥0((を指定して得られる2)^J)個のHaar scaling/wavelet係数のリスト{(Subscript[c, k]^j),(Subscript[d, k]^j)}を求める。

HaarWaveletTransfLevel[signal_List,p_Integer]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[signal]];
data=signal;
For[j=scale,j>p,j--,
For[k=0,k≤2^(j-1) -1,k++,
tmp=data[[2^{scale+1-j} k+1]];
data[[2^{scale+1-j} k+1]]=(tmp+data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
data[[2^{scale+1-j} k+2^(scale-j)+1]]=(tmp-data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
]
];
data
]

HaarWaveletTransfLevel[u,1]

HaarWaveletTransfLevel[u,2]

HaarWaveletTransfLevel[u,3]

レベルpからのwavelet逆変換

スケー((ルJの信号の長さを2)^J)としたとき、元信号{(Subscript[c, k]^J)}をlevel Jとして,レベル J≥p≥0((までの2)^J)個のHaar scaling/wavelet係数のリスト{(Subscript[c, k]^j),(Subscript[d, k]^j)}から、元信号{(Subscript[c, k]^J)}を求める

InvHaarWaveletTransfLevel[waveletleveldata_List,p_Integer]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[waveletleveldata]];
data=waveletleveldata;
For[j=p,j<scale,j++,
For[k=0,k<=2^j -1,k++,
tmp=data[[2^(scale-j) k+1]];
data[[2^(scale-j) k+1]]=tmp + data[[2^(scale-j) k+2^(scale-j-1)+1]];
data[[2^(scale-j) k+2^(scale-j -1)+1]]=tmp-data[[2^(scale-j) k+2^(scale-j-1)+1]];
]
];
data
]

InvHaarWaveletTransfLevel[HaarWaveletTransfLevel[u,1],1]

InvHaarWaveletTransfLevel[HaarWaveletTransfLevel[u,2],2]

InvHaarWaveletTransfLevel[HaarWaveletTransfLevel[u,3],3]

音声データを読み込む

音声は例えば媒質の圧縮と膨脹の波が伝わって我々に到達することで認知されるが、音圧によって板を振動させて連続的な電流を生じさせるマイクからアナログ電気信号に変換して伝達することができる。
アナログ信号をCDやDVDなにに記録されるデジタル信号に変換するためには
	アナログ信号をサンプリング（標本化）する
	連続値を取る信号値を量子化（階段化）する
という2つの処理が必要である。
CDの場合、1秒間に44100回サンプリングし（サンプリング周波数=44100Hz）、得られた最大振幅を(2^16)=65536段階(16bit=2byte)に量子化している。サンプリング周波数の1/2の周波数まで再現できる（標本化定理）ため、CDでは周波数22050Hzを上限とした音声データが収録されている。
Stereo CDデジタルデータは理論的には
	1秒間で 44100 Hz × 2byte × 2 =176400 byte ≈172KB/sec
	1分間で 172KB/sec × 60 sec = 10320KB≈10MB
のデータサイズが必要になる。CD 1枚に収納できる出来るデータ量は最大約700MBであるため、70分のStereo 音声データを記録することができる。

Wave ファイル（モノラル）を読み込む

以下ではWavが置かれたフォルダに移動し、Sound_of _Silence.wavを読み込むとしている

SetDirectory["~/Desktop"]

通常の音楽データの読み込みはオプションなしでインポート関数 Import[ ] をつかう

sound=Import["Sound_of_Silence.wav"]

wavファイルのサンプリングデータの取得

音楽数値データの読み込みは、オプション “Data” を指定して Import[ ]

dsound=Import["Sound_of_Silence.wav","Data"]

sounddata=First[dsound];

len=Length[sounddata]

2^(19)

2^(18)

level pの音声信号のHaar Wavelet近似

近似level p の範囲 1 ≤ p ≤ n

((CDのサンプリング周波数44100Hzの長さ2)^19)=524288の音声データを用意

(2^19)/44100 =11.888 秒分の音声信号

signal=Take[sounddata,2^(19)];

Log[2,Length[signal]]

サンプリング周波数を合わせてHaar Wavelet近似で音声を再生

approx18=HaarScalingCoefficient[signal,18];

Length[approx18]

approxSound1=Sound[SampledSoundList[List[approx18],44100/2]]

たいへん効率の悪い近似度音声の生成

毎回、標本データに対するwavelet変換を繰り返している。

farlevel = scale (= Log[2, Length[signal])が原音、level=1,2,.. で平均化が進む

SoundApproximationLevel[signal_List,level_Integer]:=Module[{scale,approx,j},
scale=Log[2,Length[signal]];
Print[farlevel,"-th Haar Wavelet Approximation:","  Data size=",Length[signal]/2^level];
Print["Sampling Frequency=",Floor[44100/2^(scale-level)]];
approx=HaarScalingCoefficient[signal,scale-level];
Sound[SampledSoundList[List[approx],Floor[44100/2^level]]]
]

approxSound2=SoundApproximationLevel[signal,2]

approxSound3=SoundApproximationLevel[signal,3]

approxSound4=SoundApproximationLevel[signal,4]

approxSound5=SoundApproximationLevel[signal,5]

approxSound6=SoundApproximationLevel[signal,6]