Hi Experts. I am Physics student, and i wanna write referat about Fourier transformations, also about FFT 1D real case. Hope this is best place for ask this, and here are best experts. I wanaa demonstrate how diverse window functions changes measured spectrum, how much CPU ressources take diverse FFT algorithms ... Yet i understand [i hope so] how works windowing . Is somewhere available copy-paste self contained C example functions or makros for diverse FFT algorithms FFT, QFT, Goertzel, radix-2, radix-4, split-radix, mixed-radix ... ? Which variables in FFT function must/should be defined as static, register ... for best performance ? What typical comes in to FFT function ? Pointer to already windowed array of samples ? What return FFT ? What exact physical dimensions return FFT , if input function dimensions was - voltage depends from (time) U=U(t) ? How from ...1024 or 2048 or 4096... FFT return values i calculate power magnitudes for all bands, and finally values for visual 10-20 hopping bars, like in Winamp , XMMS , QMMP ... ? If i exact know my FFT window size [for example 4096 samples] and window type , and it will be constant forever, is it possible to calculate window(,sine,cosine,) and store all values in constant array, so that in runtime it do not need be calculated , but i can just take values form array ? I have google_d about FFT and have found such proggie [see below] I have it little bit remixed , i generate pure sine frekwenz = 796.7285 HZ , and in output file i got so what : [ 71] 764.4287 Hz: Re= 0.0000000011142182 Im= 0.0000002368905824 M= 0.0000002368932028 [ 72] 775.1953 Hz: Re= 0.0000000011147694 Im= 0.0000003578625618 M= 0.0000003578642981 [ 73] 785.9619 Hz: Re= 0.0000000011164234 Im= 0.0000007207092628 M= 0.0000007207101275 [ 74] 796.7285 Hz: Re=-0.0000022785007614 Im=-0.5000000048748561 M= 0.5000000048800476 [ 75] 807.4951 Hz: Re= 0.0000000011098065 Im=-0.0000007304711756 M= 0.0000007304720186 [ 76] 818.2617 Hz: Re= 0.0000000011114605 Im=-0.0000003676273975 M= 0.0000003676290776 [ 77] 829.0283 Hz: Re= 0.0000000011120118 Im=-0.0000002466579503 M= 0.0000002466604569 ... ... ... [4019] 43270.9717 Hz: Re= 0.0000000011120118 Im= 0.0000002466579503 M= 0.0000002466604569 [4020] 43281.7383 Hz: Re= 0.0000000011114605 Im= 0.0000003676273975 M= 0.0000003676290776 [4021] 43292.5049 Hz: Re= 0.0000000011098065 Im= 0.0000007304711756 M= 0.0000007304720186 [4022] 43303.2715 Hz: Re=-0.0000022785015510 Im= 0.5000000048748419 M= 0.5000000048800334 [4023] 43314.0381 Hz: Re= 0.0000000011164234 Im=-0.0000007207092628 M= 0.0000007207101275 [4024] 43324.8047 Hz: Re= 0.0000000011147694 Im=-0.0000003578625618 M= 0.0000003578642981 [4025] 43335.5713 Hz: Re= 0.0000000011142182 Im=-0.0000002368905824 M= 0.0000002368932028 Where 43303.2715 + 796.7285 = 44100 or 44100 - 43303.2715 = 796.7285 Why Frequency 796.7285 is mirrored as Frequency 43303.2715 , and magnitude for both Frequencies is divided by 2 ???? Is here way direct calculate full magnitude and without Frequency mirroring , in band 0 Hz ... FSampl/2 ONLY , and not in full band - 0 Hz ... FSampl ?? In this case - how do i calculate corresponding frequency of each band that returns FFT, if sample frequency is 32K, 44K or 48K ? Pleaz do not point me to FFTW[3] and such libs , i must self write/combine code . Except FFTW has best short self contained 1D functions for copy-paste :) Each pointer and example is welcomed. Tnx in advance @ all. Alfs Kurmis. ==== #include <stdio.h> #include <stdlib.h> #include <string.h> #include <ctype.h> #include <math.h> #define BUFFER 4096 //2048 //#define BUFFER 256 # define M_PI_LD 3.1415926535897932384626433832795029L /* pi */ //void readwave(const char *filename, double *tr); short FFT(short int dir, long m, double *x, double *y); double fstep; int main (int argc, char *argv []) { int i; double f = 0.0; double amplitude = 0.; double samplerate = 44100., frekwenz=796.7285; double xarg=0. , argplus = 0. , pti ; double real[BUFFER], img[BUFFER]; // ==== ===== argplus = ( frekwenz *2.0*M_PI )/samplerate; //printf("Reading %d bytes from %s.rn", BUFFER, argv[1]); for(i=0;i<BUFFER;i++) { // xarg=0.0; argplus=0.2; floarArg =0.0001; pti = sin ( xarg ) /* * 32767.0*/ ; /*if ( kante>0 ) { if( pti > 0.0 ){ pti = 23767.0; }else{pti = -23767.0;} } /**/ xarg = xarg +argplus;//+floarArg; if ( xarg > (4.0*M_PI) ) xarg = xarg - (4.0*M_PI); real[i] = pti; } printf("F= %7.2f Hz nn", (samplerate*argplus)/(2*M_PI) ); memset(img, 0, sizeof(img)); /* Fill all the imaginary parts with zeros */ //fstep = (double) samplerate / (double) (BUFFER*2); fstep = (double) samplerate / (double) (BUFFER); printf("Frequency step : %10.6frn", fstep); FFT(1, /*11*/ 12, real, img); /* Fast Fourier Transform with 2^11 bins */ // FFT(1, 7, real, img); /* Fast Fourier Transform with 2^11 bins */ /* Write fourier transformed data to stdio */ i = 0; while(i < BUFFER) { amplitude = sqrt((real[i]*real[i]) + (img[i]*img[i])); //printf("(%4d) %.2f Hz: Re=%f Im=%f P=%frn", i, f, real[i], img[i], amplitude); printf("[%4d] %8.4f Hz: Re=%22.16f Im=%22.16f M=%22.16fn", i, f, real[i], img[i], amplitude); i++; f += fstep; } return 0 ; } // main /* This computes an in-place complex-to-complex FFT x and y are the real and imaginary arrays of 2^m points. dir = 1 gives forward transform dir = -1 gives reverse transform */ short FFT(short int dir, long m, double *x, double *y) { long n,i,i1,j,k,i2,l,l1,l2; double c1,c2,tx,ty,t1,t2,u1,u2,z; /* Calculate the number of points N = M ^2 */ n = 1; for (i=0;i<m;i++) n *= 2; printf("FFT -->> (n) 2 ^ %ld = %ldn", m, n); /* Do the bit reversal */ i2 = n >> 1; j = 0; for (i=0;i<n-1;i++) { if (i < j) { tx = x[i]; ty = y[i]; x[i] = x[j]; y[i] = y[j]; x[j] = tx; y[j] = ty; } k = i2; while (k <= j) { j -= k; k >>= 1; } j += k; } // for (i=0;i<n-1;i++) /* Compute the FFT */ c1 = -1.0; c2 = 0.0; l2 = 1; for (l=0;l<m;l++) { l1 = l2; l2 <<= 1; u1 = 1.0; u2 = 0.0; for (j=0;j<l1;j++) { for (i=j;i<n;i+=l2) { i1 = i + l1; t1 = u1 * x[i1] - u2 * y[i1]; t2 = u1 * y[i1] + u2 * x[i1]; x[i1] = x[i] - t1; y[i1] = y[i] - t2; x[i] += t1; y[i] += t2; } z = u1 * c1 - u2 * c2; u2 = u1 * c2 + u2 * c1; u1 = z; } c2 = sqrt((1.0 - c1) / 2.0); if (dir == 1) c2 = -c2; c1 = sqrt((1.0 + c1) / 2.0); } /* Scaling for forward transform */ if (dir == 1) { for (i=0;i<n;i++) { x[i] /= n; y[i] /= n; } } return(1); //return(TRUE); } ----