MFCC算法的JAVA实现

public final class FFT { public int logm; final int MAXLOGM = 20; /* max FFT length 2^MAXLOGM */ final double TWOPI = 6.28318530717958647692; final double SQHALF = 0.707106781186547524401; int brseed[] = new int[4048]; float tab[][]; public FFT(int nlength) { double dtemp = Math.log(nlength) / Math.log(2); if ((dtemp - (int) dtemp) != 0.0) { throw new Error("FFT length must be a power of 2."); } else { this.logm = (int) dtemp; } if (logm >= 4) { creattab(logm); } } /** * Calculates the magnitude spectrum of a real signal. The returned vector * contains only the positive frequencies. */ public float[] calculateFFTMagnitude(float x[]) { int i, n; n = 1 << this.logm; if (x.length > n) { throw new Error("Tried to use a " + n + "-points FFT for a vector with " + x.length + " samples!"); } rsfft(x); float[] mag = new float[n / 2 + 1]; mag[0] = x[0]; // DC frequency must be positive always if (n == 1) { return mag; } mag[n / 2] = (float) Math.abs(x[n / 2]); // pi (meaning: fs / 2) // System.out.println("FFT before magnitude"); // IO.DisplayVector(x); for (i = 1; i < n / 2; i++) { mag[i] = (float) Math.sqrt(x[i] * x[i] + x[n - i] * x[n - i]); // System.out.println(mag[i] + " " + x[i] + " " + x[n-i]); } // IO.DisplayVector(mag); return mag; } /** * Calculates the magnitude spectrum of a real signal. The returned vector * contains only the positive frequencies. */ public double[] calculateFFTMagnitude(double inputData[]) { int i, n; n = 1 << this.logm; if (inputData.length > n) { throw new Error("Tried to use a " + n + "-points FFT for a vector with " + inputData.length + " samples!"); } // System.out.println("magnitude test"); // double[] dtest = DSP.DFTMagnitude(inputData); // IO.DisplayVector(dtest); float[] x = new float[n]; for (i = 0; i < inputData.length; i++) { x[i] = (float) inputData[i]; } rsfft(x); // System.out.println("FFT before magnitude"); // IO.DisplayVector(x); double[] mag = new double[n / 2 + 1]; mag[0] = x[0]; // DC frequency must be positive always if (n == 1) { return mag; } mag[n / 2] = Math.abs(x[n / 2]); // pi (meaning: fs / 2) for (i = 1; i < n / 2; i++) { mag[i] = Math.sqrt(x[i] * x[i] + x[n - i] * x[n - i]); // System.out.println(mag[i] + " " + x[i] + " " + x[n-i]); } // IO.DisplayVector(mag); return mag; } /** * Calculates the power (magnitude squared) spectrum of a real signal. The * returned vector contains only the positive frequencies. */ public double[] calculateFFTPower(double inputData[]) { int i, n; n = 1 << this.logm; // System.out.println("magnitude test"); // double[] dtest = DSP.DFTMagnitude(inputData); // IO.DisplayVector(dtest); float[] x = new float[n]; for (i = 0; i < inputData.length; i++) { x[i] = (float) inputData[i]; } rsfft(x); // System.out.println("FFT before magnitude"); // IO.DisplayVector(x); double[] mag = new double[n / 2 + 1]; mag[0] = x[0]; // DC frequency must be positive always if (n == 1) { return mag; } mag[n / 2] = Math.abs(x[n / 2]); // pi (meaning: fs / 2) for (i = 1; i < n / 2; i++) { mag[i] = x[i] * x[i] + x[n - i] * x[n - i]; // mag[i] = Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]); // System.out.println(mag[i] + " " + x[i] + " " + x[n-i]); } // IO.DisplayVector(mag); return mag; } /** * In place calculation of FFT magnitude. */ public void FFTMagnitude(float x[]) { int i, n; rsfft(x); n = 1 << this.logm; if (n == 1) return; for (i = 1; i < n / 2; i++) { x[i] = (float) Math.sqrt(x[i] * x[i] + x[n - i] * x[n - i]); x[n - i] = x[i]; } // Aldebaro modification: x[n / 2] = Math.abs(x[n / 2]); } void rsfft(float x[]) { /* creat table */ // if(logm>=4) creattab(logm); /* Call recursive routine */ rsrec(x, logm); /* Output array unshuffling using bit-reversed indices */ if (logm > 1) { BR_permute(x, logm); return; } } /* * -------------------------------------------------------------------- * * Inverse transform for real inputs * * -------------------------------------------------------------------- */ void rsifft(float x[]) { int i, m; float fac; int xp; /* Output array unshuffling using bit-reversed indices */ if (logm > 1) { BR_permute(x, logm); } x[0] *= 0.5; if (logm > 0) x[1] *= 0.5; /* creat table */ // if(logm>=4) creattab(logm); /* Call recursive routine */ rsirec(x, logm); /* Normalization */ m = 1 << logm; fac = (float) 2.0 / m; xp = 0; for (i = 0; i < m; i++) { x[xp++] *= fac; } } /* * -------------------------------------------------------------------- * * Creat multiple fator table * * -------------------------------------------------------------------- */ void creattab(int logm) { int m, m2, m4, m8, nel, n, rlogm; int cn, spcn, smcn, c3n, spc3n, smc3n; double ang, s, c; tab = new float[logm - 4 + 1][6 * ((1 << logm) / 4 - 2)]; for (rlogm = logm; rlogm >= 4; rlogm--) { m = 1 << rlogm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2; nel = m4 - 2; /* Initialize pointers */ cn = 0; spcn = cn + nel; smcn = spcn + nel; c3n = smcn + nel; spc3n = c3n + nel; smc3n = spc3n + nel; /* Compute tables */ for (n = 1; n < m4; n++) { if (n == m8) continue; ang = n * TWOPI / m; c = Math.cos(ang); s = Math.sin(ang); tab[rlogm - 4][cn++] = (float) c; tab[rlogm - 4][spcn++] = (float) (-(s + c)); tab[rlogm - 4][smcn++] = (float) (s - c); ang = 3 * n * TWOPI / m; c = Math.cos(ang); s = Math.sin(ang); tab[rlogm - 4][c3n++] = (float) c; tab[rlogm - 4][spc3n++] = (float) (-(s + c)); tab[rlogm - 4][smc3n++] = (float) (s - c); } } } /* * -------------------------------------------------------------------- * * Recursive part of the RSFFT algorithm. Not externally * callable. * * -------------------------------------------------------------------- */ void rsrec(float x[], int logm) { int m, m2, m4, m8, nel, n; int x0 = 0; int xr1, xr2, xi1; int cn = 0; int spcn = 0; int smcn = 0; float tmp1, tmp2; double ang, c, s; /* Check range of logm */ try { if ((logm < 0) || (logm > MAXLOGM)) { System.err.println("FFT length m is too big: log2(m) = " + logm + "is out of bounds [" + 0 + "," + MAXLOGM + "]"); throw new OutofborderException(logm); } } catch (OutofborderException e) { throw new OutOfMemoryError(); } /* Compute trivial cases */ if (logm < 2) { if (logm == 1) { /* length m = 2 */ xr2 = x0 + 1; tmp1 = x[x0] + x[xr2]; x[xr2] = x[x0] - x[xr2]; x[x0] = tmp1; return; } else if (logm == 0) return; /* length m = 1 */ } /* Compute a few constants */ m = 1 << logm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2; /* Build tables of butterfly coefficients, if necessary */ // if ((logm >= 4) && (tab[logm-4][0] == 0)) { /* Allocate memory for tables */ // nel = m4 - 2; /* * if ((tab[logm-4] = (float *) calloc(3 * nel, sizeof(float))) == NULL) * { printf("Error : RSFFT : not enough memory for cosine tables.\n"); * error_exit(); } */ /* Initialize pointers */ // tabi=logm-4; // cn =0; spcn = cn + nel; smcn = spcn + nel; /* Compute tables */ /* * for (n = 1; n < m4; n++) { if (n == m8) continue; ang = n * * (float)TWOPI / m; c = Math.cos(ang); s = Math.sin(ang); * tab[tabi][cn++] = (float)c; tab[tabi][spcn++] = (float)(- (s + c)); * tab[tabi][smcn++] =(float)( s - c); } } *