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

Masahiro MIZUTANI
22 Nov  2011

スケールJのデータ数 L=(2^J)個の信号データ (Subscript[f, J]) を圧縮する方法として、微少なwavele係数を無視して（0として考え）データ量を縮約する方法がある。
詳細wavelet係数は局所平均からのズレである。そのwavelet係数が「無視し得る」程度に微少であれば、それら0で置き換えて修正されたwaveletデータとし、これをwavelet逆変換して「元データとして復元する」非可逆データ圧縮である。

Haar waveletを音楽データに適用した場合、レベルJからレベルp  (1 ≤p ≤ J)の平均値で置き換えた近似データ (Subscript[f, p])による圧縮法（(2^(p-J))倍に減少）とHaar waveletの微少なwavele係数を無視する方法のどちらが音楽情報の縮約方法として適切かの評価は難しい。Daubechies Wavelet変換との比較が、今後課題となる。

最大振幅に比べて相対的に微少なデータを別な記号に置き換える

リストdataの要素の絶対値が、data内の要素で絶対値が最大 maxabs の thresholdRaio*maxabs以下 Abs[x] < thresholdRaio*maxabs のであれば、その要素を記号  replacedSymbol で置き換えたリストを返す

ModifiedListNeglectingRelativeSmallElement[data, thresholdRatio, replaceSymbol]
   [引数] data: 数リスト,   thresholdRatio: 閾値の割合(実数),   replaceSymbol: 置き換え記号

ModifiedListNeglectingRelativeSmallElement[data_List,thresholdRatio_,replaceSymbol_]:=Module[{maxabs},
maxabs=Max[Abs[data]];
data /. x_ /;Abs[x]<thresholdRatio*maxabs ->replaceSymbol
]

ModifiedListNeglectingRelativeSmallElement[{2,4,-5,1,6,-6,9,8,-2},0.3,X]

高速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
]

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

FirstHaarWaveletTransform[u]

高速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
]

FirstHaarInverseWaveletTransform[FirstHaarWaveletTransform[u]]

音声データを読み込む

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

SetDirectory["~/Desktop"]

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

通常の音楽データの読み込みはオプションなしでインポート関数 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)

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

(2^18)/44100 =5.944 秒分の音声信号

計算で使用する元データの設定

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

Wavelet係数

Wavelet係数リストの先頭値は全データの平均値 （Scaling係数） であり、実際のWavelet係数は残りのリストであることに注意。

signalWaveletData=FirstHaarWaveletTransform[signal];

計算チェック（逆変換して元に戻るかを確かめる）

invWaveletSignalData=FirstHaarInverseWaveletTransform[signalWaveletData];

originalSound=Sound[SampledSoundList[List[invWaveletSignalData],44100]]

微少Waveletデータを無視して 「圧縮」 し、復元する

Wavelet係数リストの最大絶対値閾値の割合を35 % まで変化させて、未満のデータを0にして圧縮

最大絶対値の2% 未満を0にして圧縮

modified02=ModifiedListNeglectingRelativeSmallElement[signalWaveletData,0.02,X];

num=Count[modified02,X]

圧縮比

Print["Compression Raio for original volume= ", N[(Length[signal]-num)/Length[signal]]*100,"%"]

修正Waveletデータ

modifiedwavelet02=modified02/. X->0;

modifiedsignal02=FirstHaarInverseWaveletTransform[modifiedwavelet02];

Sound[SampledSoundList[List[modifiedsignal02],44100]]

最大絶対値の5% 未満を0にして圧縮

modified05=ModifiedListNeglectingRelativeSmallElement[signalWaveletData,0.05,X];

num=Count[modified05,X]

圧縮比

Print["Compression Raio for original volume= ", N[(Length[signal]-num)/Length[signal]]*100,"%"]

修正Waveletデータ

modifiedwavelet05=modified05/. X->0;

modifiedsignal05=FirstHaarInverseWaveletTransform[modifiedwavelet05];

Sound[SampledSoundList[List[modifiedsignal03],44100]]

最大絶対値の10% 未満を0にして圧縮

modified10=ModifiedListNeglectingRelativeSmallElement[signalWaveletData,0.10,X];

num=Count[modified10,X]

圧縮比

Print["Compression Raio for original volume= ", N[(Length[signal]-num)/Length[signal]]*100,"%"]

修正Waveletデータ

modifiedwavelet10=modified10/. X->0;

modifiedsignal10=FirstHaarInverseWaveletTransform[modifiedwavelet10];

Sound[SampledSoundList[List[modifiedsignal10],44100]]

最大絶対値の15% 未満を0にして圧縮

modified15=ModifiedListNeglectingRelativeSmallElement[signalWaveletData,0.15,X];

num=Count[modified15,X]

圧縮比

Print["Compression Raio for original volume= ", N[(Length[signal]-num)/Length[signal]]*100,"%"]

修正Waveletデータ

modifiedwavelet15=modified15/. X->0;

modifiedsignal15=FirstHaarInverseWaveletTransform[modifiedwavelet15];

Sound[SampledSoundList[List[modifiedsignal15],44100]]

最大絶対値の25% 未満を0にして圧縮

modified25=ModifiedListNeglectingRelativeSmallElement[signalWaveletData,0.25,X];

num=Count[modified25,X]

圧縮比

Print["Compression Raio for original volume= ", N[(Length[signal]-num)/Length[signal]]*100,"%"]

修正Waveletデータ

modifiedwavelet25=modified25/. X->0;

modifiedsignal25=FirstHaarInverseWaveletTransform[modifiedwavelet25];

Sound[SampledSoundList[List[modifiedsignal25],44100]]

最大絶対値の35% 未満を0にして圧縮

modified35=ModifiedListNeglectingRelativeSmallElement[signalWaveletData,0.35,X];

num=Count[modified35,X]

圧縮比

Print["Compression Raio for original volume= ", N[(Length[signal]-num)/Length[signal]]*100,"%"]

修正Waveletデータ

modifiedwavelet35=modified35/. X->0;

modifiedsignal35=FirstHaarInverseWaveletTransform[modifiedwavelet35];

Sound[SampledSoundList[List[modifiedsignal35],44100]]