Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
152 lines (129 sloc) 3.86 KB
/* Copyright (c) 2016 Nick Appleton
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE. */
/* A boring FFT as used in an example on www.appletonaudio.com */
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
void fft_cplx(float *inout, float *scratch, unsigned len)
{
unsigned i;
if (len == 1)
return;
/* DIF. */
for (i = 0; i < len / 2; i++) {
float re0 = inout[2*i+0];
float im0 = inout[2*i+1];
float re1 = inout[2*i+len+0];
float im1 = inout[2*i+len+1];
float twr = cosf(i * -2.0f * M_PI / len);
float twi = sinf(i * -2.0f * M_PI / len);
float sr = re0 - re1;
float si = im0 - im1;
scratch[2*i+0] = re0 + re1;
scratch[2*i+1] = im0 + im1;
scratch[2*i+len+0] = sr * twr - si * twi;
scratch[2*i+len+1] = sr * twi + si * twr;
}
/* recurse. */
fft_cplx(scratch, inout, len/2);
fft_cplx(scratch+len, inout, len/2);
/* reorder. */
for (i = 0; i < len / 2; i++) {
inout[4*i+0] = scratch[2*i+0];
inout[4*i+1] = scratch[2*i+1];
inout[4*i+2] = scratch[2*i+len+0];
inout[4*i+3] = scratch[2*i+len+1];
}
}
void fft_real_sft(float *inout, float *scratch, unsigned len)
{
unsigned i;
if (len == 1)
return;
/* data munge. */
for (i = 0; i < len/2; i++) {
float re0 = inout[i];
float im0 = -inout[i+len/2];
float twr = cos(-M_PI * i * 1.0f / len);
float twi = sin(-M_PI * i * 1.0f / len);
scratch[2*i+0] = re0 * twr - im0 * twi;
scratch[2*i+1] = re0 * twi + im0 * twr;
}
/* DFT. */
fft_cplx(scratch, inout, len/2);
/* reorder. */
if (len == 2) {
inout[0] = scratch[0];
inout[1] = scratch[1];
} else {
for (i = 0; i < len/4; i++) {
inout[4*i+0] = scratch[2*i+0];
inout[4*i+1] = scratch[2*i+1];
inout[4*i+2] = scratch[len-2*i-2];
inout[4*i+3] = -scratch[len-2*i-1];
}
}
}
void fft_real(float *inout, float *scratch, unsigned len)
{
unsigned i;
if (len == 1)
return;
/* twiddle. */
for (i = 0; i < len/2; i++) {
float re0 = inout[i];
float re1 = inout[i+len/2];
scratch[i] = re0 + re1;
scratch[i+len/2] = re0 - re1;
}
/* recurse. */
fft_real (scratch, inout, len/2);
fft_real_sft(scratch+len/2, inout+len/2, len/2);
/* reorder. */
if (len == 2) {
inout[0] = scratch[0];
inout[1] = scratch[1];
} else {
for (i = 0; i < len/4; i++) {
float re0 = scratch[2*i+0];
float im0 = scratch[2*i+1];
float re1 = scratch[2*i+len/2+0];
float im1 = scratch[2*i+len/2+1];
inout[4*i+0] = re0;
inout[4*i+1] = im0;
inout[4*i+2] = re1;
inout[4*i+3] = im1;
}
}
}
int main(int argc, char *argv[])
{
float io[32];
float scratch[32];
unsigned i;
for (i = 0; i < 32; i++) {
io[i] = sin(i*M_PI*3*2.0/32 + M_PI/4);
}
fft_real(io, scratch, 32);
for (i = 0; i < 16; i++) {
printf("%f,%f\n", io[2*i+0], io[2*i+1]);
}
}