/* * FFT.java * 用于频谱显示的快速傅里叶变换 * http://jmp123.sf.net/ */ public class FFT { public static final int FFT_N_LOG = 9; // FFT_N_LOG <= 13 public static final int FFT_N = 1 << FFT_N_LOG; private static final float MINY = (float) ((FFT_N << 2) * Math.sqrt(2)); //(*) private float[] real, imag, sintable, costable; private int[] bitReverse; public FFT() { real = new float[FFT_N]; imag = new float[FFT_N]; sintable = new float[FFT_N >> 1]; costable = new float[FFT_N >> 1]; bitReverse = new int[FFT_N]; int i, j, k, reve; for (i = 0; i < FFT_N; i++) { k = i; for (j = 0, reve = 0; j != FFT_N_LOG; j++) { reve <<= 1; reve |= (k & 1); k >>>= 1; } bitReverse[i] = reve; } double theta, dt = 2 * 3.14159265358979323846 / FFT_N; for (i = 0; i < (FFT_N >> 1); i++) { theta = i * dt; costable[i] = (float) Math.cos(theta); sintable[i] = (float) Math.sin(theta); } } /** * 用于频谱显示的快速傅里叶变换 * @param realIO 输入FFT_N个实数,也用它暂存fft后的FFT_N/2个输出值(复数模的平方)。 */ public void calculate(float[] realIO) { int i, j, k, ir, exchanges = 1, idx = FFT_N_LOG - 1; float cosv, sinv, tmpr, tmpi; for (i = 0; i != FFT_N; i++) { real[i] = realIO[bitReverse[i]]; imag[i] = 0; } for (i = FFT_N_LOG; i != 0; i--) { for (j = 0; j != exchanges; j++) { cosv = costable[j << idx]; sinv = sintable[j << idx]; for (k = j; k < FFT_N; k += exchanges << 1) { ir = k + exchanges; tmpr = cosv * real[ir] - sinv * imag[ir]; tmpi = cosv * imag[ir] + sinv * real[ir]; real[ir] = real[k] - tmpr; imag[ir] = imag[k] - tmpi; real[k] += tmpr; imag[k] += tmpi; } } exchanges <<= 1; idx--; } j = FFT_N >> 1; /* * 输出模的平方(的FFT_N倍): * for(i = 1; i <= j; i++) * realIO[i-1] = real[i] * real[i] + imag[i] * imag[i]; * * 如果FFT只用于频谱显示,可以"淘汰"幅值较小的而减少浮点乘法运算. MINY的值 * 和Spectrum.Y0,Spectrum.logY0对应. */ sinv = MINY; cosv = -MINY; for (i = j; i != 0; i--) { tmpr = real[i]; tmpi = imag[i]; if (tmpr > cosv && tmpr < sinv && tmpi > cosv && tmpi < sinv) realIO[i - 1] = 0; else realIO[i - 1] = tmpr * tmpr + tmpi * tmpi; } } }
时间: 2024-10-06 04:43:16