■６．１　波形の合成

表6-1 波形合成用共通ライブラリ（waveSynthesis.py）

from graph import * #graph.pyについては表4-1参照
import math
def dspV(canvas, R,CL,width=3):#グラフ化
    N=len(R)
    for i in range(1, N):
        X1=R[i-1][0]; Y1=R[i-1][1]; X2=R[i][0]; Y2=R[i][1]
        if (Y1>=-2 and Y1<=2)and(Y2>=-2 and Y2<=2):
            drawV(canvas, X1,Y1,X2,Y2,CL,
                  X0=240,Y0=140,width=width,Xscale=100,Yscale=50)
def exeFun(fn, T,V0,Phi):#sin関数設定(T:周期, V0:振幅, Phi:位相)
    NST=-2;NED=2 #サンプリング開始と終了時刻
    M=200; IST=NST*M;IED=NED*M;DX=1/M # M:グラフ上の分割数
    X1=IST*DX; Y1=V0*fn(X1*2*math.pi/T+Phi); R=[]
    for i in range(IST,IED+1):
        R.append([X1,Y1])
        X2=i*DX; Y2=V0*fn(X2*2*math.pi/T+Phi); X1=X2;Y1=Y2
    R.append([X1,Y1]); return R
def setSin(T,V0,Phi):#sin関数設定(T:周期, V0:振幅, Phi:位相)
    return exeFun(math.sin,T,V0,Phi)
def setCos(T,V0,Phi):#cos設定(T:周期, V0:振幅, Phi:位相)
    return exeFun(math.cos,T,V0,Phi)
def setTan(T,V0,Phi):#tan設定(T:周期, V0:振幅, Phi:位相)
    return exeFun(math.tan,T,V0,Phi)
def initScale(canvas):#初期スケールの表示
    drawScale(canvas, 40,40, 15, 10, 10, 5, -2, -2,
          0.1, 0.1, 5,10, 5, 5, 40, 40,"%.1f","%.1f")
def typeCheck(A):# int型かfloat型のチェック
    if type(A) is int or type(A) is float: return True
    return False
def mltConst(V,A):# 波形×定数
    NN=len(A);R=[]
    for i in range(NN):R.append([A[i][0],V*A[i][1]])
    return R
def mltWave(A1,A2):# 波形×波形
    if typeCheck(A1): return mltConst(A1,A2)
    if typeCheck(A2): return mltConst(A2,A1)
    N1=len(A1); N2=len(A2); NN=N1;R=[]
    if NN>N2: NN=N2# 乗算では短いほうを基準とする
    for i in range(NN):R.append([A1[i][0],A1[i][1]*A2[i][1]])
    return R
def addConst(V,A):# 定数×波形
    NN=len(A);R=[]
    for i in range(NN):R.append([A[i][0],V+A[i][1]])
    return R
def addWave(A1,A2):# 波形＋波形
    if typeCheck(A1): return addConst(A1,A2)
    if typeCheck(A2): return addConst(A2,A1)
    N1=len(A1); N2=len(A2); NN=N1;R=[]
    if NN<N2: NN=N2#加算では長いほうを基準とする
    for i in range(NN):
        V=0;X=0    #不足する側を0とする
        if i<N1: V+=A1[i][1];X=A1[i][0]
        if i<N2: V+=A2[i][1];X=A2[i][0]
        R.append([X,V])
    return R
def subConst1(V,A):# 波形－定数
    NN=len(A);R=[]
    for i in range(NN):R.append([A[i][0],V-A[i][1]])
    return R
def subConst2(A,V):# 定数－波形
    NN=len(A);R=[]
    for i in range(NN):R.append([A[i][0],A[i][1]-V])
    return R
def subWave(A1,A2):# 波形－波形
    if typeCheck(A1): return subConst1(A1,A2)
    if typeCheck(A2): return subConst2(A1,A2)
    N1=len(A1); N2=len(A2); NN=N1;R=[]
    if NN<N2: NN=N2 #減算では長いほうを基準とする
    for i in range(NN):
        V=0;X=0     #不足する側を0とする
        if i<N1: V+=A1[i][1];X=A1[i][0]
        if i<N2: V-=A2[i][1];X=A2[i][0]
        R.append([X,V])
    return R
def setErrDT(SGN,A):
    if A>=0: return SGN*1E30
    else: return -SGN*1E30
def divConst1(V,A):# 波形／定数
    NN=len(A);R=[];SGN=1
    if V<0: SGN=-1
    for i in range(NN):
        if abs(A[i][1])<1E-30:Y=setErrDT(SGN,A[i][1])
        else: Y=V/A[i][1]
        R.append([A[i][0],Y])
    return R
def divConst2(A,V):# 定数／波形
    NN=len(A);R=[]; CHK=False
    if abs(V)<1E-30:CHK=True
    for i in range(NN):
        if CHK:Y=setErrDT(1,A[i][1])
        else: Y=A[i][1]/V
        R.append([A[i][0],Y])
    return R
def divWave(A1,A2):# 波形／波形
    if typeCheck(A1): return divConst1(A1,A2)
    if typeCheck(A2): return divConst2(A1,A2)
    #print(A1,A2)
    N1=len(A1); N2=len(A2); NN=N1;R=[]
    if NN>N2: NN=N2 #除算では短いほうを基準とする
    for i in range(NN):
        if abs(A2[i][1])<1E-30:Y=setErrDT(1,A[i][1])
        else: Y=A1[i][1]/A2[i][1]
        R.append([A1[i][0],Y])
    return R
def minusWave(A):#－波形
    NN=len(A); R=[]
    for i in range(NN): R.append([A[i][0],-A[i][1]])
    return R
def dspGraphProc(canvas,ft):
    N=len(ft)
    for i in range(N):
        dspV(canvas,exeFun(ft[i][0],ft[i][1],
                           ft[i][2],ft[i][3]),ft[i][4])
def dspGraphInit():
    root=initTk("Numerical Analysis")
    canvas=initCanvas(root,500,300)
    initScale(canvas);return canvas
def dspGraph(ft):dspGraphProc(dspGraphInit(),ft)

表6-2 波形合成用共通ライブラリ（表6-1）を使った定義

from waveSynthesis import * #waveSynthesis.pyについては表6-1参照
def sqrSig(t):
    s=0
    for i in range(1,401,2):s+=math.sin(i*t)/i
    return s*4/math.pi
dspGraph([[math.sin,2,1,0,"blue"  ],
          [math.cos,2,1,0,"red"   ],
          [math.tan,2,1,0,"green" ],
          [sqrSig  ,2,1,0,"purple"]])

表6-3 波形合成確認のためのプログラム例

from waveSynthesis import * #waveSynthesis.pyについては表6-1参照
canvas=dspGraphInit()
#          周期，振幅，位相
A = setSin( 2.0, 1.0,   0      )#波形定義
B = setSin( 2.0, 0.5, math.pi/2)
C = addWave(A,B)                #波形合成
dspV(canvas, A, "green")        #信号表示
dspV(canvas, B, "blue")
dspV(canvas, C, "red")

表6-4 三角級数合成波形の確認用プログラム例

from waveSynthesis import * #waveSynthesis.pyについては表6-1参照
def triSignal(t):
    global N; NN=N*2+1; s = 0
    for i in range(1,NN,2): s+=math.cos(i*t)/(i*i)
    return s
N=1; canvas=dspGraphInit()
R=exeFun(triSignal,2,1,0)
dspV(canvas, R, "red" )

■いろいろなフーリエ級数

A. 矩形波

def sqrSignal(t):
   global N; NN=N*2+1; s = 0
   for i in range(1,NN,2): s += math.sin(i*t)/i
   return s*4/math.pi

B. のこぎり波

def sawSignal(t):
   s = 0;SG=1; global N
   for i in range(1,N+1):
     s +=SG* math.sin(i*t)/i; SG=-SG
   return 2*s/math.pi

C. 二次式（ t 2 ： －π < t < π）

def quadratic(t):#表示上の都合で1/5を返却
    s = 0; SG=1; global N
    for i in range(1, N+1):
        s += SG*math.cos(i*t)/(i*i); SG=-SG
    return (math.pi*math.pi/3 - 4*s)/5

D. 一次式と0の繰り返し

def repeatXand0(t):#表示上の都合で2/πを返却
   s = 0;SG=1 ; global N
   for i  in range(1,N+1):
     s += 2*math.cos((2*i-1)*t)/((2*i-1)*(2*i-1))/math.pi
     s -= SG* math.sin(i*t)/i; SG=-SG  
   return (math.pi / 4 - s)*2/math.pi

E. 中途半端な余弦関数

def evenSignal(t): 
    s = 0; SG=1; global N
    for i in range(1,N+1):s += SG*math.cos(i*t)/(2-i*i); SG=-SG
    return s

F. 三次式（－π ≦ t ＜ π）

def cubicSignal(t):#表示上の都合で1/(π^3)を返却
    s = 0;SG=1; Pi=math.pi; global N
    for i in range(1,N+1):
        s += (2*Pi*Pi/i-12/(i*i*i))*SG*math.sin(i*t); SG=-SG
    return s/(Pi*Pi*Pi)   

G. 二段溝形

def step2Signal(t):
    s = 0; SG=-1;global N
    for i in range(1,N+1):
        ii=2*i; ii2=ii-1;SG=-SG
        s -= 2 * SG*math.cos((ii2)*t)/(ii2)
        s -= (1 +SG)*math.sin(ii*t)/i
    return s/math.pi

■６．２　フーリエ級数（プログラム例なし）

■６．３ 複素フーリエ級数

表6-5 複素フーリエ級数用グラフ化

from waveSynthesis import * #waveSynthesis.pyについては表6-1参照
def expJ(TH):#オイラーの公式
    return complex(math.cos(TH),math.sin(TH))
def triComplex(t):#ここに関数を書く
    s = complex(0, 0); SG=1
    for i in range(-200, 201):
        ii=2*i-1; s+=expJ(ii*t)/(ii*ii)
    s=s*4/(math.pi*math.pi)
    return s.real
canvas=dspGraphInit()
R=exeFun(triComplex,2,1,0)#関数名を指定する
dspV(canvas, R, "red" )

A. のこぎり波

def sawComplex(t):
    s = complex(0, 0); SG=1
    for i in range(-400, 401):
        if i != 0: s += SG*expJ(t*i)/i
        SG=-SG
    return -s.imag/math.pi # jを乗じる替わりに虚部の符号逆転

B. 三角波（表6-5参照）

C. 矩形波（方形波）

def sqrComplex(t):
    s = complex(0, 0)
    for i in range(-400,401):
        ii=2*i-1;s += expJ(ii * t)/ ii
    return 0.5+ s.imag/math.pi #jを乗じる替わりに虚部の符号逆転

D. 階段状パルス

def stepComplex(t):#表示の都合で1/2の値を返却
    s = complex(0, 0);Pi=math.pi
    for i in range(-400,401):
        ii = i % 4
        if i == 0  : ex=complex(1.5, 0)
        elif ii==0 : ex=complex(0, 0)
        elif ii==1 : ex=complex(-2/(i*Pi), 0)
        elif ii==2 : ex=complex(0,-2/(i*Pi))
        else       : ex=complex(2/(i*Pi), 0)
        s += ex*expJ(i*t)
    return s.real/2

E. 複流RZ(Return to Zero)パルス

def doubleflowRZ(t):
    s = complex(0, 0)
    for i in range(-400,401):
        if i % 2:
            x = i % 4
            if x == 1 : ex = complex(-2, -2)
            else:       ex = complex(-2, 2)
            s += ex*expJ(i * t)/i
    return -s.imag/(2*math.pi)#jを乗じる替わりに虚部の符号逆転

F. 二次式（ t 2 ）

def quadraticComplex(t):
    s = complex(0, 0); SG=1
    for i in range(-100,101):
        if i !=0 : 
s += SG*expJ(i*t)/(i*i)
        SG=-SG
    return 1/3+2*s.real/(math.pi*math.pi)

G. 正弦関数の絶対値（ | sin πt | ）

def absSinComplex(t):
    s = complex(0, 0)
    for i in range(-100,101):
        s += expJ(i*t)/(1 - 4*i*i)
    return s.real*2/math.pi
 
H. のこぎり波（２）

def saw2Complex(t):
    s = complex(0, 0);SG=1; Pi=math.pi
    for i in range(-300,301):
        if i!=0:
            IP=i*Pi;IP2=2*IP;IP22=2*IP*IP
            ex=complex((SG-1)/IP22, SG/IP2)
            s += ex*expJ(i*t)
        SG=-SG
    return s.real + 0.25

■６．４　フーリエ変換

表6-6 フーリエ変換確認用プログラム

from waveSynthesis import * #waveSynthesis.pyについては表6-1参照
def expJ(TH):#オイラーの公式
    return complex(math.cos(TH),math.sin(TH))
def integralFun(fn, T,V0,Phi):#関数設定(T:周期, V0:振幅, Phi:位相)
    NST=-2;NED=2 #サンプリング開始と終了時刻
    M=200; IST=NST*M;IED=NED*M;DX=1/M # M:グラフ上の分割数
    X1=IST*DX; Y1=V0*inverseTran(fn,X1*2*math.pi/T+Phi); R=[]
    for i in range(IST,IED+1):
        R.append([X1,Y1])
        X2=i*DX; Y2=V0*inverseTran(fn,X2*2*math.pi/T+Phi)
        X1=X2;Y1=Y2
    R.append([X1,Y1]); return R
dw = 0.05
def inverseTran(fn, t):
    global dw; w = 0; u2 = 0; s = complex(0, 0)
    for i in range(800):
        w +=dw; u1 = u2; u2 = fn(w)# 関数値の計算
        udt = (u1 + u2) * dw / 2   # 台形の公式(1区間分)
        cudt=udt.conjugate()       # 共役複素数
        s += expJ( w * t) * udt    # プラス側
        s += expJ(-w * t) * cudt   # マイナス側
    return s.real /(2*math.pi)
def delta(w):# 理論的には無限大だが，数値積分上の幅があるので
    # 階段関数の数値微分上の値 1 / dwを返却する
    global dw
if abs(w)<0.5*dw : return 1 / dw
    return 0
def expF(w):#ここに関数を書く
    return complex(2/(1+w*w),0)
canvas=dspGraphInit()
R=integralFun(expF,2,1,0)#関数名を指定する
dspV(canvas, R, "red" )

A. 矩形波（方形波）

def sqrF(w):
    if abs(w)<1E-30: return complex(0,0)
    return complex(2*math.sin(w)/w,0)

B. ディラックのデルタ関数

def deltaF(w):# 本来は1を返却して
# δ関数になることを確認したいが，
# 合計すると大きすぎる値になるの
# で，以下の値を返却する。
    global dw
    return dw*math.pi/2

C. 余弦波信号

def cosF(w): 
A=math.pi*(delta(w-1))
B=math.pi*(delta(w+1))
return complex(A+B,0)
　　　　　　　　　

D. 正弦波信号

def sinF(w):
A=math.pi*(delta(w-1))
B=math.pi*(delta(w+1))
return complex(0,B-A)
 
D. 指数関数（その１）

def exp2F(w):
    return complex(1,-w)/(1+w*w)
 
E. 二次関数

def quadF(w):
#横軸の値 = 1/π≒±0.31831
    w2 = w * w; w3 = w2 * w
    A1 =-4 * math.cos(w)/w2
    A2 = 4 * math.sin(w)/w3
    return complex(A1+A2,0)   

E. 指数関数（その２）

def expF(w):
    return complex(2/(1+w*w),0)

 
【補足】ディラックのデルタ関数について

from waveSynthesis import * #waveSynthesis.pyについては表6-1参照
def expJ(TH):#オイラーの公式
    return complex(math.cos(TH),math.sin(TH))
def diracDeltaFun(width, T,V0,Phi):#関数設定(T:周期, V0:振幅, Phi:位相)
    NST=-2;NED=2 #サンプリング開始と終了時刻
    M=200; IST=NST*M;IED=NED*M;DX=1/M # M:グラフ上の分割数
    X1=IST*DX; Y1=V0*diracDelta(width, X1*2*math.pi/T+Phi); R=[]
    for i in range(IST,IED+1):
        R.append([X1,Y1])
        X2=i*DX; Y2=V0*diracDelta(width, X2*2*math.pi/T+Phi)
        X1=X2;Y1=Y2
    R.append([X1,Y1]); return R
def diracDelta(width, t):
    N=400; w = 0; dw =0.05; n=int(width/ dw); s = complex(0, 0)
    for i in range(n):
        w +=dw ; w2 = w *(t-0.5*math.pi)
        s += expJ(w2)*dw; s += expJ(-w2)*dw
    return s.real*10/n # nに比例して大きくなるのでnで除す
canvas=dspGraphInit()
# 最初の引数で積分範囲を指定する。
#
#               積分範囲,周期,振幅,位相
#
R=diracDeltaFun(   4    , 2  ,  1 , 0  )
dspV(canvas, R, "red" )