新しく発言する EXIT インデックスへ
FFT
　　FFT takufuji 2005/07/16 01:01:26 

　　FFT takufuji 2005/07/16 01:01:26  ツリーへ

FFT 返事を書く
takufuji 2005/07/16 01:01:26
　白石先生、十進BASICを活用している皆さん、はじめまして。
白石先生にはただ頭が下がるばかりです。大変僭越ですがFFTを作ったのでよかったら使ってください。
------ 以下FFT -------
OPEN #1: NAME "D:\いろいろ\数学、物理\BASICw32\Onsa_440.txt"
OPEN #2: NAME "D:\いろいろ\数学、物理\BASICw32\Onsa_440_FFT2.TXT"
OPTION ARITHMETIC COMPLEX !複素関数
OPTION BASE 0 !配列は0から
DIM f(10000),c(100,10000),w(10000),n(10000) !ポイント数に注意
LET l=4096 !ポイント数(2の乗数)
LET fs=22050 !サンプリング周波数[Hz]
LET m=l
LET t0=TIME
!--------------------- データ入力 --------------------------
FOR i=0 TO l-1
INPUT #1: f(i)
NEXT i
ERASE #2
!---------------------------------------------------------------
!----------------- nのビットリバースとωの作成 ---------------
FOR k=0 TO l-1
LET c(0,k)=complex(f(k),0)
LET w(k)=complex(COS(2*(k)*PI/l),-SIN(2*(k)*PI/l))
FOR j=0 TO LOG2(l)-1
LET n(k)=n(k)+MOD(INT(k/2^j),2)*2^(LOG2(l)-1-j)
NEXT j
NEXT k
!---------------------------------------------------------------
!---------------------- バタフライ演算 ----------------------
FOR j=0 TO LOG2(l)-1
LET kk=2^j-1
LET m=m/2
FOR jj=0 TO kk
FOR k=0 TO m-1
LET c(j+1,k+2*jj*m)=c(j,k+2*jj*m)+c(j,k+2*jj*m+m)*w(n(2*jj))
LET c(j+1,k+2*jj*m+m)=c(j,k+2*jj*m)+c(j,k+2*jj*m+m)*w(n(2*jj+1))
NEXT k
NEXT jj
NEXT j
!---------------------------------------------------------------
!----------------- スペクトルのビットリバース ----------------
FOR k=0 TO l-1
LET c(LOG2(l)+1,n(k))=c(LOG2(l),k)
NEXT k
!---------------------------------------------------------------
!--------------------- スペクトルの出力 ----------------------
SET WINDOW 0,fs/2,0,10
FOR k=0 TO l/2-1
PRINT #2: Re((2/l)*c(LOG2(l)+1,k));",";(2/l)*Im(c(LOG2(l)+1,k))
! PRINT k*fs/l;(2/l)*SQR(Re(c(LOG2(l)+1,k))^2+Im(c(LOG2(l)+1,k))^2)
PLOT LINES: k*fs/l,(2/l)*SQR(Re(c(LOG2(l)+1,k))^2+Im(c(LOG2(l)+1,k))^2);
NEXT k
!---------------------------------------------------------------
CLOSE #1
CLOSE #2
PRINT time-t0;"秒かかりました"
END

インデックスへ EXIT
新規発言を反映させるにはブラウザの更新ボタンを押してください。