fft.rar_fft_fft1D java_image_java FFT1D_三维 图像处理

package toolbox; import java.io.PrintStream; public class FFT1D { private boolean radix2; private double Rearg[]; private double Imarg[]; private double yReOut[]; private double yImOut[]; private int maxPrimeFactor; private int maxPrimeFactorDiv2; private int maxFactorCount; private int n; private int nFactor; private final double c3_1 = -1.5D; private final double c3_2 = 0.86602540378444004D; private final double c5_1 = -1.25D; private final double c5_2 = 0.55901699437495D; private final double c5_3 = -0.95105651629514998D; private final double c5_4 = -1.5388417685875999D; private final double c5_5 = 0.36327126400268001D; private final double c8 = 0.70710678118655002D; private int groupOffset; private int dataOffset; private int adr; private int groupNo; private int dataNo; private int blockNo; private int twNo; private double omega; private double tw_re; private double tw_im; private double twiddleRe[]; private double twiddleIm[]; private double trigRe[]; private double trigIm[]; private double zRe[]; private double zIm[]; private double vRe[]; private double vIm[]; private double wRe[]; private double wIm[]; private int sofarRadix[]; private int actualRadix[]; private int remainRadix[]; public FFT1D(int size) { radix2 = true; maxFactorCount = 20; int m = 1; for(int size1 = size; size1 > 2;) { size1 /= 2; m++; } if((int)Math.round(Math.pow(2D, m)) == size) //round()方法取最靠近的其的数值 { radix2 = true; n = 1 << m; double fact = 6.2831853071795862D / (double)n; Imarg = new double[n]; Rearg = new double[n]; for(int i = 0; i < n; i++) { double arg = fact * (double)i; Rearg[i] = Math.cos(arg); Imarg[i] = -Math.sin(arg); } } else { radix2 = false; maxPrimeFactor = (int)(double)(size + 1); maxPrimeFactorDiv2 = (maxPrimeFactor + 1) / 2; twiddleRe = new double[maxPrimeFactor]; twiddleIm = new double[maxPrimeFactor]; trigRe = new double[maxPrimeFactor]; trigIm = new double[maxPrimeFactor]; zRe = new double[maxPrimeFactor]; zIm = new double[maxPrimeFactor]; vRe = new double[maxPrimeFactorDiv2]; vIm = new double[maxPrimeFactorDiv2]; wRe = new double[maxPrimeFactorDiv2]; wIm = new double[maxPrimeFactorDiv2]; yReOut = new double[size]; yImOut = new double[size]; n = size; sofarRadix = new int[maxFactorCount]; actualRadix = new int[maxFactorCount]; remainRadix = new int[maxFactorCount]; transTableSetup(sofarRadix, actualRadix, remainRadix); } } public final void transform(double Re[], double Im[], int size, int shift) { n = size; if(radix2) { doFFT1D_CooleyTukey(Re, Im, size, shift); } else if(shift == 0) { doFFT_Mix(Re, Im, size); } else { doFFT_Mix(Re, Im, size, shift); } } public final void inverse(double Re[], double Im[], int size, int shift) { n = size; if(radix2) { doIFFT1D_CooleyTukey(Re, Im, size, shift); } else if(shift == 0) { doIFFT_Mix(Re, Im, size); } else { doIFFT_Mix(Re, Im, size, shift); } } private void doFFT1D_CooleyTukey(double Re[], double Im[], int size, int shift) { int m = 1; for(int size1 = size; size1 > 2;) { size1 /= 2; m++; } int j; for(int i = j = shift; i < (shift + n) - 1; i++) { if(i < j) { double Retmp = Re[i]; double Imtmp = Im[i]; Re[i] = Re[j]; Im[i] = Im[j]; Re[j] = Retmp; Im[j] = Imtmp; } int k; for(k = n >> 1; k + shift <= j; k /= 2) { j -= k; } j += k; } int stepsize = 1; for(int shifter = m - 1; stepsize < n; shifter--) { for(int j = shift; j < shift + n; j += stepsize << 1) { for(int i = 0; i < stepsize; i++) { int i_j = i + j; int i_j_s = i_j + stepsize; double Retmp; if(i > 0) { Retmp = Rearg[i << shifter] * Re[i_j_s] - Imarg[i << shifter] * Im[i_j_s]; Im[i_j_s] = Rearg[i << shifter] * Im[i_j_s] + Imarg[i << shifter] * Re[i_j_s]; Re[i_j_s] = Retmp; } Retmp = Re[i_j] - Re[i_j_s]; double Imtmp = Im[i_j] - Im[i_j_s]; Re[i_j] += Re[i_j_s]; Im[i_j] += Im[i_j_s]; Re[i_j_s] = Retmp; Im[i_j_s] = Imtmp; } } stepsize <<= 1; } } private void doIFFT1D_CooleyTukey(double Re[], double Im[], int size, int shift) { for(int i = shift; i < shift + size; i++) { Im[i] = -Im[i]; } transform(Re, Im, size, shift); for(int i = shift; i < shift + size; i++) { Re[i] = Re[i] / (double)size; Im[i] = -Im[i] / (double)size; } } private void factorize(int fact[], int num) { int radices[] = new int[7]; int factors[] = new int[maxFactorCount]; int nRadix = 6; radices[1] = 2; radices[2] = 3; radices[3] = 4; radices[4] = 5; radices[5] = 8; radices[6] = 10; int j; if(num == 1) { j = 1; factors[1] = 1; } else { j = 0; } for(int i = nRadix; num > 1 && i > 0;) { if(num % radices[i] == 0) { num /= radices[i]; j++; factors[j] = radices[i]; } else { i--; } } if(factors[j] == 2) { int i; for(i = j - 1; i > 0 && factors[i] != 8; i--) { } if(i > 0) { factors[j] = 4; factors[i] = 4; } } if(num > 1) { for(int k = 2; (double)k < Math.sqrt(num) + 1.0D; k++) { while(num % k == 0) { num /= k; j++; factors[j] = k; } } if(num > 1) { j++; factors[j] = num; } } for(int i = 1; i <= j; i++) { fact[i] = factors[(j - i) + 1]; } nFactor = j; } private final void transTableSetup(int sofar[], int actual[], int remain[]) { factorize(actual, n); if(actual[1] > maxPrimeFactor) { System.out.println((new StringBuilder("\nPrime factor of FFT length too large : %6d")).append(actual[1]).toString());