DSP tutorial: Calculating FFT of a sine wave
Last modified: November 16th, 2011
This example generates a sine wave, calculates the FFT transform of it and determines the sine wave’s frequency by the FFT transformed values. FFT library used is JTransforms.
FFT is: Fast Fourier transform. It’s an efficient algorithm to compute the discrete Fourier transform (DFT) and its inverse. Oh, and what that’s used for? We can get the relative strengths of the periodic components of the signal to be transformed. So, if we have a buffer filled with sound data and we calculate the DFT of it (using the FFT algorithm), we will get the relative power of each frequency in that buffer.
We can calculate the inverse DFT also. If we calculate the DFT of a signal and then inverse DFT it, we get back the original signal. If we modify the transformed signal and then we calculate the inverse, we will get a different signal. As the transformed signal data represents the power of the frequency components, if we zero them out and calculate the inverse DFT, we can make a filter (zeroing out the unwanted frequencies).
Important note: when we calculate a signal’s DFT, we get complex numbers, so we have to calculate with those. If we want the strength of a frequency component, we can get it by calculating the “length” of the complex number represented as a vector (sqrt(x^2+y^2)).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 | public class SineFFT { private final static int SAMPLERATE = 8000; public static void main(String[] args) { int frequency = 3560; // freq of our sine wave double lengthInSecs = 1; int samplesNum = (int) Math.round(lengthInSecs * SAMPLERATE); System.out.println("Samplesnum: " + samplesNum); double[] audioData = new double[samplesNum]; int samplePos = 0; // http://en.wikibooks.org/wiki/Sound_Synthesis_Theory/Oscillators_and_Wavetables for (double phase = 0; samplePos < lengthInSecs * SAMPLERATE && samplePos < samplesNum; phase += (2 * Math.PI * frequency) / SAMPLERATE) { audioData[samplePos++] = Math.sin(phase); if (phase >= 2 * Math.PI) phase -= 2 * Math.PI; } // we compute the fft of the whole sine wave DoubleFFT_1D fft = new DoubleFFT_1D(samplesNum); // we need to initialize a buffer where we store our samples as complex numbers. first value is the real part, second is the imaginary. double[] fftData = new double[samplesNum * 2]; for (int i = 0; i < samplesNum; i++) { // copying audio data to the fft data buffer, imaginary part is 0 fftData[2 * i] = audioData[i]; fftData[2 * i + 1] = 0; } // calculating the fft of the data, so we will have spectral power of each frequency component // fft resolution (number of bins) is samplesNum, because we initialized with that value fft.complexForward(fftData); try { // writing the values to a txt file BufferedWriter outputStream = new BufferedWriter(new OutputStreamWriter(new FileOutputStream("output.txt"), "UTF-8")); int max_i = -1; double max_fftval = -1; for (int i = 0; i < fftData.length; i += 2) { // we are only looking at the half of the spectrum double hz = ((i / 2.0) / fftData.length) * SAMPLERATE; outputStream.write(i + ".\tr:" + Double.toString((Math.abs(fftData[i]) > 0.1 ? fftData[i] : 0)) + " i:" + Double.toString((Math.abs(fftData[i + 1]) > 0.1 ? fftData[i + 1] : 0)) + "\t\t" + hz + "hz\n"); // complex numbers -> vectors, so we compute the length of the vector, which is sqrt(realpart^2+imaginarypart^2) double vlen = Math.sqrt(fftData[i] * fftData[i] + fftData[i + 1] * fftData[i + 1]); if (max_fftval < vlen) { // if this length is bigger than our stored biggest length max_fftval = vlen; max_i = i; } } double dominantFreq = ((max_i / 2.0) / fftData.length) * SAMPLERATE; System.out.println("Dominant frequency: " + dominantFreq + "hz (output.txt line no. " + max_i + ")"); outputStream.close(); } catch (Exception e) { e.printStackTrace(); } } } |
Projects
I just run your code snippet and don’t unterstand this:
int frequency = 3560; // freq of our sine wave
output: 7120. r:0.0 i:-4000.000000000001 1780.0hz
The output frequency is half of the original frequency. Why?
The problems are on lines 44 and 57. Since fftData is already (samplesNum * 2) length (due to the Real/Imaginary-pairs layout required by JTransforms – check http://incanter.org/docs/parallelcolt/api/edu/emory/mathcs/jtransforms/fft/DoubleFFT_1D.html#complexForward%28double%5B%5D%29 ) and the for-loop is already incrementing i by 2, there is no need to divide i by 2 on said lines (each pair of elements on fftData represent a single frequency).
So, to fix it, cast i and max_i to double on lines 44 or 57, and remove the division by 2. Both the output frequency and the max frequency read (which should be close to SAMPLERATE value of 8000, and which was also halved) will be correct.