■７．１　離散時間信号（プログラム例なし）

■７．２　離散時間フーリエ変換

表7-1  DTFT用共通関数（ファイル名：DTFT.py）

from graph import * #graph.pyについては表4-1参照
def initScaleComplex(title="Numerical Analysis"):#初期スケールの表示
    root=initTk(title)
    canvas=initCanvas(root,500,300)
    drawScale(canvas, 40,40, 15, 10, 10, 2, -1, -3,
          0.1, 0.1, 5,10, 5, 10, 20, 60,"%.1f","%.1f")
    canvas.create_text(60,30,text="実部",font=('Courier',10))
    drawScale(canvas, 280,40, 15, 10, 10, 2, -1, -3,
          0.1, 0.1, 5,10, 5, 10, 20, 60,"%.1f","%.1f")
    canvas.create_text(300,30,text="虚部",font=('Courier',10))
    return canvas
def dspComplex(canvas, R,XSC,YSC,CL,width=3):#複素数のグラフ化
    N=len(R)
    for i in range(1, N):
        X1=R[i-1][0]*XSC;X2=R[i][0]*XSC
        if (X1>=-1 and X1<=1)or(X2>=-0.99 and X2<=1.01):
            Y1=R[i-1][1].real*YSC;Y2=R[i][1].real*YSC
            if (Y1>=-3 and Y1<=3)and(Y2>=-3 and Y2<=3):
                drawV(canvas, X1,Y1,X2,Y2,CL,X0=140,Y0=100,
                      width=width,Xscale=100,Yscale=20)
            Y1=R[i-1][1].imag*YSC; Y2=R[i][1].imag*YSC
            if (Y1>=-3 and Y1<=3)and(Y2>=-3 and Y2<=3):
                drawV(canvas, X1,Y1,X2,Y2,CL,X0=380,Y0=100,
                      width=width,Xscale=100,Yscale=20)
def expJ(TH):#オイラーの公式
    return complex(math.cos(TH),math.sin(TH))
def DTFT(L,n0,w):#DTFT(離散信号をリストで指定)
    NUM=len(L); s=complex(0,0)
    for N in range(0,NUM):
        s+= expJ(-w * (N-n0))*L[N]
    return s
def fnDTFT(fn, w0,w,NUM):#永久信号用（関数指定）
    s=complex(0,0);n0=int(NUM/2);dw=1.0/n0
    for N in range(0,NUM+1):
        s+= expJ(-w * (N-n0))*fn((N-n0) * w0)
        s+= expJ( w * (N-n0))*fn(-(N-n0) * w0)
    return s/2
def listDTFT(L,n0):# -π～πまでDTFT
    DP=math.pi*2/2000; RL=[]
    for i in range(2001):
        TH=(i-1000)*DP; X=DTFT(L,n0,TH)
        RL.append([TH/math.pi,X])
    return RL
def listFnDTFT(fn,w0,NUM):# -π～πまでDTFT(関数指定)
    DP=math.pi*2/NUM; RL=[];n0=NUM/2
    for i in range(NUM+1):
        TH=(i-n0)*DP; X=fnDTFT(fn,w0,TH,NUM)
        RL.append([TH/math.pi,X])
        #print(TH/math.pi,X)
    return RL
def IDTFT(L,n):#逆DTFT
    NUM=len(L); s=complex(0,0); 
    for i in range(1,NUM):
        # 定義上はdwにπを乗ずべきだが計算量
        # を減らすため，ここではπを乗ぜず，
        # また，最後にπで除さないこととする。
        dw = L[i][0]-L[i-1][0]
        V1 = L[i-1][1]*expJ(L[i-1][0]*n)
        V2 = L[i  ][1]*expJ(L[i  ][0]*n)
        s += (V1+V2)*dw/2
    return s /2
def listIDTFT(L):# -π～πまでIDTFT
    NUM=len(L); R=[]; NN=int((NUM-1)/2)
    for i in range(-NN,NN+1):
        R.append([i,IDTFT(L,i)])
    return R
def adjustUnit(L):# 0からNを-π～πに単位調整
    NUM=len(L); NN=(NUM-1)/2; dt=1/NN; R=[]
    for i in range(NUM):
        R.append([L[i][0]*dt, L[i][1]])
    return R

（６）離散時間フーリエ変換の実際

①単位インパルス信号

from DTFT import * #DFT.pyについては表7-1参照（以下同じ）
canvas=initScaleComplex()
L = [0,1,0]; # 中央を0番目とする。リストの外側は0とする。
R = listDTFT(L,(len(L)-1)/2); dspComplex(canvas, R, 1,1,"red" )
IR= adjustUnit(listIDTFT(R)); dspComplex(canvas, IR,1,1,"blue" )

②余弦波信号

from DTFT import *
canvas=initScaleComplex()
m = 50 #この値を変更する。
R=listFnDTFT(math.cos, math.pi/2, m)
dspComplex(canvas, R, 1, 4/m,"red" )# Y方向スケール2倍(本来は2/m)

（分割数1000）
from DTFT import * 
canvas=initScaleComplex()
m=1000; R=listFnDTFT(math.cos, math.pi/2,m)
dspComplex(canvas, R, 1,2/m,"red" )
IR=adjustUnit(listIDTFT(R))
dspComplex(canvas, IR,math.pi*12.5, 1,"blue" )

正弦波（余弦波の場合の第3行目を以下のように変更）
m=1000; R=listFnDTFT(math.sin, math.pi/2,m)
	
合成波（分割数1,000）

from DTFT import * 
def funSynth(w0):
    return math.sin(w0/2)+0.8*math.sin(w0/3)+math.cos(w0/2)
canvas=initScaleComplex()
m=1000; R=listFnDTFT(funSynth, math.pi,m)
dspComplex(canvas, R, 1, 2/m,"red" )
IR=adjustUnit(listIDTFT(R))
dspComplex(canvas, IR,math.pi*5, 1,"blue" )

④パルス信号

from DTFT import * 
canvas=initScaleComplex()
L=[1,1,1]# 中央を0番目とする。リストの外側は0とする。
R=listDTFT(L,(len(L)-1)/2)
dspComplex(canvas, R, 1,1,"red" )
IR=adjustUnit(listIDTFT(R))
dspComplex(canvas, IR,50,1,"blue" )

⑤のこぎり波信号

　上記④の3行目を以下のように変更
L=[-1,0,1]# 中央を0番目とする。リストの外側は0とする。

■７．３　離散フーリエ変換（プログラム例なし）

■７．４　離散フーリエ変換のプログラム

表7-2  離散的フーリエ変換共通（DFTcom.py）

from graph import * #graph.pyについては表4-1参照
def initScaleDFT(title="Numerical Analysis",
                 XST=-1, YST=-3, DX=0.1,DY=0.1, BX=5,BY=10,
                 SX=10,SY=2,NX=20,NY=60,LX=5,LY=10,
                 Xfm="%.0f",Yfm="%.1f"):#初期スケールの表示
    root=initTk(title)
    canvas=initCanvas(root,500,300)
    drawScale(canvas, 40,40, 15, 10, SX, SY, XST, YST,
          DX, DY, BX,BY, LX, LY, NX, NY,Xfm,Yfm)
    canvas.create_text(60,30,text="実部",font=('Courier',10))
    drawScale(canvas, 280,40, 15, 10, SX, SY, XST, YST,
          DX, DY, BX,BY, LX, LY, NX, NY,Xfm,Yfm)
    canvas.create_text(300,30,text="虚部",font=('Courier',10))
    return canvas
def dspXY(canvas, R,XSC,YSC,CL,width=3,X0=40,Y0=100,Xscale=100,Yscale=20):
    N=len(R); X2=0
    for i in range(1, N):
        X1=X2;X2=i*XSC
        Y1=R[i-1].real*YSC;Y2=R[i].real*YSC
        if (Y1>=-3 and Y1<=3)and(Y2>=-3 and Y2<=3):
            drawV(canvas, X1,Y1,X2,Y2,CL,X0=X0,Y0=Y0,
                  width=width,Xscale=Xscale,Yscale=Yscale)
        Y1=R[i-1].imag*YSC; Y2=R[i].imag*YSC
        if (Y1>=-3 and Y1<=3)and(Y2>=-3 and Y2<=3):
            drawV(canvas, X1,Y1,X2,Y2,CL,X0=X0+240,Y0=Y0,
                  width=width,Xscale=Xscale,Yscale=Yscale)
def setRange(Y):
    if Y< -3: return -3
    if Y>  3: return  3
    return Y
def dspDTF(canvas, R,XSC,YSC,CL,width=3,X0=40, Y0=100,
           Xscale=100,Yscale=20,Nmax=100000,Nmin=0):#複素数のグラフ化
    N=Nmax; NN=len(R)
    if NN<Nmax: N=NN
    for i in range(Nmin, N):
        X=(i-Nmin)*XSC; Y=setRange(R[i].real*YSC)
        drawV(canvas, X,0,X,Y,CL,X0=X0,Y0=Y0,
              width=width,Xscale=Xscale,Yscale=Yscale)
        Y = setRange(R[i].imag*YSC)
        drawV(canvas, X,0,X,Y,CL,X0=X0+240,Y0=Y0,
              width=width,Xscale=Xscale,Yscale=Yscale)
def nomalScale(SX,DX,title="DFT"):#通常スケール表示
    return initScaleDFT(XST=0,SX=SX,DX=DX,Yfm="%.0f",title=title)
def centerScale(N,SX,DX,title="DFT"):#中央スケール表示
    return initScaleDFT(XST=-N/2,SX=SX,DX=DX,Yfm="%.0f",title=title)
def expandScale(title="DFT"):#先頭拡大スケール表示
    return initScaleDFT(XST= 0,SX=20,DX=1,NX=10,LX=1,
                        YST=-1.5,SY=5,DY=0.1,NY=30,
                        LY=5,BY=5,title=title)
def expandCenterScale(title="DFT"):#中央拡大スケール表示
    #NC=N/2
    return initScaleDFT(XST=-5,SX=20,DX=1,NX=10,LX=1,
                        YST=-1.5,SY=5,DY=0.1,NY=30,
                        LY=5,BY=5,title=title)
def dspDTFexCenter(canvas,S,N,color="blue"):#中央拡大表示
    dspDTF(canvas, S,2/N*10,2/N,color,
           X0=40,Y0=115,Yscale=50,Nmin=int(N/2-5))
def dspDTFexpand(canvas,S,N,color="blue"):#先頭拡大表示
    dspDTF(canvas, S,2/N*10,2/N,color,
           X0=20,Y0=115,Yscale=50)
def dspDTFnormal(canvas,S,N,color="blue"):#通常表示
    dspDTF(canvas, S,2/N,2/N,color)
def dspDTFcenter(canvas,S,N,color="blue"):#中央表示
    dspDTF(canvas, S,2/N,2/N,color,width=3,Nmin=int(N/2-10),X0=120)
    
def printDT(S, X):#データ表示
    print("■ ",S);num=len(X)
    for i in range(num): print(i+1, X[i].real,",", X[i].imag)
def expJ(TH):#オイラーの公式
    return complex(math.cos(TH),math.sin(TH))
def centerDT(X):# 中央を0位置とする並び替え
    N=len(X); NC=int((N-1)/2); L=[]
    for i in range(1,NC): L.append(X[NC+i])
    for i in range(NC+2):L.append(X[i])
    return L
def setDT(N, S):# データ設定
    # Sの形式  [[振幅, 周波数, 関数名], ・・・]
    X=[]; DTH=2.0*math.pi/N; ND=len(S)
    for i in range(N+1):
        DT=0
        for j in range(ND):
            DT+=S[j][0]*S[j][2](DTH*i*S[j][1])
        X.append(complex(DT,0))
    return X


表7-3  離散的フーリエ変換（DFT.py）
	
from DFTcom import *# 表7-2参照
def DFT(X):
    N=len(X);num=N-1; DN=-2.0*math.pi/num; Y=[]
    for k in range(num):
        YD=complex(0)
        for N in range(num):YD+=X[N]*expJ(DN*N*k)
        Y.append(YD)
    Y.append(Y[0])
    return Y
def IDFT(X):
    N=len(X);num=N-1; DN=2.0*math.pi/num; Y=[]
    for k in range(num):
        YD=complex(0)
        for N in range(num):YD+=X[N]*expJ(DN*N*k)
        Y.append(YD)
    for k in range(num): Y[k]/=num
    Y.append(Y[0])
    return Y

（DFTとIDFTの共通化)

def DFT_IDFT(X,iflag):#iflag=TrueのときIDFT
    N=len(X);num=N-1; DN=-2.0*math.pi/num; Y=[]
    if iflag: DN=-DN
    for k in range(num):
        YD=complex(0)
        for N in range(num):YD+=X[N]*expJ(DN*N*k)
        Y.append(YD)
    if iflag:
        for k in range(num): Y[k]/=num
    Y.append(Y[0])
    return Y

（４）合成波のDFTとIDFT

from DFT import * #DFT.pyについては表7-2参照
DT=[[1.0, 1, math.cos],[0.3, 3, math.cos],
    [0.6, 4, math.sin],[0.8, 2, math.sin]]
ND=20;N=100;SX=int(200/ND);DX=int(N/ND)
canvas1=nomalScale(SX,DX,title="データおよびDFT結果")
canvas2=nomalScale(SX,DX,title="IDFT結果")
X=setDT(N,DT); dspXY(canvas1, X,2/N,1,"red")      #データ表示
Y1=DFT(X); dspDTFnormal(canvas1,Y1,N,color="blue")#DTF結果表示
Y2=IDFT(Y1);dspXY(canvas2, Y2,2/N,1,"green")      #IDTF結果表示

【追加行】DFT結果の中央移動と中央近辺の拡大
canvas4=expandCenterScale()
S=centerDT(Y1);
dspDTFexCenter(canvas4,S,N,color="indigo")

（５）２次元フーリエ変換


表7-4  2次元離散的フーリエ変換
	
from tkinter import *
import math
def initDTFdraw(Title="Simple Image"):
    root=Tk();root.title(Title);return root
def initDTFcanvas(root,width=400, height=400, bg='white'):
    canvas=Canvas(root,width=width,height=height,bg=bg)
    canvas.pack();return canvas
def array(N1,N2=0):
    if N2==0: return [0 for i in range(N1)]
    return [[0 for i in range(N2)] for j in range(N1)]
#■以下，画像表示用
def Cdata(R):
    RR=int(R+0.5)
    if RR>255:return 255
    if RR<0  :return 0
    return RR
def RGB(R, G, B):return "#%02X%02X%02X" % (Cdata(R),Cdata(G),Cdata(B))
def setPM(V):#正負を設定
    if V>=0: return 1
    return -1
def dspSpe(canvas, DT, flag=0,X=0,Y=0):#スペクトラムの表示
    H=len(DT); W=len(DT[0]);spec=array(H,W)
    log10=math.log(10.0)#極端なスペクトル値を対数値で緩和する
    for i in range(H):
        for j in range(W):
            DR=DT[i][j].real; DI=DT[i][j].imag
            IR=setPM(DR); II=setPM(DI);DR=abs(DR); DI=abs(DI);           
            if   flag==0: spec[i][j] = IR*math.log(DR+1)/log10
            elif flag==1: spec[i][j] = II*math.log(DI+1)/log10
            else        : spec[i][j] = IR*II*math.log(DR+DI+1)/log10
    maxS=spec[0][0]; minS=spec[0][0] # 最大・最小を求める
    for i in range(H):
        for j in range(W):
            if spec[i][j] > maxS:maxS = spec[i][j];
            if spec[i][j] < minS:minS = spec[i][j];
    for i in range(H):#比例配分で表示
        for j in range(W):
            C =(spec[i][j]-minS)/(maxS-minS)*255
            CL=RGB(C,C,C); Xj=X+j; Yi=Y+i
            canvas.create_rectangle(Xj,Yi,Xj,Yi,fill=CL,outline='')
def drawDT(cv,R,flag, X=0,Y=0,scale=50, base=128):
    height=len(R);width=len(R[0])
    for j in range(height):
        for i in range(width):
            if flag: V=R[j][i].imag*scale+base
            else:    V=R[j][i].real*scale+base
            CL=RGB(V,V,V)
            cv.create_rectangle(X+i,Y+j,X+i,Y+j,fill=CL,outline='')
#■以下，DFT2D
def expJ(TH):return complex(math.cos(TH), math.sin(TH))
def modS(K):#奇数のとき-1. 偶数のとき1
    if K & 1: return -1
    else: return 1
def moveCenter(X):#中央へ移動(高速化目的なら不要。表示のとき行う)
    H=len(X);W=len(X[0]); Y=array(H,W)
    for k in range(H):#中央への移動
        for N in range(W):
            Y[k][N] = X[k][N]*modS(k+N)
    return Y
def Xtran(Y,ARH):#X方向変換
    H=len(Y);W=len(Y[0]);tm=array(H,W)
    for i in range(H):
        for k in range(W):
            tm[i][k]=complex(0,0)
            for N in range(H):
                tm[i][k] += Y[i][N]*expJ(ARH*N*k);
    return tm
def Ytran(tm,ARW):#Y方向変換
    H=len(tm);W=len(tm[0]);Y=array(H,W)
    for i in range(W):#
        for k in range(H):
            Y[k][i]=complex(0,0)
            for N in range(H):
                Y[k][i] += tm[N][i]*expJ(ARW*N*k)
    return Y
def divArray(Z,H,W,V):#j除算
    for i in range(H):
        for k in range(H): Z[i][k]/= V
    return Z
def moveAround(Z,tm,W,H):
    for k in range(H): #周辺への移動
        for N in range(W):Z[k][N]=tm[k][N]*modS(k+N)
    return Z
def DFT2D(X): #■フーリエ変換
    H=len(X);W=len(X[0]);P2=-2*math.pi;ARH=P2/H;ARW=P2/W
    Y=moveCenter(X)       #中央への移動
    tm=Xtran(Y,ARH)       #X方向変換
    return Ytran(tm,ARW)  #Y方向変換
def IDFT2D(Y):#■逆フーリエ変換
    H=len(Y);W=len(Y[0]);P2=2*math.pi
    tm =array(H,W);  Z=array(H,W); ARH=P2/H;ARW=P2/W
    Z=Xtran(Y,ARH)        #X方向変換
    Z=divArray(Z,H,W,H)   #Hで除算
    tm=Ytran(Z,ARW)       #Y方向変換
    tm=divArray(tm,H,W,W) #Wで除算
    return moveAround(Z,tm,W,H)#周辺に移動
#■以下データ設定用
def waveSynthe(L,i,N):#波の合成
    V=0
    for k in range(len(L)):
        V += L[k][0]*L[k][2](L[k][1]*i*math.pi/N)
    return V
def set2D(L,NX,NY):
    NN=NX
    if NY>NN: NN=NY #現時点では正方行列とする
    R=array(NN,NN)
    for i in range(NX):
        VX =waveSynthe(L[0],i,NX)
        for j in range(NY):
            VY =waveSynthe(L[1],j,NY)
            R[j][i]=complex(VX+VY)
    return R
root=initDTFdraw()#メインプログラム
canvas=initDTFcanvas(root,width=400,height=200)
L=[[[1,4,math.cos],[1, 8,math.cos],[1, 16,math.cos]],
   [[1,6,math.sin],[1, 12,math.sin],[1, 24,math.sin]]]
B=set2D(L,100,100); drawDT(canvas,B,False,20,20,scale=20)
A=DFT2D(B)        ; dspSpe(canvas,A,flag=2,X=130,Y=20)
Z=IDFT2D(A)       ;drawDT(canvas,Z,False,240,20,scale=20)
	
■７．５　高速フーリエ変換

表7-3 時間間引きアルゴリズム（FFT.pyとして保存する）

from DFT import * #DFTのグラフ関数を借用する
import copy
#■以下グラフ表示用
def Scale1(title="FFT"):
    return initScaleDFT(XST= 0,SX=25,DX=1,NX=8,BX=4,LX=1,
                        YST=-4,SY=4,DY=0.2,NY=40,
                        LY=5,BY=5,title=title)
def Scale2(title="データ"):
    return initScaleDFT(XST= 0,SX=25,DX=1,NX=8,BX=4,LX=1,
                        YST= 0,SY=8,DY=0.1,NY=20,
                        LY=10,BY=5,title=title)
def Scale3(SX,DX,title="データ"):
    return initScaleDFT(XST=0,SX=SX,DX=DX,Yfm="%.0f",title=title)
def dspFFT1(canvas,S, color="blue"):
    dspDTF(canvas, S,1,1,color,Xscale=25,
           X0=40,Y0=120,Yscale=20,Nmin=0)
def dspXYFFT(canvas, R,XSC,YSC,CL,width=3):
    dspXY(canvas,R,XSC,YSC,CL,Y0=200,Yscale=80)
#■以下，高速フーリエ変換
def FFT(X):#FFT（時間間引きアルゴリズム）
    m=len(X)-1; Y=[]; DH=-2*math.pi/m
    for k in range(m+1):Y.append(X[k])
    for j in range(m):#ビット逆順の並び替え
        jp = j; jx = 0; k = 1
        while k < m:
            BT = jp & 1 ;jx = (jx <<1) | BT
            jp >>= 1 ; k <<= 1
        if j < jx :YY = Y[j]; Y[j] = Y[jx]; Y[jx] = YY
    n_half = 2;  ne = m >>1	# FFTの計算
    while ne >= 1:
        n_half2 = n_half >> 1
        for jp in range(0, m, n_half):
            jx = 0
            for j  in range(jp, jp+n_half2):
                jnh = j + n_half2
                yTemp = Y[jnh]*expJ(DH*jx)#バタフライ計算
                Y[jnh] = Y[j]-yTemp;Y[j] = Y[j]+yTemp
                jx += ne
        n_half <<= 1; ne >>= 1
    Y[m]=Y[0]
    return Y
def IFFT(X):#FFT（時間間引きアルゴリズム）
    m=len(X)-1; Y=[]; DH=2*math.pi/m
    for k in range(m+1):Y.append(X[k])
    for j in range(m):#ビット逆順の並び替え
        jp = j; jx = 0; k = 1
        while k < m:
            BT = jp & 1 ;jx = (jx <<1) | BT
            jp >>= 1 ; k <<= 1
        if j < jx :
            print(j,jx)
            YY = Y[j]; Y[j] = Y[jx]; Y[jx] = YY
    n_half = 2;  ne = m >>1	# FFTの計算
    while ne >= 1:
        n_half2 = n_half >> 1
        for jp in range(0, m, n_half):
            jx = 0
            for j  in range(jp, jp+n_half2):
                jnh = j + n_half2
                yTemp = Y[jnh]*expJ(DH*jx)#バタフライ計算
                Y[jnh] = Y[j]-yTemp;Y[j] = Y[j]+yTemp
                jx += ne
        n_half <<= 1; ne >>= 1
    for i in range(m):Y[i]/=m
    Y[m]=Y[0]
    return Y

■結果の確認

X=[0,0,0,1,1,1,0,0,0] # 配列の大きさはデータ数+1とする。
canvas1=Scale2(title="FFT")
canvas2=Scale1(title="データ")
canvas3=Scale2(title="IFFT")
dspXYFFT(canvas1, X,2/8,1,"red")
Y=FFT(X) ; dspFFT1(canvas2, Y)
Z=IFFT(Y); dspXYFFT(canvas3, Z,2/8,1,"green")

表7-4 時間間引きアルゴリズムの改良

def FFT_IFFT(X, IFFTflag):#FFT（時間間引きアルゴリズム）
    m=len(X)-1; Y=[]; DH=-2*math.pi/m
    if IFFTflag:DH=-DH
    for k in range(m+1):Y.append(X[k])
    for j in range(m):#ビット逆順の並び替え
        jp = j; jx = 0; k = 1
        while k < m:
            BT = jp & 1 ;jx = (jx <<1) | BT
            
            jp >>= 1 ; k <<= 1
        if j < jx :
            print(j,jx)
            YY = Y[j]; Y[j] = Y[jx]; Y[jx] = YY
    n_half = 2;  ne = m >>1	# FFTの計算
    while ne >= 1:
        n_half2 = n_half >> 1
        for jp in range(0, m, n_half):
            jx = 0
            for j  in range(jp, jp+n_half2):
                jnh = j + n_half2
                yTemp = Y[jnh]*expJ(DH*jx)#バタフライ計算
                Y[jnh] = Y[j]-yTemp;Y[j] = Y[j]+yTemp
                jx += ne
        n_half <<= 1; ne >>= 1
    if IFFTflag:
        for i in range(m):Y[i]/=m
    Y[m]=Y[0]
    return Y

表7-5 周波数間引きアルゴリズム

def FFT_IFFTfreq(X, IFFTflag):#FFT（周波数間引きアルゴリズム）
    n=len(X)-1; Y=[]; DH=-2*math.pi/n; n_half= n>>1
    if IFFTflag:DH=-DH
    ne = 1	# FFTの計算
    while ne<n:
        n_half2 = n_half <<1
        for jp in range(0,n,n_half2):
            jx = 0
            for j in range(jp,jp + n_half):
                jnh = j+n_half; yTemp = X[j]; X[j] = yTemp+X[jnh]
                X[jnh]=(yTemp-X[jnh])*expJ(DH * jx)
                jx += ne
        n_half >>= 1;  ne <<= 1
    for j in range(n):	# ビット逆順の並び替え
        jp = j; jx = 0; k = 1
        while k < n:
            BT = jp & 1; jx = (jx<<1)|BT; jp>>=1; k<<=1
        if j < jx : yTemp=X[j]; X[j]=X[jx]; X[jx] = yTemp
    X[n] = X[0]
    if IFFTflag:
        for i in range(n):X[i]/=n
    X[n] = X[0]
    return X

表7-6 回転因子とビット逆順の値をあらかじめ計算しておく方法

W=128
WH=int(W>>1)
w_tbl=[0 for i in range(WH)] # 回転因子計算用テーブル
b_tbl=[0 for i in range(W)] # ビット逆順交換用テーブル
n_old=0  # テーブルが作成されているかどうかの確認用変数
def FFT_tbl(xx, n_FFT):
    global n_old, w_tbl, b_tbl
    num=abs(n_FFT);
    x=copy.copy(xx) 
    if(num !=n_old): #1回目の呼び出しのときのみテーブル作成
        pow2=1       #要素数が2のべき乗かどうかを判定
        for j in range(32):
            if num==pow2 : break
            else :pow2<<=1
        if j>=32:
            print("\n データ数が2のべき乗になっていません")
            return []
        argV = -2.0 * math.pi /num #回転因子計算用テーブル作成
        for j in range(int(num/2)): w_tbl[j]=expJ(argV*j)
        n_half=num>>1; b_tbl[0]=0 #ビット逆順交換用テーブル作成
        ne=1
        while ne<num:
            for jp in range(ne):
                b_tbl[jp+ne]=b_tbl[jp]+n_half
            n_half>>=1
            ne<<=1
        n_old=num;
    n_half=num>>1 #FFT計算開始
    ne=1
    while ne<num:
        n_half2=n_half<<1
        jp=0
        while jp<num:
            jx=0
            for j in range(jp, jp+n_half):
                jnh=j + n_half; xtmp = x[j]; x[j] = xtmp+x[jnh]
                x[jnh]= (xtmp-x[jnh])* w_tbl[jx]
                jx += ne
            jp += n_half2
        n_half>>=1
        ne<<=1
    for j in range(num):#ビット逆順入れ換え
        if j<b_tbl[j]:
            xx=x[j]; x[j]=x[b_tbl[j]]; x[b_tbl[j]]=xx
    if n_FFT<0:
        for j in range(num): x[j]/= num #逆FFT
    return x

表7-7  2次元高速フーリエ変換

def fft_2d(beforDFT, n_fft, W, H):
    afterDFT=[[0 for j in range(H+1)] for i in range(W+1)]
    num=abs(n_fft)
    u=[0 for i in range(W+1)]
    for n1 in range(num):
        for n2 in range(num):
            if (n1+n2) & 1 : DD = -1
            else           : DD =  1
            afterDFT[n1][n2]= beforDFT[n1][n2]*DD
        afterDFT[n1]=FFT_tbl(afterDFT[n1], n_fft)
        if afterDFT[n1]==[]: return
    for n2 in range(num):
        for n1 in range(num): u[n1]=afterDFT[n1][n2]
        u=FFT_tbl(u, n_fft)
        for n1 in range(num): afterDFT[n1][n2]=u[n1];
    return afterDFT

表7-8  高速化2次元FFT

def fft_2dX(beforDFT, n_fft, W, H):
    num=abs(n_fft)
    u=[0 for i in range(W+1)]
    afterDFT=copy.copy(beforDFT)
    for n1 in range(num):
        afterDFT[n1]=FFT_tbl(afterDFT[n1], n_fft)
        if(afterDFT[n1]==[]) :return
    for n2 in range(num):
        for n1 in range(num): u[n1]=afterDFT[n1][n2]
        u=FFT_tbl(u,n_fft)
        for n1 in range(num): afterDFT[n1][n2]=u[n1]
    return afterDFT

表7-9  スペクトル位置の移動

def setDspSpec(DT, W, H):
    TM=[[0 for j in range(H+1)] for i in range(W+1)]
    w2=int(W/2); h2=int(H/2)
    for x in range(W):
        for y in range(H):
            x2=x+w2; y2=y+h2
            if x2>W: x2-= W
            if y2>H: y2-= H
            TM[x2][y2]=DT[x][y]
    return TM

表7-10  高速化2次元FFT

def rev_fft_2d(afterDFT, n_fft, W, H):
    num=abs(n_fft);num2=int(num/2)
    rDT=fft_2d(afterDFT,-num, W, H)
    for k in range(num+1):
        N=0; N2=num-1
        while N<num2:
            rr = rDT[k][N]; rDT[k][N]=rDT[k][N2];rDT[k][N2]=rr
            N+=1; N2-=1
    for k in range(num+1):
        N=0; N2=num-1
        while N<num2:
            rr=rDT[N][k];rDT[N][k]=rDT[N2][k]; rDT[N2][k]=rr
            N+=1; N2-=1
    for i in range(int(num2)):
        for j in range(int(num2)):
            ii=i+num2; jj=j+num2
            rr=rDT[i][j];rDT[i][j]=rDT[ii][jj];rDT[ii][jj]=rr
            rr=rDT[ii][j];rDT[ii][j]=rDT[i][jj]; rDT[i][jj]=rr
    return rDT
 
■７．６　離散的コサイン変換

表7-11  DCTプログラム

import math
def array(N1, N2=0):
    if N2!=0: return [[0 for i in range(N2)] for i in range(N1)]
    return [0 for i in range(N1)]
def DCT(X):  #DCT(離散的コサイン変換)
    Num=len(X);Y=array(Num);Pn=math.pi/(2*Num)
    C0 = 1 / math.sqrt(Num); Cn0 = math.sqrt(2)*C0
    for k in range(Num):
        for j in range(Num): Y[k]+=X[j]*math.cos((2*j+1)*k*Pn)
        if k == 0 : Y[k] = C0 * Y[k]
        else      : Y[k] = Cn0 * Y[k]
    return Y
def IDCT(Y): #IDCT(離散的逆変換)
    Num=len(X);Z=array(Num);Pn=math.pi/(2*Num)
    C0 = 1 / math.sqrt(2); C = math.sqrt(2/Num)
    for j in range(Num):
        Z[j] = C0 * Y[0]
        for k in range(1,Num): Z[j]+=Y[k]*math.cos((2*j+1)*k*Pn)
        Z[j] = C * Z[j]
    return Z
def printXYZ(X,Y,Z):
    print("\n  No   Data     DCT      IDCT")
    print("-"*30); N=len(X)
    for i in range(N):
        print("%3d %8.5f %8.5f %8.5f"%(i,X[i],Y[i],Z[i]))
X=[0,0,0,1,1,1,0,0,0]
Y=DCT(X);Z=IDCT(Y)
printXYZ(X,Y,Z)


表7-12  Chenのアルゴリズムによるプログラム

def FDCT(X): #FDCT(Chenのアルゴリズム)
    PI=math.pi; Num=len(X)
    Cos4  =math.cos(PI/ 4); Cos8  =math.cos(PI/8)
    Cos16 =math.cos(PI/16); Cos163=math.cos(3*PI/16)
    Sin8  =math.sin(PI/ 8); Sin16 =math.sin(PI/16)
    Sin163=math.sin(3*PI/16)
    Y=[]
    for j in range(Num): Y.append(X[j])
    for J in range(4):
        tmp = Y[J]+Y[7-J];Y[7-J]=Y[J]-Y[7-J];Y[J]=tmp
    tmp = Y[0]+Y[3]; Y[3]=Y[0]-Y[3]; Y[0]=tmp
    tmp = Y[1]+Y[2]; Y[2]=Y[1]-Y[2]; Y[1]=tmp
    tmp = Cos4*(-Y[5]+Y[6])
    Y[6]= Cos4*( Y[5]+Y[6]);Y[5]= tmp
    tmp = Cos4*(Y[0]+Y[1])
    Y[1]= Cos4*(Y[0]-Y[1]);Y[0]= tmp
    tmp =  Sin8*Y[2]+Cos8*Y[3]
    Y[3]= -Cos8*Y[2]+Sin8*Y[3];Y[2]= tmp
    tmp = Y[4]+Y[5];Y[5]=Y[4]-Y[5];Y[4]=tmp
    tmp =-Y[6]+Y[7];Y[7]=Y[6]+Y[7];Y[6]=tmp
    tmp = Sin16*Y[4]+Cos16*Y[7]
    Y[7]=-Cos16*Y[4]+Sin16*Y[7];Y[4]=tmp
    tmp = Cos163*Y[5]+Sin163*Y[6]
    Y[6]=-Sin163*Y[5]+Cos163*Y[6];Y[5]=tmp

    tmp =Y[1]; Y[1]=Y[4]; Y[4]=tmp
    tmp =Y[3]; Y[3]=Y[6]; Y[6]=tmp
    for j in range(Num): Y[j] *= 0.5
    return Y

表7-13  2次元DCT

from tkinter import *
import math
#■画像表示用
def Cdata(R):
    RR=int(R+0.5)
    if RR>255:return 255
    if RR<0  :return 0
    return RR
def RGB(R, G, B):
return "#%02X%02X%02X" % (Cdata(R),Cdata(G),Cdata(B))
def drawDTreal(cv,R,X=0,Y=0,scale=50, base=128):
    height=len(R);width=len(R[0])
    for j in range(height):
        for i in range(width):
            V=R[j][i]*scale+base;CL=RGB(V,V,V)
            cv.create_rectangle(X+i,Y+j,X+i,Y+j,fill=CL,outline='')
#■2次元DCT
def array(N1,N2):
    return [[0 for i in range(N2)] for j in range(N1)]
def DCT2D(dt, inv ):
    W=len(dt);H=len(dt[0]); out=array(W,H)
    trt=array(W,H); wrt=array(W,W); hrt=array(H,H) #作業用配列

    if inv : nm1 =  0.5; nm2 = 2.0/W ; stID   = 1
    else   : nm1 =  0.0; nm2 = 1.0   ; stID   = 0
    M2    = math.pi/W
    for i in range(W):
        for j in range(stID,W):
            if inv: wrt[i][j]=math.cos(M2*(i+0.5)*j)
            else  : wrt[i][j]=math.cos(M2*(j+0.5)*i)
    for y in range(H):
        for x in range(W):
            trt[x][y]=0.0
            for xx in range(stID,W):trt[x][y]+=dt[xx][y]*wrt[x][xx]
            trt[x][y] += dt[0][y]*nm1; trt[x][y] *= nm2
    M2 = math.pi/H
    if inv : nm2=2.0/H
    else   : nm=1
    for i in range(W):
        for j in range(stID,H):
            if inv: hrt[i][j]=math.cos(M2*(i+0.5)*j)
            else  : hrt[i][j]=math.cos(M2*(j+0.5)*i)
    for x in range(W):
        for y in range(H):
            out[x][y]=0.0
            for yy in range(stID,H):
                out[x][y]+=trt[x][yy]*hrt[y][yy]
            out[x][y] += trt[x][0]*nm1;	out[x][y] *= nm2;
    return out
#■データ設定用
def waveSynthe(L,i,N):#波の合成
    V=0
    for k in range(len(L)):
        V += L[k][0]*L[k][2](L[k][1]*i*math.pi/N)
    return V
def set2Dreal(L,NX,NY):
    NN=NX
    if NY>NN: NN=NY #現時点では正方行列とする
    R=array(NN,NN)
    for i in range(NX):
        VX =waveSynthe(L[0],i,NX)
        for j in range(NY):
            VY =waveSynthe(L[1],j,NY)
            R[j][i]=VX+VY
    return R

#■メインプログラム
root=Tk();root.title("2D DCT")
canvas=Canvas(root,width=400,height=200,bg="white")
canvas.pack()
L=[[[1,4,math.cos],[1, 8,math.cos],[1, 16,math.cos]],
   [[1,6,math.sin],[1, 12,math.sin],[1, 24,math.sin]]]
B=set2Dreal(L,100,100); drawDTreal(canvas,B,20,20,scale=20)
A=DCT2D(B,False);drawDTreal(canvas,A,140,20,scale=20)
Z=DCT2D(A,True );drawDTreal(canvas,Z,260,20,scale=20)

■７．７　アダマール変換

表7-14  アダマール行列の生成

def array(N1, N2=0):
    if N2!=0: return [[0 for i in range(N2)] for i in range(N1)]
    return [0 for i in range(N1)]
def HT_Mat(N):
    MX = (1<<N) # 2のN乗 
    HMat=array(MX, MX); HMat[0][0] = 1
    K = 1
    for L in range(1,N+1):
        for i in range(K):
            for j in range(K):
                HMat[K + i][j] = HMat[i][j]
                HMat[i][K + j] = HMat[i][j]
                HMat[K + i][K + j] = -HMat[i][j]
        K <<=1 # K *= 2の替わり
    return HMat
 
表7-15  符号変化回数によるアダマール行列の行入替え

def CountChange(HMAT,K):
    MX=len(HMAT);CC = 0
    ID = HMAT[K][0]
    for j in range(1,MX):
        if ID != HMAT[K][j]: CC+=1;ID=HMAT[K][j]
    return CC
def HT_Sort(HMAT):
    MX = len(HMAT)
    for i in range(MX-1):
        MXC = CountChange(HMAT, i);ID = i
        for j in range(i+1,MX):
            MXP = CountChange(HMAT,j)
            if MXP < MXC : ID=j; MXC=MXP
        if ID!=i : T=HMAT[ID]; HMAT[ID]=HMAT[i]; HMAT[i]=T

表7-16  アダマール変換プログラム

---（前略：先頭に表7-14, 表7-15のプログラムを付加する）---
def HT(HMAT,X): #Hadamard変換
    Num=len(X); Y=array(Num)
    for i in range(Num):
        for j in range(Num):
            if HMAT[i][j]>0: Y[i]+=X[j] #行列値は1と-1のみ
            else           : Y[i]-=X[j]
    return Y
def IHT(HMAT, X):#逆Hadamard変換
    Num=len(X); Y=HT(HMAT,X)
    for i in range(Num):Y[i]/=Num
    return Y
def dspHmat(Mat):
    N=len(Mat)
    for i in range(N):
        S=""
        for j in range(N): S+= "%3d" % Mat[i][j]
        print(S)
def printXYZ(X,Y,Z):
    print("\n  No    X         Y         Z")
    print("-"*30); N=len(X)
    for i in range(N):
        print("%3d %8.5f %8.5f %8.5f"%(i,X[i],Y[i],Z[i]))
HMT=HT_Mat(3); HT_Sort(HMT); dspHmat(HMT)
X=[0,2,0,1,-1,0,1,0]
Y=HT(HMT,X); 
Z=IHT(HMT,Y);
printXYZ(X,Y,Z)

表7-17  高速アダマール変換プログラム

def FHT(X):  #高速Hadamard変換
    Num=len(X); Y=[]
    for i in range(Num):Y.append(X[i])
    #ビット逆順の並び替え
    for j in range(Num):
        JP = j; JX = 0; k = 1
        while k<Num:
            MD = JP & 1 # Mod 2の替わり
            JX = (JX <<1) + MD
            JP >>= 1; k<<= 1  # >>1は2で除算,<<1は2を乗算
        if j<JX: T= Y[j];Y[j]= Y[JX];Y[JX]=T
 
    #FHTの計算（乗算がないことに注意）
    n_half = Num >>1; Ne = 1
    while Ne < Num:
        SG=True
        n_half2 = n_half<<1
        for JP in range(0,Num,n_half2):
            JX = 0
            for J in range(JP, JP + n_half):
                jnh = J + n_half
                if SG: xtmp= Y[jnh]
                else : xtmp=-Y[jnh]
                Y[jnh] = Y[J] - xtmp
                Y[J  ] = Y[J] + xtmp
                JX += Ne
            SG = not SG
        n_half >>=1;  Ne <<= 1
    return Y
def IFHT(X):#■逆高速Hadamard変換
#画像データのようにすぺてのデータが整数値のとき
#Numの値が2のN乗となるので，除算もシフト演算で
#置き換えることができる。
    Num=len(X); Y=FHT(X)
    for i in range(Num):Y[i]/=Num
    return Y
def printXYZ(X,Y,Z):
    print("\n  No    X         Y         Z")
    print("-"*30); N=len(X)
    for i in range(N):
        print("%3d %8.5f %8.5f %8.5f"%(i,X[i],Y[i],Z[i]))
X=[0,2,0,1,-1,0,1,0]
Y=FHT(X)
Z=IFHT(Y)
printXYZ(X,Y,Z)