[LAD] FFT 1D real

Previous message: [thread] [date] [author]
Next message: [thread] [date] [author]
To: <linux-audio-dev@...>
Date: Monday, March 14, 2011 - 12:11 pm

---moq13001046607cb8cd571861071d4828323d9bae2c5a
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: quoted-printable

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=3DU(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=20
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 =3D
796.7285 HZ ,
and in output file i got so what :
[ 71] 764.4287 Hz: Re=3D 0.0000000011142182 Im=3D 0.0000002368905824 M=3D
0.0000002368932028
[ 72] 775.1953 Hz: Re=3D 0.0000000011147694 Im=3D 0.0000003578625618 M=3D
0.0000003578642981
[ 73] 785.9619 Hz: Re=3D 0.0000000011164234 Im=3D 0.0000007207092628 M=3D
0.0000007207101275
[ 74] 796.7285 Hz: Re=3D-0.0000022785007614 Im=3D-0.5000000048748561 M=3D
0.5000000048800476
[ 75] 807.4951 Hz: Re=3D 0.0000000011098065 Im=3D-0.0000007304711756 M=3D
0.0000007304720186
[ 76] 818.2617 Hz: Re=3D 0.0000000011114605 Im=3D-0.0000003676273975 M=3D
0.0000003676290776
[ 77] 829.0283 Hz: Re=3D 0.0000000011120118 Im=3D-0.0000002466579503 M=3D
0.0000002466604569
...
...
...
[4019] 43270.9717 Hz: Re=3D 0.0000000011120118 Im=3D 0.0000002466579503
M=3D 0.0000002466604569
[4020] 43281.7383 Hz: Re=3D 0.0000000011114605 Im=3D 0.0000003676273975
M=3D 0.0000003676290776
[4021] 43292.5049 Hz: Re=3D 0.0000000011098065 Im=3D 0.0000007304711756
M=3D 0.0000007304720186
[4022] 43303.2715 Hz: Re=3D-0.0000022785015510 Im=3D 0.5000000048748419
M=3D 0.5000000048800334
[4023] 43314.0381 Hz: Re=3D 0.0000000011164234 Im=3D-0.0000007207092628
M=3D 0.0000007207101275
[4024] 43324.8047 Hz: Re=3D 0.0000000011147694 Im=3D-0.0000003578625618
M=3D 0.0000003578642981
[4025] 43335.5713 Hz: Re=3D 0.0000000011142182 Im=3D-0.0000002368905824
M=3D 0.0000002368932028
Where
43303.2715 + 796.7285 =3D 44100 or 44100 - 43303.2715 =3D 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.
=3D=3D=3D=3D
#include
#include
#include
#include
#include
#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 =3D 0.0;
double amplitude =3D 0.; double samplerate =3D 44100.,
frekwenz=3D796.7285; double xarg=3D0. , argplus =3D 0. , pti ;
double real[BUFFER], img[BUFFER];
// =3D=3D=3D=3D =3D=3D=3D=3D=3D
argplus =3D ( frekwenz *2.0*M_PI )/samplerate;
//printf("Reading %d bytes from %s.rn", BUFFER, argv[1]);
for(i=3D0;i0 ) { if( pti > 0.0 ){ pti =3D 23767.0;
}else{pti =3D -23767.0;} } /**/
xarg =3D xarg +argplus;//+floarArg;
if ( xarg > (4.0*M_PI) ) xarg =3D xarg - (4.0*M_PI);
real[i] =3D pti;
}
printf("F=3D %7.2f Hz nn", (samplerate*argplus)/(2*M_PI) );
memset(img, 0, sizeof(img)); /* Fill all the imaginary parts with
zeros */
//fstep =3D (double) samplerate / (double) (BUFFER*2);
fstep =3D (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 =3D 0;
while(i < BUFFER)
{
amplitude =3D sqrt((real[i]*real[i]) + (img[i]*img[i]));
//printf("(%4d) %.2f Hz: Re=3D%f Im=3D%f P=3D%frn", i, f, real[i],
img[i], amplitude);
printf("[%4d] %8.4f Hz: Re=3D%22.16f Im=3D%22.16f M=3D%22.16fn", i,
f, real[i], img[i], amplitude);
i++;
f +=3D 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 =3D 1 gives forward transform
dir =3D -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 =3D M ^2 */
n =3D 1; for (i=3D0;i> (n) 2 ^ %ld =3D %ldn", m, n);
/* Do the bit reversal */
i2 =3D n >> 1;
j =3D 0;
for (i=3D0;i>=3D 1;
}
j +=3D k;
} // for (i=3D0;iI am Physics student, and i wanna write referat abou=
t Fourier transformations, also about FFT 1D real case.Hope this is b=
est place for ask this, and here are best experts.I wanaa demon=
strate how diverse window functions changes measured spectrum,how muc=
h CPU ressources take diverse FFT algorithms ...Yet i understan=
d [i hope so] how works windowing .Is somewhere available copy-=
paste self contained C example functions or makros for diverse FFT algorith=
msFFT, QFT, Goertzel, radix-2, radix-4, split-radix, mixed-radix ... =
?Which variables in FFT function must/should be defined as stat=
ic, register ... for best performance ?What typical comes in to FFT f=
unction ? Pointer to already windowed array of samples ?What return F=
FT ?What exact physical dimensions return FFT , if input function dim=
ensions was - voltage depends from (time) U=3DU(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 valu=
es 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 =3D 796.7285 HZ , and in out=
put file i got so what :[ 71] 764.4287 Hz: Re=3D 0.000000001114=
2182 Im=3D 0.0000002368905824 M=3D 0.0000002368932028[ 72] 775.1953 H=
z: Re=3D 0.0000000011147694 Im=3D 0.0000003578625618 M=3D 0.000000357864298=
1[ 73] 785.9619 Hz: Re=3D 0.0000000011164234 Im=3D 0.0000007207092628=
M=3D 0.0000007207101275[ 74] 796.7285 Hz: Re=3D-0.0000022785007614 I=
m=3D-0.5000000048748561 M=3D 0.5000000048800476[ 75] 807.4951 Hz: Re=
=3D 0.0000000011098065 Im=3D-0.0000007304711756 M=3D 0.0000007304720186[ 76] 818.2617 Hz: Re=3D 0.0000000011114605 Im=3D-0.0000003676273975 M=3D=
0.0000003676290776[ 77] 829.0283 Hz: Re=3D 0.0000000011120118 Im=3D-=
0.0000002466579503 M=3D 0.0000002466604569.........=
[4019] 43270.9717 Hz: Re=3D 0.0000000011120118 Im=3D 0.0000002466579503 M=
=3D 0.0000002466604569[4020] 43281.7383 Hz: Re=3D 0.0000000011114605 =
Im=3D 0.0000003676273975 M=3D 0.0000003676290776[4021] 43292.5049 Hz:=
Re=3D 0.0000000011098065 Im=3D 0.0000007304711756 M=3D 0.0000007304720186<=
br />[4022] 43303.2715 Hz: Re=3D-0.0000022785015510 Im=3D 0.500000004874841=
9 M=3D 0.5000000048800334[4023] 43314.0381 Hz: Re=3D 0.00000000111642=
34 Im=3D-0.0000007207092628 M=3D 0.0000007207101275[4024] 43324.8047 =
Hz: Re=3D 0.0000000011147694 Im=3D-0.0000003578625618 M=3D 0.00000035786429=
81[4025] 43335.5713 Hz: Re=3D 0.0000000011142182 Im=3D-0.000000236890=
5824 M=3D 0.0000002368932028Where43303.2715 + 796.7285 =
=3D 44100 or 44100 - 43303.2715 =3D 796.7285Why Frequency 796.7285 is=
mirrored as Frequency 43303.2715 , and magnitude for both Frequencies is d=
ivided by 2 ????Is here way direct calculate full magnitude and witho=
ut Frequency mirroring , in band 0 Hz ... FSampl/2 ONLY ,and not in f=
ull band - 0 Hz ... FSampl ??In this case - how do i calculate corres=
ponding 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 con=
tained 1D functions for copy-paste :)Each pointer and example i=
s welcomed.Tnx in advance @ all.Alfs Kurmis.<=
br />=3D=3D=3D=3D#include <std=
io.h>#include <stdlib.h>#include <string.h>#include <ctype.h>#include <math.h>#define B=
UFFER 4096 //2048//#define BUFFER 256# define M_PI_LD &#16=
0;      3.1415926535897932384626433832795029L  /* =
pi *///void readwave(const char *filename, double *tr);short FF=
T(short int dir, long m, double *x, double *y);double fstep;int main (int argc, char *argv []){  int i;=
  double f =3D 0.0;  double amplitude =3D 0.;   d=
ouble  samplerate =3D 44100., frekwenz=3D796.7285; double xarg=3D0. , =
argplus =3D 0. , pti ;  double real[BUFFER], img[BUFFER];/=
/  =3D=3D=3D=3D  =3D=3D=3D=3D=3D  argplus =3D ( =
frekwenz *2.0*M_PI )/samplerate;  //printf("Reading %d bytes fro=
m %s.\r\n", BUFFER, argv[1]);   for(i=3D0;i<BUFFER;i++)<=
br />    {  //  xarg=3D0.0;  argplus=3D0.2;&#=
160;  floarArg =3D0.0001;       pt=
i =3D sin ( xarg ) /* * 32767.0*/ ;     &#16=
0; /*if ( kante>0 ) { if( pti > 0.0 ){ pti =3D 23767.0; }else{pti =3D=
-23767.0;}  }  /**/       xa=
rg =3D xarg +argplus;//+floarArg;      =
   if ( xarg > (4.0*M_PI) ) xarg =3D xarg - (4.0*M_PI);&=
#160;      real[i] =3D pti;   &#16=
0; }  printf("F=3D  %7.2f Hz \n\n", (samplerate*argplus)/(2=
*M_PI) );  memset(img, 0, sizeof(img)); /* Fill all the im=
aginary parts with zeros */   //fstep =3D (double) samplera=
te / (double) (BUFFER*2);   fstep =3D (double) samplerate /=
(double) (BUFFER);   printf("Frequency step : %10.6f\r\n",=
fstep);   FFT(1, /*11*/ 12, real, img); /* Fast Four=
ier Transform with 2^11 bins */   // FFT(1,  7, real, =
img); /* Fast Fourier Transform with 2^11 bins */  /* Writ=
e fourier transformed data to stdio */  i =3D 0;  whi=
le(i < BUFFER)   {    amplitude =3D sqrt=
((real[i]*real[i]) + (img[i]*img[i]));    //printf("(%=
4d) %.2f Hz: Re=3D%f Im=3D%f P=3D%f\r\n", i, f, real[i], img[i], amplitude)=
;    printf("[%4d] %8.4f Hz: Re=3D%22.16f  Im=3D%=
22.16f  M=3D%22.16f\n", i, f, real[i], img[i], amplitude); =
   i++;    f +=3D 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 =3D&#160=
; 1 gives forward transform   dir =3D -1 gives reverse tran=
sform */short FFT(short int dir, long m, double *x, double *y)<=
br />{   long n,i,i1,j,k,i2,l,l1,l2;   doub=
le c1,c2,tx,ty,t1,t2,u1,u2,z;   /* Calculate the numb=
er of points  N =3D M ^2 */   n =3D 1;  &#16=
0; for (i=3D0;i<m;i++)  n *=3D 2;    &#16=
0;             =
printf("FFT -->> (n) 2 ^ %ld =3D  %ld\n", m, n); &#16=
0; /* Do the bit reversal */   i2 =3D n >> 1;&#=
160;  j =3D 0;   for (i=3D0;i<n-1;i++) {&#160=
;     if (i < j) {    &#16=
0;    tx =3D x[i];       =
  ty =3D y[i];         x=
[i] =3D x[j];       y[i] =3D y[j];&#160=
;        x[j] =3D tx;   &=
#160;     y[j] =3D ty;    &#1=
60; }      k =3D i2;   =
   while (k <=3D j) {     &#160=
;   j -=3D k;       &#16=
0; k >>=3D 1;      }  =
    j +=3D k;   } // for (i=3D0;i<n-1;i++=
)   /* Compute the FFT */   c1 =3D -1=
.0;    c2 =3D 0.0;   l2 =3D 1; &=
#160; for (l=3D0;l<m;l++) {      l1 =3D l=
2;      l2 <<=3D 1;  &=
#160;   u1 =3D 1.0;       u2 =3D 0=
.0;      for (j=3D0;j<l1;j++) {&#16=
0;        for (i=3Dj;i<n;i+=3Dl2) {            i1 =
=3D i + l1;         &#16=
0;  t1 =3D u1 * x[i1] - u2 * y[i1];  t2 =3D u1 * y[i1] + u2 * x[i=
1];           =
x[i1] =3D x[i] - t1;         =
    y[i1] =3D y[i] - t2;     =
       x[i] +=3D t1;    &=
#160;           &#16=
0;   y[i] +=3D t2;      &#160=
;  }         z =3D =
u1 * c1 - u2 * c2;         u=
2 =3D u1 * c2 + u2 * c1;       &#1=
60; u1 =3D z;      }   =
   c2 =3D sqrt((1.0 - c1) / 2.0);    &#=
160; if (dir =3D=3D 1)        &#16=
0; c2 =3D -c2;      c1 =3D sqrt((1.0 + c1) /=
2.0);   }   /* Scaling for forward t=
ransform */   if (dir =3D=3D 1) {   &#=
160;  for (i=3D0;i<n;i++) {     &#16=
0;   x[i] /=3D n;    y[i] /=3D n; &#160=
;    }   }    &#1=
60; return(1); //return(TRUE);} ----  
=0A
---moq13001046607cb8cd571861071d4828323d9bae2c5a--

Previous message: [thread] [date] [author]
Next message: [thread] [date] [author]

Messages in current thread:
[LAD] FFT 1D real, Alfs Kurmis, (Mon Mar 14, 12:11 pm)
Re: [LAD] FFT 1D real, Fons Adriaensen, (Mon Mar 14, 4:02 pm)