Petitverse

ガウシアンフィルタとは

信号処理の世界ではFIRフィルタ、画像処理の世界では線形フィルタと呼ばれるフィルタの一種です(厳密には画像処理のそれは因果性というのが成り立たないので、信号処理のそれとは別物)

ガウシアンフィルタは以下の優れた特徴を持ちます
・理想的な直線位相特性を持つ(十分周波数が小さい信号を入力したとき、その形状が保持される)
・オーバーシュートやギブス現象が一切生じない
・レンズのピンぼけはガウシアンフィルタによって数学的に近似できる。
・全ての係数を有理数で近似することができる

以下の欠点もあります
・長さNのガウシアンフィルタではN/2の遅延を生じる
・周波数特性があまり急峻でない

奇数次のガウシアンフィルタは以下ような係数で近似することが出来ます
参考 ttps://t2r2.star.titech.ac.jp/rrws/file/CTT100586956/ATD100000413/
//TODO 偶数次のときの係数も求める

F[i]=(N-1)!/i!/(N-i-1)!/POW(2,N-1)


画像処理の場合は二次元で行うため、一次元の縦ベクトルと一次元の横ベクトルの行列積を計算します。

ただし、真面目に計算しようとすると値が大きくなりすぎるため、一般的には以下の漸化式を用いて計算します。

F[0]=1/POW(2,N-1)
F[i]=F[N-i-1]=F[i-1]*(N-i)/i
(i=0 to N>>1)

もう少し言うと、F[0]=1で計算すると整数で完結するので、F[0]=1で計算した後、ARYOPで1<<N-1を割ると誤差0で計算できます。


DEF GAUSS0 A[],N
　A[0]=1/(1<<N-1)
　FOR I=1 TO (N>>1)-1
　　A[i]=A[N-i-1]=A[i-1]*(N-i)/i
　NEXT
END


DIM A%[1024]
DEF GAUSS1 A[],N
　A%[0]=1
　FOR I=1 TO (N>>1)-1
　　A%[i]=A%[N-i-1]=A%[i-1]*(N-i)/i
　NEXT
　ARYOP 3,A,A%,1<<N-1
END

ガウシアンフィルタ…というか、FIRフィルタは畳み込みと呼ばれる計算のため、次数が大きい場合は高速フーリエ変換や高速数論変換を用いて高速に計算することができます。(画像処理の場合は次数が小さいことが一般的なので、普通に計算した方が速いことが多い)

入力をx,フィルタ係数をfとすると以下のように計算すると高速に計算できます。(プチコンの場合。他環境だとIFFTとFFTは逆なことが多い)
FFT(IFFT(x)*IFFT(f))

//TODO ソースコードをここに記入

ガウシアンフィルタのフーリエ変換,逆フーリエ変換
ガウシアンをフーリエ変換するとガウシアンになるという性質があります。
同様にFFT,IFFTの結果もガウシアンとして近似できるはずです。

http://zakii.la.coocan.jp/enumeration/79a_binomial_gaussian.htm

F=SQRT(2/PI()/(N+1))exp(-2/N*(k-N/2)^2)

…なんか色々やってたら、配列長をL,フィルタ長をNとしたら長さN'=L/2Nのガウシアンでフーリエ変換になるっぽい？

画像処理だとガウシアンフィルタ5x5とかですからね。
スマベだと二次元FFT/NTTがないのもあってARYOPの方が遥かに速いと思います。
信号処理だとフィルタ長が2〜3桁とかになるので、ARYOPはだいぶ無理がある。(それもプチコンだからってのは大きいですね。FORで回すより組み込みのIFFT/FFT使った方がマシって感じ。C言語とかだと全然FORで回す)

そしてFFTの畳み込みはクオリティ的にも微妙でした(誤差が大きかったから？あるいはIFFT→FFTでやるべきなのをFFT→IFFTでやったから？)
http://petitverse.hosiken.jp/community/petitcom/topic/?read=1821&ukey=1
ガウシアンフィルタのIFFTを直接求められれば速度的にもクオリティ的にも多少はマシになると思うんですけど、単なるぼかしでやる価値はないと思います。(まず二次元IFFTの時点で時間が厳しい)
…フィルタのサイズが3〜4桁とかになるなら話は別ですけど。(アナログ調の画像をリサイズするときとかは3〜4桁のフィルタ使うかな)

あった。IFFTの畳み込み
http://petitverse.hosiken.jp/community/petitcom/diary/?read=5657&ukey=0

とまぁ、そんなわけで、実際に1次元データに対してガウシアンフィルタを計算してみます。

ACLS
XSCREEN 1
OPTION STRICT
OPTION DEFINT

VAR N%=512,L%=63
DIM SIGNAL#[2048],FILTER#[2048] '信号長+フィルタ長は2048以内
DIM BLURED0#[2048],BLURED1#[2048] 'フィルタリング結果

DIM FFTR#[2048],FFTJ#[2048] 'FFT用バッファ
DIM CONVOLR#[2048],CONVOLJ#[2048] '畳み込み用バッファ
DIM T#[2048]

'ガウシアンフィルタの係数を計算(因果的でない)
DEF GAUSS A[],N
　VAR I%
　A[0]=POW(2,-N+1) 'N>31でオーバーフローするため修正
　FOR I%=1 TO (N>>1)-1
　　A[I%]=A[I%-1]*(N-I%)/I%
　A[N-I%-1]=A[I%]
　NEXT I%
END

'愚直に計算
DEF FIR0 O,IN,FIL,N,L
　VAR I%,J%
　VAR SIZE%=LEN(FIL)
　FOR I%=0 TO N-1
　　VAR TMP#=0
　　FOR J%=0 TO L-1
　　　TMP#=TMP#+IN[(I%-J%+SIZE%) MOD SIZE%]*FIL[(J%+SIZE%)MOD SIZE%]
　　NEXT J%
　　O[I%]=TMP#
　NEXT I%
END

'IFFTで計算
DEF FIR1 O,IN,FIL
　FILL FFTJ#,0:FILL CONVOLJ#,0 '虚部を0で初期化

　IFFT FFTR#,FFTJ#,IN,FFTJ#
　IFFT CONVOLR#,CONVOLJ#,FIL,CONVOLJ#

　'要素ごとに乗算
　COPY T#,FFTJ#
　ARYOP 2,FFTJ#,T#,CONVOLR#
　ARYOP 4,FFTJ#,FFTR#,CONVOLJ#,FFTJ#
　ARYOP 2,FFTR#,FFTR#,CONVOLR#
　ARYOP 2,T#,T#,CONVOLJ#
　ARYOP 1,FFTR#,FFTR#,T#

　FFT O,FFTJ#,FFTR#,FFTJ#
END


'入力,フィルタ初期化
FILL SIGNAL#,79,N%>>1,N%>>1
GAUSS FILTER#,L%

'フィルタ計算
FIR0 BLURED0#,SIGNAL#,FILTER#,N%,L%
FIR1 BLURED1#,SIGNAL#,FILTER#

VAR I%

'描画
FOR I%=1 TO 399
　GLINE I%-1,79-ROUND(SIGNAL#[I%-1]),I%,79-ROUND(SIGNAL#[I%])
　GLINE I%-1,159-ROUND(BLURED0#[I%-1]),I%,159-ROUND(BLURED0#[I%])
　GLINE I%-1,239-ROUND(BLURED1#[I%-1]),I%,239-ROUND(BLURED1#[I%])
NEXT I%

GPUTCHR 0,32,"入力",2,2
GPUTCHR 0,112,"出力(愚直)",2,2
GPUTCHR 0,192,"出力(IFFT)",2,2

WAIT 1E5