Hi all.

I'm very new to the DSP world, and I've been spending
some time puzzling over the Fast Fourier Transform. I hope you can
entertain some questions and let me know if I'm on the
right track.

I have based the majority of my understanding on the
book Numerical Recipes in C, 2nd ED, Chapter 12, and the
four() algorithm listed therein.

For other beginners I recommend:
Fourier Transforms and the Fast Fourier Transform (FFT) Algorithm
by Paul Heckbert, Feb. 1995 and
Dave's Short Course on Complex Numbers
http://www.clarku.edu/~djoyce/complex/

I decided to try and create a Discrete
Fourier Transform algorithm, as a way of explaining
to myself exactly what's going on. Here it is.
Comments welcome.

(NOTE that the algorithm is not designed to be fast,
but to be as obvious as possible. I understand that the
mathematics is very redundant and operations could be
reduced via Danielson/Lanczos divide & conquer. Also,
some of the complex math can be taken out because I only
use the real components of the input/output buffers.
Also, my implimentation uses library trigonometry, which
could be sped up by lookup tables)

-Andrew Schwerin


//AOS_DFT illustration
//Andrew Schwerin, March 2003
//
//Discrete Fourier Transform without speedups.
//Let me know if I get something wrong.
//
//see also four1() from Numerical Recipes in C, 2nd Ed., page 507-508
//
//First off, N is the number of data points for input and output, 8 in this example.
//It does not need to be (but is anyway) a power of 2, because we do not
//use the Danielson/Lanczos Lemma of chopping the problem into
//smaller pieces.
//Second, I prefer the input to be buffered into real and
//imaginary arrays simply to clarify some of the pointer arithmetic.
//
//Third, I have tried to express the general problem in terms of trigonometry.
//
//Start with the textbook definition:
//
// W = e^((2*PI*i)/N)
//
// H[n] = SUM from k=0 to k=N-1 of (W^(nk))*h[k]
//
// Since e^(i*angle) = cos (angle) + i*sin(angle)
//
// W is a constant and can be interpreted as the angle 2*PI/N
// which is, essentially, dividing the circle (2PI radians) into 1/Nth the size.
// This means that any angle that moves by W, will repeat itself every N times.
// This angle is used to generate basis sin and cos with which the input data is correlated.
// In DSP, to detect correlation of two signals you multiply the signals together
// and add up each of the multiplications to arrive at a correlation number.
//
// To use some clumsy vernacular, it is very much like moving a wheel along the data.
// That is, you say, if I had a wheel spinning with this so and so angular velocity,
// And I ran it across the input data, what would the total correlation be?
//
// H[n] = SUM from k=0 to k=N-1 of TRIGONOMETRY(W*n*k) * h[k]
//
// H[n] = SUM correlation of SPINNING_WHEEL * input_signal, over entire span of input
// and for any n, the wheel spins at W*n angles/sample
//
// QUESTIONS FOR THE AUDIENCE:
// Can you get more data by asking for more frequencies in the answer,
// like with non-integer n n=0,0.2,0.4,0.6...
//
// Can the data for n >= N/2 be safely ignored? In my model it represents
// a wheel spinning so fast that it seems to spin backwards.
//
// Is there any more acurracy to be had by interpolating along the input data
// by non-integer k as in t[0]=t[0],t[0.5] = (t[0]+t[1])/2
//
//

#include <math.h>

#define DATA_SIZE 8 //Number of input data points in this example
#define TWO_PI 6.28318530717959

void AOS_DFT(float data_in__r[],
float data_in_i[],
float data_out_r[],
float data_out_i[],
unsigned long N);//length of the arrays

int main()
{
printf("Illustration of AOS_DFT\n\n");

float input_r[DATA_SIZE];//real
float input_i[DATA_SIZE];//imaginary
float output_r[DATA_SIZE];//real
float output_i[DATA_SIZE];//imaginary

/* load data */
input_r[0] = 1.0;input_i[0] = 0.0;
input_r[1] = 0.0;input_i[1] = 0.0;
input_r[2] = -7.0;input_i[2] = 0.0;
input_r[3] = 0.0;input_i[3] = 0.0;
input_r[4] = 1.0;input_i[4] = 0.0;
input_r[5] = 0.0;input_i[5] = 0.0;
input_r[6] = -1.0;input_i[6] = 0.0;
input_r[7] = 0.0;input_i[7] = 0.0;

/* do discrete fourier transform */
AOS_DFT(input_r, input_i, output_r, output_i, DATA_SIZE);

/* display results */
for(int k=0;k<DATA_SIZE;k++)
{ printf("output[%d]: {%f,%f}\n", k, output_r[k],output_i[k]); }

printf("\ndone.\n");
return 0;
}

void AOS_DFT(float data_in_r[],
float data_in_i[],
float data_out_r[],
float data_out_i[],
unsigned long N)
{
float delta_angle = TWO_PI/N;
printf("Compute e^((2*PI*i)/N) as an angle: %f degrees\n", (delta_angle*360)/TWO_PI);

for (unsigned long Hn=0;Hn<N;Hn++)
{
//COMPUTE H[Hn], the correllation for some given frequency

float SUM_r,SUM_i;
SUM_r = 0.0;
SUM_i = 0.0;

float WHEEL_SPEED = Hn * delta_angle;//frequency of basis function. Roughly like W^Hn
float WHEEL_ANGLE = 0; //all basis functions start at phase of 0 angles

for(int k=0;k<N;k++)
{
//Do corellation with input data h[k]

float SAMPLE_r,SAMPLE_i;//Get the h[k]
SAMPLE_r = data_in_r[k];
SAMPLE_i = data_in_i[k];
float WHEEL_r,WHEEL_i;//Calculate the basis function produced by cos(W^n^k) and sin(W^n^k)
WHEEL_r = cos(WHEEL_ANGLE);
WHEEL_i = sin(WHEEL_ANGLE);
float GROOVE_r,GROOVE_i;//do the DSP multiply with complex numbers
GROOVE_r = (SAMPLE_r * WHEEL_r) - (SAMPLE_i * WHEEL_i);
GROOVE_i = (SAMPLE_r * WHEEL_i) + (SAMPLE_i * WHEEL_r);

SUM_r += GROOVE_r;//Add the correllation into the result buffer
SUM_i += GROOVE_i;

float display = (WHEEL_ANGLE*360)/TWO_PI;while(display>=360){display-=360;}

printf("Hn=%d k=%d wheel_angle=%.1f GROOVINESS: (%.2f,%.2f) = (%.2f,%.2f)*(%.2f,%.2f)\n", Hn, k, display, GROOVE_r, GROOVE_i, SAMPLE_r, SAMPLE_i, WHEEL_r, WHEEL_i);

WHEEL_ANGLE += WHEEL_SPEED;//move the basis function forward in time by turning the wheel
}

data_out_r[Hn] = SUM_r;//Put computed H[n] into output buffer
data_out_i[Hn] = SUM_i;
printf(" SOLVED H[%d]: (%f,%f)\n", Hn, SUM_r, SUM_i);
}
}