Petitverse

'FFT/IFFTの要素数が引数配列の要素数から自動決定するのやめてほしい。さすがに扱いづらすぎる


'グロヴァ(ローカル変数だとメモリが断片化するので)
'要素数はてきとー(大きすぎると速度が著しく低下するので、大きすぎないサイズ)
DIM AR#[16384],BR#[16384]
DIM AJ#[16384],BJ#[16384]
DIM T#[16384]

'畳み込みのプロ・畳小宮の存在だ
'AとBを畳み込んだ結果をOに返す
DEF CONVOL O[],A[],B[]
　FILL AJ#,0:FILL BJ#,0
　'A,Bの要素数が16384じゃないとき用
　FILL AR#,0,LEN(A)-1
　FILL BR#,0,LEN(B)-1
　COPY AR#,A:COPY BR#,B
　IFFT AR#,AJ#,AR#,AJ#
　IFFT BR#,BJ#,BR#,BJ#
　'T=AR*BJ+BR*AJ
　ARYOP 2,T#,BR#,AJ#
　ARYOP 4,T#,AR#,BJ#,T#
　'AJ=AR*BR-AJ*BJ
　ARYOP 2,BJ#,AJ#,BJ#
　ARYOP 2,AJ#,AR#,BR#
　ARYOP 1,AJ#,AJ#,BJ#
　FFT AJ#,T#,AJ#,T#
　COPY O,AJ#,0,LEN(O)
END

'INのパワースペクトルをOに返す
'TODO なんかバグってそうなので修正する
DEF FFTPOW O[],IN[]
　FILL AJ#,0
　FILL AR#,0,LEN(IN)-1
　COPY AR#,IN
　FFT T#,AJ#,AR#,AJ#
　ARYOP 2,T#,T#,T#
　ARYOP 4,AJ#,AJ#,AJ#,T#
　COPY O,AJ#,0,LEN(O)
END

'INの自己相関関数をOに返す
DEF AUTOC O,IN
　'TODO 書き方を調べておく
END

'要素数Nのハミング窓
'FFTWFN A,1の要素数制限なし版
DEF HAMMING A,N
　VAR I%
　FOR I%=0 TO N-1
　　A[I%]=0.54-COS(2*PI()*(I%)/N)*0.46
　NEXT I%
END

'要素数Nで正規化カットオフ周波数Fの因果的なLPF
DEF MKLPF A,N,F
　HIMMING T#,N '今回はハミング窓
　VAR I%
　FOR I%=0 TO N-1
　　VAR J%=I%-(N>>1)
　　IF J% THEN
　　　A[I%]=SIN(2*PI()*J%*F)*B#{I%]/J%/PI()
　　ELSE 'J%=0ではゼロ除算が生じるので、極限値を入れる
　　　A[I%]=2*F%*B#[I%]
　　ENDIF
　NEXT I%
END