『FFT によるノイズ除去』

Maple 9 で利用可能になった『DiscreteTransform』(離散データ用変換)パッケージで提供されている FFT 関数を利用して、データのノイズ除去を行うサンプルです。
数式処理を用いれば、窓関数の構築も直観的な定義に基づいて行うことができるようになります。

>   

>    restart;

>   

 サンプル・データ


このワークシートでは、次のような1つのピークを持ったデータを用いることにします。

>    peakFunc := unapply(eval(a+b/(c+(t-d)^2), {a=10, b=40000, c=380, d=128}), t);

peakFunc := proc (t) options operator, arrow; 10+40000/(380+(t-128)^2) end proc

>    noisyData := [seq(peakFunc(i)*(0.6+0.8*evalf(rand()/10^12)), i=1..256)]:

なお、データには乱数を用いて適当なノイズを与えています。

>    with(plots):
with(plottools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined


元の関数は赤線で描画しています。

>    grFunc := plot(peakFunc(t), t=1..256, color=red, thickness=2):
grData := listplot(noisyData, color=khaki):
display([grFunc, grData]);

[Maple Plot]

>   

>   

 フーリエ変換による成分分解


与えられたデータ集合(信号データ)に対して、各周波数成分に分解するために高速フーリエ変換を施します。
Maple 8 以前の離散フーリエ変換では、データ長が2のベキ乗に限られていましたが、Maple 9 で提供された DiscreteTransform パッケージの FFT 関数から任意の長さ(次元)のデータに対する高速フーリエ変換が可能となっています。

>    with(DiscreteTransforms);

[FourierTransform, InverseFourierTransform]

まず、フーリエ変換の関数を利用するために、データを Array 型へと変更します。あわせて、実数成分と虚数成分とをそれぞれデータとして定義しておきます。

>    dataReal := Array([seq(j, j=noisyData)]);

dataReal := Array(%id = 12667420)

>    dataImg := Array(1..256, 0);

dataImg := Array(%id = 12630108)


このデータに対して、高速フーリエ変換を施します。

>    fourierReal, fourierImg := FourierTransform(dataReal, dataImg);

fourierReal, fourierImg := Array(%id = 12748280), Array(%id = 35572876)

>   

 窓関数と逆フーリエ変換


変換されたデータ fourierReal, fourierImg に対して、フィルタ関数を適用しましょう。数学的な表現をそのまま記述できる Maple では、窓関数を簡単に用意することができます。
まず、ヘヴィサイド関数を定義しておきます。(注1参照)

>    heaviside := proc (x) options operator, arrow; piecewise(0 <= x,1,x < 0,0) end proc

heaviside := proc (x) options operator, arrow; piecewise(0 <= x,1,x < 0,0) end proc

続いて、いくつかのフィルタ関数を用意します。

>    N := 20:

>    # Rectanglar window
filter1 := x -> evalf(heaviside(x)*heaviside(N-x)):
plot(filter1(x), x=0..2*N);

[Maple Plot]

>    filter2 := x -> evalf(heaviside(x)*heaviside(N-x)*(1-x/N)):
plot(filter2(x), x=0..2*N);

[Maple Plot]

>    # Hanning window
filter3 := x -> evalf(heaviside(x)*heaviside(N-x)*(0.5-0.5*cos(2*Pi*x/N))):
plot(filter3(x), x=0..2*N);

[Maple Plot]

>    # Hamming window
filter4 := x -> evalf(heaviside(x)*heaviside(N-x)*(0.54-0.46*cos(2*Pi*x/N))):
plot(filter4(x), x=0..2*N);

[Maple Plot]

これらフィルタを適用したデータを作ります。あらかじめ、フィルタ処理をひとつの関数としてまとめておきます。

>    applyFilter := proc(fil)
  local filReal, filImg;

  filReal := Array( [seq(fourierReal[i]*fil(i), i=1..256)] ):
  filImg := Array( [seq(fourierImg[i]*fil(i), i=1..256)] ):
  return( filReal, filImg );
end proc:

>    filterdReal, filterdImg := applyFilter(filter1);

filterdReal, filterdImg := Array(%id = 12649028), Array(%id = 12786940)

フィルタ処理された信号データを、逆フーリエ変換します。
(実部のみ抜き出します)

>    denoisedData := InverseFourierTransform(filterdReal, filterdImg)[1];

denoisedData := Array(%id = 12835888)


グラフを作成し、元の信号データの振る舞いと比較してみましょう。

>    grDenoised := listplot(denoisedData, color=green, thickness=2):
display([grData, grDenoised], title="Denosied data", color=grey);

[Maple Plot]

今回定義したすべてのフィルタを用いた結果を見てみましょう。

>    makeDenoising := proc(fil, grTitle::string)
  local gr, denoised;

  denoised := InverseFourierTransform(applyFilter(fil))[1];
  gr := listplot(denoised, color=green, thickness=2);
  display([grData, gr], title=grTitle, tickmarks=[2,0]);
end proc:

>    gr := array(1..2, 1..2):

>    gr[1,1] := makeDenoising(filter1, "filter1"):
gr[1,2] := makeDenoising(filter2, "filter2"):
gr[2,1] := makeDenoising(filter3, "filter3"):
gr[2,2] := makeDenoising(filter4, "filter4"):

>    display(gr, scaling=constrained);

[Maple Plot]

>   

閉じる