/*---------------------------------------------------------------------------- 
 * 
 * NAME 
 * fft0  - Warm up and utilities for the FFT package. 
 * USAGE 
 * 
 * 
 * EXPLICIT INPUTS 
 * 
 * 
 * IMPLICIT INPUTS 
 * 
 * 
 * EXPLICIT OUTPUTS 
 * 
 * 
 * IMPLICIT OUTPUTS 
 * 
 * 
 * PROCEDURE 
 * 
 * 
 * HISTORY 
 * 18-Nov-82  John Schlag (jfs) at Carnegie-Mellon University 
 * Created. 
 * 25-Sep-88  Robert de Vries (rdv) at University of Technology Delft 
 * Added error handling. 
 */ 
 
#include <stdio.h> 
#include <math.h> 
#include "window.h" 
#include "display.h" 
#include "bblib.h" 
 
int 
grtocx (WINDOW *wp, WINDOW *wpc, int shift) 
{ 
    register short i, j; 
    register PIXEL     *pixel; 
    register short flip = 1; 
    COMPLEX        *data; 
    PIXEL        *pix_org; 
    int   xlen, ylen, res; 
    long  org_size; 
 
    pix_org = pixel = WindowData (wp); 
    org_size = WindowSize(wp); 
    xlen = WindowHeight(wp); 
    ylen = WindowWidth (wp); 
 
    if (wp == wpc)   /* allocate COMPLEX memory    */ 
 WindowData(wpc) = 0L; 
 
    res = set_impar( wpc, xlen, ylen, CPLX); 
    if (res != NO_ERROR) 
    { 
 if (wp == wpc) 
 { 
     WindowData(wp) = pix_org; 
     WindowSize(wp) = org_size; 
 } 
 return res; 
    } 
    data = (COMPLEX *) WindowData(wpc); 
 
    for (j = 0; j < xlen; j++) 
    { 
 if (shift) flip *= -1; 
 
 for (i = 0; i < ylen; i++) 
 { 
     if (shift) flip *= -1; 
     data->real = ((double) *pixel++) * flip; /* real value */ 
     data->imag = 0.; 
     data++; 
 } 
    } 
    if (wp == wpc ) free(pix_org); /* free the pixel memory  */ 
    return NO_ERROR; 
} 
 
int 
cxtogr (WINDOW *wpc, WINDOW *wpd) 
{ 
    int   i, j; 
    register PIXEL     *pixel; 
    register COMPLEX   *data; 
    double  scale_real = 0. , 
   scale_imag = 0.; 
    register double areal,aimag; 
    int   xlen, ylen, res; 
 
    data = (COMPLEX *) WindowData(wpc); 
    xlen = WindowWidth(wpc); 
    ylen = WindowHeight(wpc); 
 
    res = set_impar( wpd, xlen, ylen, GREY); 
    if (res != NO_ERROR) return res; 
    pixel = WindowData (wpd); 
 
    for (j = 0; j < ylen; j++) 
    { 
 for (i = 0; i < xlen; i++) 
 { 
     if (fabs(data->real) > scale_real) 
  scale_real = fabs(data->real); 
     if (fabs(data->imag) > scale_imag) 
  scale_imag = fabs(data->imag); 
     data++; 
 } 
    } 
 
    scale_real = 255.0 / sqrt(scale_real*scale_real + scale_imag*scale_imag); 
    data = (COMPLEX *) WindowData(wpc); 
 
    for (j = 0; j < ylen; j++) 
    { 
 for (i = 0; i < xlen; i++) 
 { 
     areal = data->real; 
     aimag = data->imag; 
     *pixel++ = (PIXEL)(sqrt(areal*areal + aimag*aimag) * scale_real); 
     data++; 
 } 
    } 
    return NO_ERROR; 
} 
 
void 
fft_filter (WINDOW *wp_in, WINDOW *wp_out) 
{ 
    register double areal, aimag; 
    register int i, j; 
    COMPLEX        *array_in, *array_out; 
    int   xlen, ylen; 
 
    xlen = WindowWidth(wp_in); 
    ylen = WindowHeight(wp_in); 
    array_in = (COMPLEX*) WindowData(wp_in); 
    array_out = (COMPLEX*) WindowData(wp_out); 
 
    for (i = 0; i < ylen; i++) 
    { 
 for (j = 0; j < xlen; j++) 
 { 
     areal = array_out->real; 
     aimag = array_out->imag; 
     array_out->real = array_in->real * areal - array_in->imag * aimag; 
     array_out->imag = array_in->imag * areal + array_in->real * aimag; 
     array_in++; 
     array_out++; 
 } 
    } 
} 
 
void 
fdis (WINDOW *wp_cplx, WINDOW *wp_mag, WINDOW *wp_real, WINDOW *wp_imag) 
{ 
    int   xlen, ylen; 
    short  i,j; 
    register PIXEL     *pixel_mag, 
         *pixel_real, 
         *pixel_imag; 
    register COMPLEX   *array; 
    double  scale_real = 0., 
   scale_imag = 0., 
   scale_mag; 
    register double areal, aimag; 
 
    array = (COMPLEX *) WindowData(wp_cplx); 
    xlen = WindowWidth(wp_cplx); 
    ylen = WindowHeight(wp_cplx); 
 
    for (i = 0; i < ylen; i++) 
    { 
 for (j = 0; j < xlen; j++) 
 { 
     if ((i != ylen/2) || (j != xlen/2)) 
     { 
  if (fabs(array->real) > scale_real) 
      scale_real = fabs(array->real); 
  if (fabs(array->imag) > scale_imag) 
      scale_imag = fabs(array->imag); 
  array++; 
     } 
 } 
    } 
 
    scale_mag  = 127.0 / log(scale_real * scale_real + scale_imag * scale_imag); 
    scale_real = 127.0 / log(scale_real); 
    scale_imag = 127.0 / log(scale_imag); 
 
    pixel_mag  = (PIXEL *) WindowData(wp_mag); 
    pixel_real = (PIXEL *) WindowData(wp_real); 
    pixel_imag = (PIXEL *) WindowData(wp_imag); 
    array      = (COMPLEX *) WindowData(wp_cplx); 
 
    printw (" mag %g real %g imag %g\n",scale_mag,scale_real,scale_imag); 
 
    for (i = 0; i < ylen; i++) 
    { 
 for (j = 0; j < xlen; j++) 
 { 
     areal = array->real; 
     aimag = array->imag; 
     array++; 
     *pixel_mag++ = log(areal*areal + aimag*aimag) * scale_mag + 127.0; 
     *pixel_real++ = log(fabs(areal)) * scale_real + 127.0; 
     *pixel_imag++ = log(fabs(aimag)) * scale_imag + 127.0; 
 } 
    } 
} 
