●伝達関数による周波数解析
OPTION ARITHMETIC COMPLEX
LET j=SQR(-1) !虚数単位
LET f=60 !周波数[Hz]
DEF w=2*PI*f !角周波数ω
SUB DispS(z) !複素数をS表示する ※スタインメッツ（Steinmetz）
PRINT ABS(z);
IF ABS(z)<>0 THEN
IF arg(z)<>0 THEN PRINT "∠";DEG(arg(z));"°";
END IF
PRINT
END SUB
!-------------------- ここまでがサブルーチン
!---------- ↓↓↓↓↓ ----------
LET xmax=6 !<----- ※要調整
LET ymin=-50
LET ymax=5
!●回路図 CRローパス・フィルタ
!vi・─R1┬─・vo
! C1
! │
! ≡
! 参考サイト http://sim.okawa-denshi.jp/CRlowkeisan.htm
LET R1=1.6e3 !1.6k[Ω]
LET C1=0.1e-6 !0.1μ[F]
DEF G(s)=1/(s*C1*R1+1) !入出力システムの伝達関数 G(s)=Vo/Vi=(1/(C1*R1))/(s+1/(C1*R1))
!---------- ↑↑↑↑↑ ----------
!!!SET bitmap SIZE 600,600 !画面を大きくする
SET WINDOW -0.5,xmax+0.5, ymin,ymax !表示領域
DRAW grid(1,5) !左端の目盛り
FOR f=1 TO xmax !ｘ軸が対数
PLOT TEXT ,AT f-0.3,-0.15: mid$("10 100 1k 10k 100k1M 10M 100M",4*(f-1)+1,4)
NEXT f
FOR xx=0 TO xmax STEP 0.025 !周波数[Hz]
LET f=10^xx !xx=LOG10(f)
LET t=ABS(G(j*w))
PLOT LINES: xx,20*LOG10(t); !利得[dB]
NEXT xx
PLOT LINES
SET TEXT COLOR 2
FOR k=ymax TO ymin STEP -5 !右端の縦軸目盛り
PLOT TEXT ,AT xmax,k: STR$(k*2)&"°" !※利得のグラフに合わせるために2倍する
NEXT k
SET LINE COLOR 2
FOR xx=0 TO xmax STEP 0.05 !周波数[Hz]
LET f=10^xx
LET th=arg(G(j*w))
IF th>0 THEN LET th=th-2*PI !0～-2πへ補正する <----- ※要調整
PLOT LINES: xx,DEG(th)/2; !位相θ[deg] ※利得のグラフに合わせるために1/2倍する
NEXT xx
PLOT LINES
END
●伝達関数による過渡解析
OPTION ARITHMETIC COMPLEX
DECLARE EXTERNAL SUB IFFT !逆フーリエ変換
DEF C(s)=G(s)*R(s) !応答関数
DEF R(s)=1/s !指令関数 ※ステップ関数
LET Vi=1 !1∠0°[V] ※仮の電圧源
!----- ↓↓↓↓↓ -----
LET T=2e-3 !時間区間 [0,T]
!●回路図 CRローパス・フィルタ
!vi・─R1┬─・vo
! C1
! │
! ≡
! 参考サイト http://sim.okawa-denshi.jp/CRlowkeisan.htm
LET R1=1.6e3 !1.6k[Ω]
LET C1=0.1e-6 !0.1μ[F]
DEF G(s)=1/(s*C1*R1+1) !入出力システムの伝達関数 G(s)=Vo/Vi=(1/(C1*R1))/(s+1/(C1*R1))
!----- ↑↑↑↑↑ -----
!逆ラプラス変換
! f(t)=lim [ω→∞] { 1/(2π*j)∫F(s)*Exp(s*t)ds [γ-j*ω,γ+j*ω] }、jは虚数単位、γ>0とωは実数。
!s=γ+j*2πωとすると、ds=j*2πdω。
!これを上式に代入して、変形すると、
! f(t)=Exp(γt)∫F(γ,ω)*Exp(j*2πω)dω [-∞,∞]
!右辺の積分部分は逆フーリエ変換。
LET N=1024*8 !データ総数 ※２のべき乗、大きいほど精度がよい
DIM d(0 TO N-1) !入出力用配列
LET Gamma=5 !※3～7
!γを大きくすれば精度は上がりそうだが、
!あとでExp(γt)をかけるので、tが大きいところで発散する。
LET rr=Gamma/T !γを決める
FOR k=0 TO N/2 !変換データの作成
LET Cs=C( COMPLEX(rr, 2*PI*k/T) ) !sk=γ+j*2π*k/T、k=0～n-1
LET Hw=(COS(2*PI*k/N)+1)/2 !データは離散なため、ハニング関数をかけて平滑化する ※精度の向上
LET d(k)=N/T*Cs*Hw !n/T*Y(k) * Hw 0～N/2
IF NOT(k=0 OR k=N/2) THEN LET d(N-k)=COMPLEX(Re(d(k)),-Im(d(k))) !N-1～N/2+1は、共役複素数
NEXT k
CALL IFFT(d) !高速逆フーリエ変換
SET WINDOW -T/8,T,-Vi*1.2,Vi*1.2 !表示領域
DRAW grid(T/4,Vi/5)
PLOT TEXT ,AT T*7/8,0: "[秒]"
LET dt=T/N !時間刻み幅 ⊿t
FOR k=0 TO N-1 STEP 8 !結果の表示 [0,T]
LET tt=k*dt !経過時間
PLOT LINES: tt, Re(d(k))*EXP(rr*tt); !Exp(γt)をかける ※Exp(γ*k*T/n)、k=0～n-1
PRINT Re(d(k))*EXP(rr*tt) !debug
NEXT k
END
EXTERNAL SUB IFFT(x()) !高速逆フーリエ変換 x() : 入力/出力データ
OPTION ARITHMETIC COMPLEX
DECLARE EXTERNAL SUB FFTMAIN
LET nx=SIZE(x)
LET theta=2*PI/nx ! W = Exp(-j * 2π/N * -1) = Exp(j * theta)とする
CALL FFTMAIN(x, theta)
MAT x=(1/nx)*x
END SUB
EXTERNAL SUB FFTMAIN(x(), theta)
OPTION ARITHMETIC COMPLEX
LET nx=SIZE(x)
IF MOD(nx, 2)<>0 THEN !DFTの計算
DIM w(0 TO nx-1), xtmp(0 TO nx-1)
MAT xtmp=x
FOR k=0 TO nx-1
LET tmp=theta*k
FOR n=0 TO nx-1
LET w(n)=EXP( COMPLEX(0, tmp*n) )
NEXT n
LET x(k)=DOT(w, xtmp)
NEXT k
ELSE !２分して再帰呼出し
LET hnx=nx/2
DIM x0(0 TO hnx-1), x1(0 TO hnx-1)
FOR k=0 TO hnx-1
LET x0(k)=x(k)+x(k+hnx)
LET wk=EXP( COMPLEX(0, theta*k) )
LET x1(k)=wk*(x(k)-x(k+hnx))
NEXT k
CALL FFTMAIN(x0, 2*theta)
CALL FFTMAIN(x1, 2*theta)
FOR k=0 TO hnx-1
LET x(2*k)=x0(k)
LET x(2*k+1)=x1(k)
NEXT k
END IF
END SUB
投稿日：2009年 2月17日(火)12時23分48秒