1: //========================================================================================
     2: //  FFT.java
     3: //    en:Fast Fourier Transform
     4: //    ja:高速フーリエ変換
     5: //  Copyright (C) 2003-2019 Makoto Kamada
     6: //
     7: //  This file is part of the XEiJ (X68000 Emulator in Java).
     8: //  You can use, modify and redistribute the XEiJ if the conditions are met.
     9: //  Read the XEiJ License for more details.
    10: //  https://stdkmd.net/xeij/
    11: //========================================================================================
    12: 
    13: //----------------------------------------------------------------------------------------
    14: //  フーリエ変換
    15: //    X[m];0<=m<Nがx[n];0<=n<Nのフーリエ変換のとき
    16: //    X[m]=sum[0<=n<N]{ω[N]^(-m*n)*x[n]}
    17: //
    18: //  回転子
    19: //    ω[N]=e^(2*pi*i/N)=cos(2*pi/N)+i*sin(2*pi/N)
    20: //    ω[N]^n=ω[N/n]
    21: //    ω[N]^(-n)=conj(ω[N]^n)=conj(ω[N/n])
    22: //    ω[N]^N=ω[1]=1
    23: //    ω[N]^(N/2)=ω[2]=-1
    24: //    ω[N]^(N/4)=ω[4]=i
    25: //
    26: //  Sande-Tukeyの周波数間引き型FFT(基数K)
    27: //    NはKの倍数とする
    28: //    X[m];0<=m<Nをx[n];0<=n<Nのフーリエ変換とする
    29: //    X[m]=sum[0<=n<N]{ω[N]^(-m*n)*x[n]}
    30: //    周波数のインデックスmをKで割った商をj、余りをuとする
    31: //    m=K*j+u    0<=j<N/K, 0<=u<K
    32: //    時間のインデックスnをN/Kで割った商をv、余りをkとする
    33: //    n=k+N/K*v  0<=k<N/K, 0<=v<K
    34: //    X[m]=sum[0<=n<N]{ω[N]^(-m*n)*x[n]}
    35: //    X[K*j+u]=sum[0<=k<N/K]{sum[0<=v<K]{ω[N]^(-(K*j+u)*(k+N/K*v))*x[k+N/K*v]}}
    36: //            =sum[0<=k<N/K]{sum[0<=v<K]{ω[N]^(-K*j*k-N*j*v-u*k-N/K*u*v)*x[k+N/K*v]}}
    37: //            =sum[0<=k<N/K]{ω[N]^(-K*j*k)*ω[N]^(-u*k)*sum[0<=v<K]{ω[N]^(-N/K*u*v)*x[k+N/K*v]}}
    38: //            =sum[0<=k<N/K]{ω[N/K]^(-j*k)*ω[N]^(-u*k)*sum[0<=v<K]{ω[K]^(-u*v)*x[k+N/K*v]}}
    39: //    ここで
    40: //    X'[j,u]=X[K*j+u]                                         0<=j<N/K, 0<=u<K
    41: //    x'[k,u]=ω[N]^(-u*k)*sum[0<=v<K]{ω[K]^(-u*v)*x[k+N/K*v]}  0<=k<N/K, 0<=u<K
    42: //    とおくと、X'[j,u]はx'[k,u]のフーリエ変換に等しい
    43: //----------------------------------------------------------------------------------------
    44: 
    45: package xeij;
    46: 
    47: import java.lang.*;  //Boolean,Character,Class,Comparable,Double,Exception,Float,IllegalArgumentException,Integer,Long,Math,Number,Object,Runnable,SecurityException,String,StringBuilder,System
    48: 
    49: public class FFT {
    50: 
    51:   private int fftN;
    52:   private double[] fftSinCos;  //サインコサインテーブル。0..5*N/4
    53:   private int[] fftSwap;  //スワップリスト。0..N-1の範囲の整数でできるビット反転のペアを小さい順に列挙した配列。ビット反転しても値が変わらない数を含まない
    54: 
    55:   //コンストラクタ
    56:   //  Nは4以上の2の累乗に限る
    57:   public FFT (int n) {
    58:     fftN = n;
    59:     int n2 = n >> 1;  //N/2
    60:     int n4 = n >> 2;  //N/4
    61:     //サインコサインテーブルを作る
    62:     fftSinCos = new double[5 * n4 + 1];
    63:     double d = 2.0 * Math.PI / (double) n;  //Nは2の累乗なので先に割っても誤差が増えることはない
    64:     for (int i = 1; i < n4; i++) {
    65:       //          |<-------------cos------------->|
    66:       //  |<-------------sin------------->|
    67:       //  +-----*****-----+-------+-------+-----***
    68:       //  |  ***  |  ***  |       |       |  ***  |
    69:       //  |**     |     **|       |       |**     |
    70:       //  *-------+-------*-------+-------*-------+
    71:       //  0      N/4     N/2*   3*N/4   **N     5*N/4
    72:       //  |       |       |  ***  |  ***  |       |
    73:       //  +-------+-------+-----*****-----+-------+
    74:       //   ------> <------ ------> <------ ------>
    75:       //     (1)     (2)     (3)     (4)     (5)
    76:       //         (4)                 (3)                   (2)                (5)              (1)
    77:       fftSinCos[n - i] = fftSinCos[n2 + i] = -(fftSinCos[n2 - i] = fftSinCos[n + i] = fftSinCos[i] = Math.sin (d * (double) i));
    78:     }
    79:     fftSinCos[n] = fftSinCos[n2] = fftSinCos[0] = 0.0;  //代入が2方向から来て重複するのは構わないが0と±1には誤差が入らないようにしたい
    80:     fftSinCos[5 * n4] = fftSinCos[n4] = 1.0;  //[5*N/4]は現在使用していないが入れておく
    81:     fftSinCos[3 * n4] = -1.0;
    82:     //スワップリストを作る
    83:     int b = Integer.numberOfTrailingZeros (n);
    84:     fftSwap = new int[n - (1 << (b + 1 >> 1))];
    85:     for (int i = 0, k = 0; i < n; i++) {
    86:       int j = Integer.reverse (i << -b);
    87:       if (i < j) {
    88:         fftSwap[k] = i;
    89:         fftSwap[k + 1] = j;
    90:         k += 2;
    91:       }
    92:     }
    93:   }  //new FFT(int)
    94: 
    95:   //fft.fftSandeTukey2 (x, y)
    96:   //  Sande-Tukeyの周波数間引き型FFT(基数2)
    97:   //    X[K*j+u]=sum[0<=k<N/K]{ω[N/K]^(-j*k)*ω[N]^(-u*k)*sum[0<=v<K]{ω[K]^(-u*v)*x[k+N/K*v]}}
    98:   //    X[2*j]=sum[0<=k<N/2]{ω[N/2]^(-j*k)*{x[k]+x[k+N/2]}}
    99:   //    X[2*j+1]=sum[0<=k<N/2]{ω[N/2]^(-j*k)*ω[N]^(-k)*(x[k]-x[k+N/2])}
   100:   public void fftSandeTukey2 (double[] xr, double[] xi) {
   101:     int n4 = fftN >> 2;  //N/4
   102:     int half = fftN >> 1;  //N/2
   103:     for (int step = 1; half >= 1; half >>= 1, step <<= 1) {
   104:       for (int k = 0, theta = 0; k < half; k++, theta += step) {
   105:         double or = fftSinCos[theta + n4];  //ω[N]^(-k)
   106:         double oi = -fftSinCos[theta];
   107:         for (int start = 0; start < fftN; start += half << 1) {
   108:           int k0 = start + k;  //k
   109:           int k1 = k0 + half;  //k+N/2
   110:           double pr = xr[k0];  //x[k]
   111:           double pi = xi[k0];
   112:           double qr = xr[k1];  //x[k+N/2]
   113:           double qi = xi[k1];
   114:           xr[k0] = pr + qr;  //X[2*j]=x[k]+x[k+N/2]
   115:           xi[k0] = pi + qi;
   116:           pr -= qr;  //x[k]-x[k+N/2]
   117:           pi -= qi;
   118:           xr[k1] = or * pr - oi * pi;  //X[2*j+1]=ω[N]^(-k)*(x[k]-x[k+N/2])
   119:           xi[k1] = or * pi + oi * pr;
   120:         }
   121:       }
   122:     }
   123:     //インデックスのビット反転
   124:     for (int k = 0, l = fftSwap.length; k < l; k += 2) {
   125:       int i = fftSwap[k];
   126:       int j = fftSwap[k + 1];
   127:       double t = xr[i];
   128:       xr[i] = xr[j];
   129:       xr[j] = t;
   130:       t = xi[i];
   131:       xi[i] = xi[j];
   132:       xi[j] = t;
   133:     }
   134:   }  //fft.fftSandeTukey2(double[],double[])
   135: 
   136:   //fft.fftSandeTukey4 (x, y)
   137:   //  Sande-Tukeyの周波数間引き型FFT(基数4)
   138:   //    X[4*j+u]=sum[0<=k<N/4]{ω[N/4]^(-j*k)*ω[N]^(-u*k)*sum[0<=v<4]{i^(-u*v)*x[k+N/4*v]}}
   139:   //    X[4*j]=sum[0<=k<N/4]{ω[N/4]^(-j*k)*((x[k]+x[k+N/2])+(x[k+N/4]+x[k+3*N/4]))}
   140:   //    X[4*j+1]=sum[0<=k<N/4]{ω[N/4]^(-j*k)*ω[N]^(-k)*((x[k]-x[k+N/2])-i*(x[k+N/4]-x[k+3*N/4]))}
   141:   //    X[4*j+2]=sum[0<=k<N/4]{ω[N/4]^(-j*k)*ω[N]^(-2*k)*((x[k]+x[k+N/2])-(x[k+N/4]+x[k+3*N/4]))}
   142:   //    X[4*j+3]=sum[0<=k<N/4]{ω[N/4]^(-j*k)*ω[N]^(-3*k)*((x[k]-x[k+N/2])+i*(x[k+N/4]-x[k+3*N/4]))}
   143:   public void fftSandeTukey4 (double[] xr, double[] xi) {
   144:     int mask = fftN - 1;
   145:     int n4 = fftN >> 2;  //N/4
   146:     int quarter = n4;  //N/4
   147:     for (int step = 1; quarter >= 2; quarter >>= 2, step <<= 2) {
   148:       for (int k = 0, theta1 = 0; k < quarter; k++, theta1 += step) {
   149:         int theta2 = theta1 + theta1 & mask;
   150:         int theta3 = theta2 + theta1 & mask;
   151:         double or1 = fftSinCos[theta1 + n4];  //ω[N]^(-k)
   152:         double oi1 = -fftSinCos[theta1];
   153:         double or2 = fftSinCos[theta2 + n4];  //ω[N]^(-2*k)
   154:         double oi2 = -fftSinCos[theta2];
   155:         double or3 = fftSinCos[theta3 + n4];  //ω[N]^(-3*k)
   156:         double oi3 = -fftSinCos[theta3];
   157:         for (int start = 0; start < fftN; start += quarter << 2) {
   158:           int k0 = start + k;  //k
   159:           int k1 = k0 + quarter;  //k+N/4
   160:           int k2 = k1 + quarter;  //k+N/2
   161:           int k3 = k2 + quarter;  //k+3*N/4
   162:           double xr0 = xr[k0];  //x[k]
   163:           double xi0 = xi[k0];
   164:           double xr1 = xr[k1];  //x[k+N/4]
   165:           double xi1 = xi[k1];
   166:           double xr2 = xr[k2];  //x[k+N/2]
   167:           double xi2 = xi[k2];
   168:           double xr3 = xr[k3];  //x[k+3*N/4]
   169:           double xi3 = xi[k3];
   170:           if (false) {
   171:             xr[k0] = (xr0 + xr2) + (xr1 + xr3);  //X[4*j]=(x[k]+x[k+N/2])+(x[k+N/4]+x[k+3*N/4])
   172:             xi[k0] = (xi0 + xi2) + (xi1 + xi3);
   173:             double tr = (xr0 + xr2) - (xr1 + xr3);  //(x[k]+x[k+N/2])-(x[k+N/4]+x[k+3*N/4])
   174:             double ti = (xi0 + xi2) - (xi1 + xi3);
   175:             xr[k1] = or2 * tr - oi2 * ti;  //X[4*j+2]=ω[N]^(-2*k)*((x[k]+x[k+N/2])-(x[k+N/4]+x[k+3*N/4]))
   176:             xi[k1] = or2 * ti + oi2 * tr;
   177:             tr = (xr0 - xr2) + (xi1 - xi3);  //(x[k]-x[k+N/2])-i*(x[k+N/4]-x[k+3*N/4])
   178:             ti = (xi0 - xi2) - (xr1 - xr3);
   179:             xr[k2] = or1 * tr - oi1 * ti;  //X[4*j+1]=ω[N]^(-k)*((x[k]-x[k+N/2])-i*(x[k+N/4]-x[k+3*N/4]))
   180:             xi[k2] = or1 * ti + oi1 * tr;
   181:             tr = (xr0 - xr2) - (xi1 - xi3);  //(x[k]-x[k+N/2])+i*(x[k+N/4]-x[k+3*N/4])
   182:             ti = (xi0 - xi2) + (xr1 - xr3);
   183:             xr[k3] = or3 * tr - oi3 * ti;  //X[4*j+3]=ω[N]^(-3*k)*((x[k]-x[k+N/2])+i*(x[k+N/4]-x[k+3*N/4]))
   184:             xi[k3] = or3 * ti + oi3 * tr;
   185:           } else {
   186:             double pr, pi, qr, qi;
   187:             xr[k0] = (pr = xr0 + xr2) + (qr = xr1 + xr3);  //X[4*j]=(x[k]+x[k+N/2])+(x[k+N/4]+x[k+3*N/4])
   188:             xi[k0] = (pi = xi0 + xi2) + (qi = xi1 + xi3);
   189:             double tr = pr - qr;  //(x[k]+x[k+N/2])-(x[k+N/4]+x[k+3*N/4])
   190:             double ti = pi - qi;
   191:             xr[k1] = or2 * tr - oi2 * ti;  //X[4*j+2]=ω[N]^(-2*k)*((x[k]+x[k+N/2])-(x[k+N/4]+x[k+3*N/4]))
   192:             xi[k1] = or2 * ti + oi2 * tr;
   193:             tr = (pr = xr0 - xr2) + (qr = xi1 - xi3);  //(x[k]-x[k+N/2])-i*(x[k+N/4]-x[k+3*N/4])
   194:             ti = (pi = xi0 - xi2) - (qi = xr1 - xr3);
   195:             xr[k2] = or1 * tr - oi1 * ti;  //X[4*j+1]=ω[N]^(-k)*((x[k]-x[k+N/2])-i*(x[k+N/4]-x[k+3*N/4]))
   196:             xi[k2] = or1 * ti + oi1 * tr;
   197:             tr = pr - qr;  //(x[k]-x[k+N/2])+i*(x[k+N/4]-x[k+3*N/4])
   198:             ti = pi + qi;
   199:             xr[k3] = or3 * tr - oi3 * ti;  //X[4*j+3]=ω[N]^(-3*k)*((x[k]-x[k+N/2])+i*(x[k+N/4]-x[k+3*N/4]))
   200:             xi[k3] = or3 * ti + oi3 * tr;
   201:           }
   202:         }
   203:       }
   204:     }
   205:     if (quarter == 1) {  //N=4
   206:       for (int k0 = 0; k0 < fftN; k0 += 4) {  //0
   207:         int k1 = k0 + 1;  //1
   208:         int k2 = k0 + 2;  //2
   209:         int k3 = k0 + 3;  //3
   210:         double xr0 = xr[k0];  //x[0]
   211:         double xi0 = xi[k0];
   212:         double xr1 = xr[k1];  //x[1]
   213:         double xi1 = xi[k1];
   214:         double xr2 = xr[k2];  //x[2]
   215:         double xi2 = xi[k2];
   216:         double xr3 = xr[k3];  //x[3]
   217:         double xi3 = xi[k3];
   218:         if (false) {
   219:           xr[k0] = (xr0 + xr2) + (xr1 + xr3);  //X[0]=(x[0]+x[2])+(x[1]+x[3])
   220:           xi[k0] = (xi0 + xi2) + (xi1 + xi3);
   221:           xr[k1] = (xr0 + xr2) - (xr1 + xr3);  //X[2]=(x[0]+x[2])-(x[1]+x[3])
   222:           xi[k1] = (xi0 + xi2) - (xi1 + xi3);
   223:           xr[k2] = (xr0 - xr2) + (xi1 - xi3);  //X[1]=(x[0]-x[2])-i*(x[1]-x[3])
   224:           xi[k2] = (xi0 - xi2) - (xr1 - xr3);
   225:           xr[k3] = (xr0 - xr2) - (xi1 - xi3);  //X[3]=(x[0]-x[2])+i*(x[1]-x[3])
   226:           xi[k3] = (xi0 - xi2) + (xr1 - xr3);
   227:         } else {
   228:           double t0p2r = xr0 + xr2;
   229:           double t0p2i = xi0 + xi2;
   230:           double t0m2r = xr0 - xr2;
   231:           double t0m2i = xi0 - xi2;
   232:           double t1p3r = xr1 + xr3;
   233:           double t1p3i = xi1 + xi3;
   234:           double t1m3r = xr1 - xr3;
   235:           double t1m3i = xi1 - xi3;
   236:           xr[k0] = t0p2r + t1p3r;  //X[0]=(x[0]+x[2])+(x[1]+x[3])
   237:           xi[k0] = t0p2i + t1p3i;
   238:           xr[k1] = t0p2r - t1p3r;  //X[2]=(x[0]+x[2])-(x[1]+x[3])
   239:           xi[k1] = t0p2i - t1p3i;
   240:           xr[k2] = t0m2r + t1m3i;  //X[1]=(x[0]-x[2])-i*(x[1]-x[3])
   241:           xi[k2] = t0m2i - t1m3r;
   242:           xr[k3] = t0m2r - t1m3i;  //X[3]=(x[0]-x[2])+i*(x[1]-x[3])
   243:           xi[k3] = t0m2i + t1m3r;
   244:         }
   245:       }
   246:     } else {  //N=2
   247:       for (int k0 = 0; k0 < fftN; k0 += 2) {  //0
   248:         int k1 = k0 + 1;  //1
   249:         double xr0 = xr[k0];  //x[0]
   250:         double xi0 = xi[k0];
   251:         double xr1 = xr[k1];  //x[1]
   252:         double xi1 = xi[k1];
   253:         xr[k0] = xr0 + xr1;  //X[0]=x[0]+x[1]
   254:         xi[k0] = xi0 + xi1;
   255:         xr[k1] = xr0 - xr1;  //X[1]=x[0]-x[1]
   256:         xi[k1] = xi0 - xi1;
   257:       }
   258:     }
   259:     //インデックスのビット反転
   260:     for (int k = 0, l = fftSwap.length; k < l; k += 2) {
   261:       int i = fftSwap[k];
   262:       int j = fftSwap[k + 1];
   263:       double t = xr[i];
   264:       xr[i] = xr[j];
   265:       xr[j] = t;
   266:       t = xi[i];
   267:       xi[i] = xi[j];
   268:       xi[j] = t;
   269:     }
   270:   }  //fft.fftSandeTukey4(double[],double[])
   271: 
   272: }  //class FFT
   273: 
   274: 
   275: