using System;using System.Diagnostics;using System.Numerics;namespace Fft
{classProgram{staticvoidWarmup(){var random =newRandom();int N =1048576;
Complex[] x =newComplex[N];// original datafor(int i =0; i < N; i++){
x[i]=newComplex(i,0);
x[i]=newComplex(-2* random.Next()+1,0);}// FFT of original data
Complex[] y = FFTImpl.fft(x);// take inverse FFT
Complex[] z = FFTImpl.ifft(y);// circular convolution of x with itself
Complex[] c = FFTImpl.cconvolve(x, x);// linear convolution of x with itself
Complex[] d = FFTImpl.convolve(x, x);}staticvoidMain(string[] args){for(int i =0; i <10; i++)Warmup();// warming upfor(int j =0; j <5; j++){var st =newStopwatch();
st.Start();var random =newRandom();int N =1048576;
Complex[] x =newComplex[N];// original datafor(int i =0; i < N; i++){
x[i]=newComplex(i,0);
x[i]=newComplex(-2* random.Next()+1,0);}// FFT of original data
Complex[] y = FFTImpl.fft(x);// take inverse FFT
Complex[] z = FFTImpl.ifft(y);// circular convolution of x with itself
Complex[] c = FFTImpl.cconvolve(x, x);// linear convolution of x with itself
Complex[] d = FFTImpl.convolve(x, x);
st.Stop();
Console.WriteLine(st.ElapsedMilliseconds);}}}}
FFTImpl.cs
using System;using System.Numerics;namespace Fft
{publicclassFFTImpl{// compute the FFT of x[], assuming its length is a power of 2publicstatic Complex[]fft(Complex[] x){int N = x.Length;// base caseif(N ==1)returnnewComplex[]{ x[0]};// radix 2 Cooley-Tukey FFTif(N %2!=0){thrownewException("N is not a power of 2");}// fft of even terms
Complex[] even =newComplex[N/2];for(int k =0; k < N /2; k++){
even[k]= x[2* k];}
Complex[] q = FFTImpl.fft(even);// fft of odd terms
Complex[] odd = even;// reuse the arrayfor(int k =0; k < N /2; k++){
odd[k]= x[2* k +1];}
Complex[] r = FFTImpl.fft(odd);// combine
Complex[] y =newComplex[N];for(int k =0; k < N /2; k++){double kth =-2* k * Math.PI / N;Complex wk =newComplex(Math.Cos(kth), Math.Sin(kth));
y[k]= q[k]+(wk * r[k]);
y[k + N /2]= q[k]-(wk * r[k]);}return y;}// compute the inverse FFT of x[], assuming its length is a power of 2publicstatic Complex[]ifft(Complex[] x){int N = x.Length;
Complex[] y =newComplex[N];// take conjugatefor(int i =0; i < N; i++){
y[i]= Complex.Conjugate(x[i]);}// compute forward FFT
y = FFTImpl.fft(y);// take conjugate againfor(int i =0; i < N; i++){
y[i]= Complex.Conjugate(y[i]);}// divide by Nfor(int i =0; i < N; i++){
y[i]= y[i]*newComplex(1.0/ N,0);}return y;}// compute the circular convolution of x and ypublicstatic Complex[]cconvolve(Complex[] x, Complex[] y){// should probably pad x and y with 0s so that they have same length// and are powers of 2if(x.Length != y.Length){thrownewException("Dimensions don't agree");}int N = x.Length;// compute FFT of each sequence
Complex[] a = FFTImpl.fft(x);
Complex[] b = FFTImpl.fft(y);// point-wise multiply
Complex[] c =newComplex[N];for(int i =0; i < N; i++){
c[i]= a[i]* b[i];}// compute inverse FFTreturn FFTImpl.ifft(c);}// compute the linear convolution of x and ypublicstatic Complex[]convolve(Complex[] x, Complex[] y){Complex ZERO =newComplex(0,0);
Complex[] a =newComplex[2* x.Length];for(int i =0; i < x.Length; i++) a[i]= x[i];for(int i = x.Length; i <2* x.Length; i++) a[i]= ZERO;
Complex[] b =newComplex[2* y.Length];for(int i =0; i < y.Length; i++) b[i]= y[i];for(int i = y.Length; i <2* y.Length; i++) b[i]= ZERO;return FFTImpl.cconvolve(a, b);}}}
Java(复数类使用 Complex.java)
Main.java:
package com.company;publicclassMain{publicstaticvoidwarmup(){int N =1048576;
Complex[] x =newComplex[N];// original datafor(int i =0; i < N; i++){
x[i]=newComplex(i,0);
x[i]=newComplex(-2*Math.random()+1,0);}// FFT of original data
Complex[] y = FFT.fft(x);// take inverse FFT
Complex[] z = FFT.ifft(y);// circular convolution of x with itself
Complex[] c = FFT.cconvolve(x, x);// linear convolution of x with itself
Complex[] d = FFT.convolve(x, x);}publicstaticvoidmain(String[] args){for(int i =0; i <10; i++)warmup();// warming upfor(int j =0; j <5; j++){long begintime = System.currentTimeMillis();int N =1048576;
Complex[] x =newComplex[N];// original datafor(int i =0; i < N; i++){
x[i]=newComplex(i,0);
x[i]=newComplex(-2*Math.random()+1,0);}// FFT of original data
Complex[] y = FFT.fft(x);// take inverse FFT
Complex[] z = FFT.ifft(y);// circular convolution of x with itself
Complex[] c = FFT.cconvolve(x, x);// linear convolution of x with itself
Complex[] d = FFT.convolve(x, x);long endtime = System.currentTimeMillis();
System.out.println(endtime - begintime);}}}
FFT.java:
package com.company;publicclassFFT{// compute the FFT of x[], assuming its length is a power of 2publicstatic Complex[]fft(Complex[] x){int N = x.length;// base caseif(N ==1)returnnewComplex[]{x[0]};// radix 2 Cooley-Tukey FFTif(N %2!=0){thrownewRuntimeException("N is not a power of 2");}// fft of even terms
Complex[] even =newComplex[N /2];for(int k =0; k < N /2; k++){
even[k]= x[2* k];}
Complex[] q =fft(even);// fft of odd terms
Complex[] odd = even;// reuse the arrayfor(int k =0; k < N /2; k++){
odd[k]= x[2* k +1];}
Complex[] r =fft(odd);// combine
Complex[] y =newComplex[N];for(int k =0; k < N /2; k++){double kth =-2* k * Math.PI / N;
Complex wk =newComplex(Math.cos(kth), Math.sin(kth));
y[k]= q[k].add(wk.mult(r[k]));
y[k + N /2]= q[k].minus(wk.mult(r[k]));}return y;}// compute the inverse FFT of x[], assuming its length is a power of 2publicstatic Complex[]ifft(Complex[] x){int N = x.length;
Complex[] y =newComplex[N];// take conjugatefor(int i =0; i < N; i++){
y[i]= x[i].conj();}// compute forward FFT
y =fft(y);// take conjugate againfor(int i =0; i < N; i++){
y[i]= y[i].conj();}// divide by Nfor(int i =0; i < N; i++){
y[i]= y[i].mult(newComplex(1.0/ N,0));}return y;}// compute the circular convolution of x and ypublicstatic Complex[]cconvolve(Complex[] x, Complex[] y){// should probably pad x and y with 0s so that they have same length// and are powers of 2if(x.length != y.length){thrownewRuntimeException("Dimensions don't agree");}int N = x.length;// compute FFT of each sequence
Complex[] a =fft(x);
Complex[] b =fft(y);// point-wise multiply
Complex[] c =newComplex[N];for(int i =0; i < N; i++){
c[i]= a[i].mult(b[i]);}// compute inverse FFTreturnifft(c);}// compute the linear convolution of x and ypublicstatic Complex[]convolve(Complex[] x, Complex[] y){
Complex ZERO =newComplex(0,0);
Complex[] a =newComplex[2* x.length];for(int i =0; i < x.length; i++) a[i]= x[i];for(int i = x.length; i <2* x.length; i++) a[i]= ZERO;
Complex[] b =newComplex[2* y.length];for(int i =0; i < y.length; i++) b[i]= y[i];for(int i = y.length; i <2* y.length; i++) b[i]= ZERO;returncconvolve(a, b);}}