14 Mar
2011
14 Mar
'11
1:11 p.m.
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);
}
----