home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libfft / fft2 / fast.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  9.0 KB  |  429 lines
  1. /* *****************************************************************************
  2. *
  3. * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  4. * All Rights Reserved.
  5. *
  6. * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  7. * the contents of this file may not be disclosed to third parties, copied or
  8. * duplicated in any form, in whole or in part, without the prior written
  9. * permission of Silicon Graphics, Inc.
  10. *
  11. * RESTRICTED RIGHTS LEGEND:
  12. * Use, duplication or disclosure by the Government is subject to restrictions
  13. * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  14. * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  15. * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  16. * rights reserved under the Copyright Laws of the United States.
  17. *
  18. ***************************************************************************** */
  19. #include <gl.h>
  20. #include <device.h>
  21. #include <math.h>
  22. #include "fast.h"
  23. #include "event.h"
  24. #include "fft.h"
  25.  
  26. #define MIN(a,b) ( ((a) < (b)) ? (a) : (b))
  27.  
  28. void display();
  29. void InitFast();
  30. void matgen();
  31. void dir_fft(float *), inv_fft(float *);
  32. void inifil();
  33. void propagate();
  34.  
  35. double second();
  36.  
  37.  
  38. extern int main_menu, forward_menu, backward_menu, quit_menu;
  39.  
  40. int do_exit;
  41. int ask_quit;
  42.  
  43. int do_propagate;
  44.  
  45. int x_size, y_size;
  46. int this_time, this_deltat;
  47. int this_proc;
  48.  
  49. float *mat, *save, *filter, *power, *copy;
  50. int * mycolors;
  51. int *ipvt, lda, info, x_0, y_0;
  52.  
  53. double this_t, this_flop;
  54.  
  55. char title[128];
  56.  
  57.  
  58. fast (progname)
  59. char *progname;
  60. {
  61.     InitFast();
  62.     Init(progname);
  63.  
  64.     do_exit = FALSE;
  65.  
  66.     while(! do_exit) {
  67.         event();
  68.         if( do_propagate)
  69.             while( qtest() == 0)
  70.             propagate();
  71.     }
  72. }
  73.  
  74. do_frame(){
  75.     
  76.     char str[32];
  77.     float v[2];
  78.  
  79.    viewport(0,WIN_SIZE+WIN_RIGHT,0,WIN_SIZE);
  80.  
  81.     cpack (oBLACK);
  82.     clear ();
  83.  
  84.      cpack (oGREEN);
  85.     ortho2(0,WIN_SIZE+WIN_RIGHT,0, WIN_SIZE);
  86.  
  87.     sprintf( str, "SIZE = %d x %d", x_size, y_size);
  88.     cmov2i( 250, WIN_SIZE - 30);
  89.     charstr(str);
  90.  
  91.     cmov2i( 100, 20);
  92.     charstr(title);
  93.  
  94.     bgnpolygon();
  95.     v[0] = WIN_SIZE;
  96.     v[1] = WIN_SIZE/2;
  97.     cpack(oBLACK);
  98.     v2f(v);
  99.     v[1] -= 100;
  100.     cpack(0x007fff);
  101.     v2f(v);
  102.     v[0] += 50;
  103.     v2f(v);
  104.     v[1] += 100;
  105.     cpack(oBLACK);
  106.     v2f(v);
  107.     endpolygon();
  108.     v[1] -= 100;
  109.     bgnpolygon();
  110.     cpack(0x007fff);
  111.     v2f(v);
  112.     v[1] -= 100;
  113.     cpack(0x00ffff);
  114.     v2f(v);
  115.     v[0] -= 50;
  116.     v2f(v);
  117.     v[1] += 100;
  118.     cpack(0x007fff);
  119.     v2f(v);
  120.     endpolygon();
  121.  
  122.     bgnpolygon();
  123.     v[0] = WIN_SIZE;
  124.     v[1] = WIN_SIZE/2;
  125.     cpack(oBLACK);
  126.     v2f(v);
  127.     v[1] += 100;
  128.     cpack(0xff7f00);
  129.     v2f(v);
  130.     v[0] += 50;
  131.     v2f(v);
  132.     cpack(0xff7f00);
  133.     v2f(v);
  134.     v[1] -= 100;
  135.     cpack(oBLACK);
  136.     v2f(v);
  137.     endpolygon();
  138.     v[1] += 100;
  139.     bgnpolygon();
  140.     cpack(0xff7f00);
  141.     v2f(v);
  142.     v[1] += 100;
  143.     cpack(oCYAN);
  144.     v2f(v);
  145.     v[0] -= 50;
  146.     v2f(v);
  147.     v[1] -= 100;
  148.     cpack(0xff7f00);
  149.     v2f(v);
  150.     endpolygon();
  151.     cpack(oWHITE);
  152.     rect(WIN_SIZE,WIN_SIZE/2-200,WIN_SIZE+50,WIN_SIZE/2+200);
  153.  
  154.     cmov2i( WIN_SIZE+55, WIN_SIZE/2 - 200 -4);
  155.     charstr(" -Max");
  156.  
  157.     cmov2i( WIN_SIZE+55, WIN_SIZE/2 + 200 -4);
  158.     charstr(" +Max");
  159.  
  160.     cmov2i( WIN_SIZE+55, WIN_SIZE/2 -4);
  161.     charstr(" 0.0 ");
  162.  
  163. }
  164. ShowFast()
  165. {
  166.     x_0 = (WIN_SIZE/2)-(x_size/2);
  167.     y_0 = (WIN_SIZE/2)-(y_size/2);
  168.     do_frame();
  169.     viewport( y_0, y_0 + y_size + 2, x_0, x_0 + x_size + 2);
  170.     ortho2(-1, y_size+1, -1, x_size+1);
  171.     DrawFast();
  172. }
  173.  
  174. void do_quit()
  175. {
  176.     if(mat)
  177.         free(mat);
  178.     gexit();
  179.     exit(0);
  180. }
  181.  
  182.  
  183. void do_rightmouse()
  184. {
  185.     ask_quit = 0;
  186.  
  187.     dopup(main_menu);
  188.  
  189.     if(ask_quit)
  190.     dopup(quit_menu);
  191.  
  192. }
  193.  
  194. void display()
  195. {
  196.     ShowFast();
  197.     swapbuffers();
  198. }
  199.  
  200. void do_resize(){
  201.     
  202.     lda = y_size + 1;
  203.  
  204.     if( mat){ 
  205.     mat  = (float *)realloc( mat, lda * x_size * sizeof(float));
  206.     copy  = (float *)realloc( copy, lda * x_size * sizeof(float));
  207.     power  = (float *)realloc( power, lda * x_size * sizeof(float));
  208.     filter  = (float *)realloc( filter, lda * x_size * sizeof(float));
  209.     mycolors  = (int *)realloc( mycolors, lda * x_size * sizeof(int));
  210.     save  = (float *)realloc( save, (y_size + 3*x_size + 4 * FACTOR_SPACE) * sizeof(float));
  211.     }
  212.     else{ 
  213.     mat   = (float *)malloc( lda * x_size * sizeof(float));
  214.     copy   = (float *)malloc( lda * x_size * sizeof(float));
  215.     power   = (float *)malloc( lda * x_size * sizeof(float));
  216.     filter   = (float *)malloc( lda * x_size * sizeof(float));
  217.     mycolors   = (int *)malloc( lda * x_size * sizeof(int));
  218.     save  = (float *)malloc( (y_size + 3*x_size + 4 * FACTOR_SPACE) * sizeof(float));
  219.     }
  220.     this_deltat = 1;
  221.     this_time = 0;
  222.  
  223.     matgen();
  224.     inifil();
  225.     main_menu = forward_menu;
  226.     sprintf(title,"Time T = %d", this_time);
  227.     display();
  228. }
  229. void add_noise() {
  230.  
  231.     matnoise_( mat, &lda, &y_size, &x_size);
  232.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  233.     sprintf(title,(main_menu == forward_menu) ? "XY domain" : "FOURIER domain");
  234.     display();
  235. }
  236.  
  237.  
  238. void size256_256(){
  239.  
  240.     x_size = 256;
  241.     y_size = 256;
  242.     do_resize();
  243. }
  244.  
  245. void size288_300(){
  246.  
  247.     x_size = 288;
  248.     y_size = 300;
  249.     do_resize();
  250. }
  251.  
  252. void size512_512(){
  253.  
  254.     x_size = 512;
  255.     y_size = 512;
  256.     do_resize();
  257. }
  258.  
  259. void dir_fft(float *mat){
  260.     float alpha = 1./(float)(x_size*y_size);
  261.  
  262.     this_t = second();
  263.     sfft2d( -1, y_size, x_size, mat, lda, save);
  264.     this_t = second() - this_t;
  265. }
  266. void inv_fft(float *mat){
  267.  
  268.     this_t = second();
  269.     sfft2d( 1, y_size, x_size, mat, lda, save);
  270.     this_t = second() - this_t;
  271. }
  272.  
  273. void this_quit(){
  274.  
  275.     ask_quit = 1;
  276. }
  277. void InitFast(){
  278.  
  279.     this_proc = mp_numthreads_();
  280.     x_size = 512;
  281.     y_size = 512;
  282. }
  283. #define SCALE        20
  284. #define XLENGTH        19*SCALE
  285. #define YLENGTH        11*SCALE
  286.  
  287. void myrec( float *pt, int xl, int yl, double val){
  288.   int i, j;
  289.     for( j = 0 ; j < xl; j++, pt += lda)
  290.     for( i = 0; i < yl ; i++)
  291.         pt[i] += val;
  292. }
  293. myF( float *pt, double val){
  294.     myrec( pt, SCALE, 7*SCALE, val);
  295.     myrec( pt+SCALE*lda+3*SCALE, 2*SCALE, SCALE, val);
  296.     myrec( pt+SCALE*lda+6*SCALE, 3*SCALE, SCALE, val);
  297. }
  298. myT( float *pt, double val){
  299.     myrec( pt+0*SCALE*lda+6*SCALE, 5*SCALE, 1*SCALE, val);
  300.     myrec( pt+2*SCALE*lda+0*SCALE, 1*SCALE, 6*SCALE, val);
  301. }
  302. void mywrite( float *pt){
  303.     myrec( pt, XLENGTH, YLENGTH, 1.);
  304.     myF( pt + 2 * SCALE +  2 * SCALE*lda, -2.);
  305.     myF( pt + 2 * SCALE +  7 * SCALE*lda, -2.);
  306.     myT( pt + 2 * SCALE + 12 * SCALE*lda, -2.);
  307. }
  308.  
  309. void add_stick(){
  310.     int xor, yor;
  311.     xor = x_size/2 - XLENGTH/2;
  312.     yor = y_size/2;
  313.     myrec(mat+yor+xor*lda, XLENGTH, SCALE, 1.);
  314. }
  315. void add_box(){
  316.     int xor, yor;
  317.     xor = x_size/2 - XLENGTH/2;
  318.     yor = y_size/2 - YLENGTH/2;
  319.     mywrite( mat+yor+xor*lda);
  320. }
  321. void add_circle(){
  322.     int i,j;
  323.     int xor, yor;
  324.     float radius, width, dx, dy, z;
  325.  
  326.     radius = MIN( x_size, y_size) / 7;    
  327.     width = radius/4;
  328.     width = width * width;
  329.  
  330.     xor = x_size/2;
  331.     yor = y_size/2;
  332.  
  333.     for(i = 0; i < x_size; i++){ 
  334.     dx = i - xor;
  335.     for(j = 0; j < y_size; j++){
  336.         dy = j - yor;
  337.         z = dx * dx + dy * dy;
  338.         z = sqrt(z) - radius;
  339.         z = z * z;
  340.         mat[i*lda+j] = width / (z + width);
  341.     }
  342.     }        
  343. }
  344. #define X_OFFSET 7
  345. #define Y_OFFSET 13
  346.  
  347. void inifil(){
  348.     float alpha = 1./(float)(x_size*y_size);
  349.     bzero( filter, lda * x_size * sizeof(float));
  350.     filter[X_OFFSET*lda + Y_OFFSET] = 1.0;
  351.     sfft2d( -1, y_size, x_size, filter, lda, save);
  352. /*
  353.     inifil_(filter, &lda, &x_size, &y_size);
  354.     sscal2d_( &y_size, &x_size, &alpha, mat, &lda);
  355. */
  356. }
  357. void matgen(){ 
  358.     float alpha = 1./(float)(x_size*y_size);
  359.     bzero( mat,  lda * x_size * sizeof(float));
  360. /*
  361.     add_box();
  362.     add_stick();
  363.     inv_fft(mat);
  364.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  365.     display();
  366. */
  367.     add_circle_( &y_size, &x_size, mat, &lda);
  368.  
  369.     sfft2di( y_size, x_size, save);
  370.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  371.     display();
  372.     sscal2d_( &y_size, &x_size, &alpha, mat, &lda);
  373.     dir_fft(mat);
  374.  
  375.     do_propagate = 0;
  376. }
  377.  
  378. void propagate(){
  379.     sprod2d( y_size, x_size, mat, lda, filter, lda);
  380.     bcopy( mat, copy, x_size*lda*sizeof(float));
  381.     inv_fft(copy);
  382.     matcolor_(&y_size, &x_size, copy, mycolors, &lda);
  383.     this_time += this_deltat;
  384.     sprintf(title,"Time T = %d, 2D FFT done in %.2f Sec.", this_time, this_t);
  385.     display();
  386.     do_propagate ++;
  387. }
  388.  
  389. void faster(){
  390.     if( this_time % (2*this_deltat))
  391.     propagate();
  392.     sprod2d( y_size, x_size, filter, lda, filter, lda);
  393.     this_deltat *= 2.;
  394. }
  395.  
  396. void reverse(){
  397.     reverse_(filter, &lda, &y_size, &x_size);
  398.     this_deltat *= -1;
  399.     propagate();
  400. }
  401.  
  402. void do_power(){
  403.     matpower_(&y_size, &x_size, mat, power, &lda);
  404.     matlog_(&y_size, &x_size, power, &lda);
  405.     matcolor_(&y_size, &x_size, power, mycolors, &lda);
  406.     sprintf(title,"Power Spectrum (Log Scale)");
  407.     display();
  408. }
  409.  
  410. void do_fourier(){
  411.     matpower_(&y_size, &x_size, mat, power, &lda);
  412.     matcolor_(&y_size, &x_size, power, mycolors, &lda);
  413.     sprintf(title,"Power Spectrum");
  414.     display();
  415. }
  416. void do_filter(){
  417.     bcopy( filter, mat, lda*x_size*sizeof(float));
  418.     matcolor_(&y_size, &x_size, mat, mycolors, &lda);
  419.     sprintf(title,"Power Spectrum");
  420.     display();
  421. }
  422. void do_filpow(){
  423.     bcopy( filter, power, lda*x_size*sizeof(float));
  424.     matpower_(&y_size, &x_size, power, power, &lda);
  425.     matcolor_(&y_size, &x_size, power, mycolors, &lda);
  426.     sprintf(title,"Power Spectrum");
  427.     display();
  428. }
  429.