1 /* dft.c - discrete fourier transformation program */
4 * author: frank.zirkelbach@physik.uni-augsburg.de
12 int main(int argc,char **argv) {
16 fftw_complex *in,*out;
35 strcpy(src.file,argv[1]);
36 strcpy(dst.file,argv[2]);
39 fourier_init(&fourier,1);
46 bmp_cut_grab_bottom(&cut,&src,src.info.width,GRAB);
49 dst.height=cut.height;
52 size=cut.width*cut.height;
55 fourier.data_len[0]=cut.width;
56 fourier.data_len[1]=cut.height;
57 fourier_alloc_data(&fourier);
59 in=fftw_malloc(sizeof(fftw_complex)*size);
60 memset(in,0,size*sizeof(fftw_complex));
61 out=fftw_malloc(sizeof(fftw_complex)*size);
62 memset(out,0,size*sizeof(fftw_complex));
63 if(in==NULL||out==NULL) {
64 printf("fftw alloc error!\n");
67 plan=fftw_plan_dft_2d(cut.width,cut.height,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
70 mag=(double *)malloc(size*sizeof(double));
72 printf("unable to alloc mag memory\n");
79 in[x][0]=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
82 fourier.data[x].r=(src.map[x].r+src.map[x].g+src.map[x].b)/3;
86 // do the fourier trafo
90 fourier_dft_2d(&fourier);
94 printf("fourier done!\n");
96 // copy back the data, intensity = sqrt(r^2+i^2)
99 for(y=0;y<cut.height;y++) {
100 for(x=0;x<cut.width;x++) {
102 mag[offt]=sqrt(out[offt][0]*out[offt][0]+out[offt][1]*out[offt][1]);
104 mag[offt]=sqrt(fourier.ftdata[offt].r*fourier.ftdata[offt].r+fourier.ftdata[offt].i*fourier.ftdata[offt].i);
106 if(mag[offt]>max) max=mag[offt];
111 printf("found max: %f\n",max);
113 if(max!=0) scale=255.0/max;
117 printf("scaling image intensity: %f\n",scale);
119 for(y=0;y<dst.height;y++) {
120 for(x=0;x<dst.width;x++) {
121 // usual image processing order
122 offy=((y<dst.height/2)?y+dst.height/2:y-dst.height/2);
124 offy+=((x<dst.width/2)?x+dst.width/2:x-dst.width/2);
126 dst.map[offy].r=scale*mag[offt];
127 dst.map[offy].g=dst.map[offy].r;
128 dst.map[offy].b=dst.map[offy].r;
136 bmp_write_file(&dst);
139 fftw_destroy_plan(plan);
140 fftw_free(in); fftw_free(out);
142 fourier_shutdown(&fourier);
148 printf("done ...\n");