#include #include #include #include "utils.h" static void hartley(float array[], float work[], int num_el); /*....................................................................... Evaluate the Fast Hartley transform of a two dimensional array, overwriting the data array with the transform. The reader is referred to the comments at the start of the 1-D fast_hartley_transform() function for the advantages and organisation of the FHT. The data array, data_array[] should actually be a 1-D array arranged as a contiguous 2-D array in FORTRAN sort order. The number of elements along each axis are sent as xnum and ynum. Both must be a power of two ie: xnum=2^I and ynum=2^J where I and J are integers. */ int two_dim_FHT(float data_array[], int xnum, int ynum, int forward) { int i,j,y_element, inc_ad, inc_bc, half_xnum, half_ynum; float *work_array, *temp_array; float *ta,*tb,*tc,*td, temp; /* Signal an error if the number of elements on either axis of the array to be transformed is not a power of two. */ if(!is_pow_of_two(xnum) || !is_pow_of_two(ynum) ) { lprintf(stderr, "Illegal array size (%d,%d) - not a powers of two - sent to\nthe Fast Hartley Transform function.\n", xnum, ynum); return -1; }; /* Create a work array for the 1-D algorithm - equal to the length of the longest axis of the data array. */ if((work_array = (float *) calloc(((xnum>ynum)?xnum:ynum)+1, sizeof(float))) == NULL) { lprintf(stderr, "Insufficient memory available for work array of Hartley transform\n"); return -1; }; /* Since the second dimension of the data array is formed of non-contiguous elements, create another work array into which each y-array will be copied. */ if((temp_array= (float *) calloc(ynum+1, sizeof(float))) == NULL) { lprintf(stderr, "Insufficient memory available for work array of Hartley transform\n"); free(work_array); return -1; }; /* Transform the 2nd dimension of the 2D data array first - if there is one. */ if(ynum > 1) { for(i=0; i 0. The resulting transform is returned in data_array[] and is organised as follows: If the original data array had time increments of, dt, per channel: elements i=0,1,2..num_el/2 correspond to frequencies i/(dt*num_el), and subsequent elements correspond to frequencies (i-num_el)/(dt*num_el). For example if num_el=8 and dt=1, the frequencies of the channels would be: Element i: 0 1 2 3 4 5 6 7 Frequency: 0 1/8 2/8 3/8 4/8 -3/8 -2/8 -1/8 The real value corresponding to frequency f is the even part of the transform, H ie: [H(f)+H(-f)]/2, the imaginary value is minus the odd part: -[H(f)-H(-f)]/2 and its modulus is: [H(f)^2+H(-f)^2]/2. In terms of array elements these correspond to: real[i] = (data_array[i]+data_array[i-num_el])/2 imag[i] = (data_array[i-num_el]-data_array[i])/2 modulus[i] = (data_array[i]*data_array[i]+data_array[i-num_el]*data_array[i-num_el])/2 Except for i=0 which is wholly real: real[0]=data_array[0], imag[0]=0. */ int fast_hartley_transform(float data_array[], float work_array[], int num_el, int do_forward_transform) { int i; /* If the data array has only one element then its transform is the same as the data array so no operations are required. */ if(num_el == 1) return 0; /* Signal an error if the number of elements in the array to be transformed is not a power of two. */ if(!is_pow_of_two(num_el)) { lprintf(stderr, "Illegal array length (%d) - not a power of two - sent to\nthe Fast Hartley Transform function.\n", num_el); return -1; }; /* Perform the transform. */ hartley(data_array, work_array, num_el); /* Divide the transform by the number of elements in the array. This is only required on the forward transform. */ if(do_forward_transform) for(i=0;i 4) { /* Sort even and odd elements into the lower and upper halves of the work array. */ half_num = num_el/2; for(i=0,j=half_num; i < half_num; i++,j++) { work[i] = array[i+i]; work[j] = array[i+i+1]; }; /* Now transform these two new sub-arrays. */ hartley(&work[half_num], &array[half_num], num_el/2); hartley(work, array, num_el/2); /* Having transformed the two new arrays, we need to combine them into one array, thus forming the transform of the data array that was passed to this routine. If we denote the upper and lower sub-arrays as A(i) and B(i) respectively, the combination follows as: H(i) = A(i) + B(i)*cos(2.PI.j/num_el) + C(i)*sin(2.PI.j/num_el) With C(i) = B(i)*(num_el-i) (except that B(0) when i=0). and where j=0,1,2,...num_el-1 and i increments by one as j does, but restarts from zero when it reaches num_el/2. */ half_num = num_el/2; two_pi_div_n = 6.2831853072/num_el; for(i=0,j=0; j