00001 #include "Transform.h"
00002 #include <stdio.h>
00003 #include <iostream>
00004 #include <fftw.h>
00005
00006 using namespace std;
00007 namespace SpecFVRNS
00008 {
00009
00010
00011 #define WISDOM "fvr_fftw.wis"
00012
00013 static fftwnd_plan g_plan_2D;
00014 static bool g_plan_exists = false;
00015 static dword g_x3d;
00016 static dword g_y3d;
00017 static dword g_z3d;
00018
00019 void initTransforms(void)
00020 {
00021 FILE* wisdom_file = fopen(WISDOM, "a+");
00022 if(!wisdom_file) {
00023 cout<<"Could not load wisdom."<<endl;
00024 return;
00025 }
00026 rewind(wisdom_file);
00027 fftw_import_wisdom_from_file(wisdom_file);
00028 fclose(wisdom_file);
00029 }
00030
00031 void destroyTransforms(void)
00032 {
00033 FILE* wisdom_file = fopen(WISDOM, "w");
00034 if(!wisdom_file) {
00035 cout << "Could not save wisdom."<<endl;
00036 return;
00037 }
00038 fftw_export_wisdom_to_file(wisdom_file);
00039 fclose(wisdom_file);
00040 }
00041
00042 void initTransform2D(dword XSize, dword YSize)
00043 {
00044
00045 g_plan_2D = fftw2d_create_plan(YSize, XSize, FFTW_BACKWARD, FFTW_MEASURE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00046 g_plan_exists=true;
00047 }
00048
00049 void transform2D(float* x)
00050 {
00051 if(g_plan_exists)
00052 fftwnd_one(g_plan_2D, (fftw_complex*) x, 0);
00053 else cout << "no plan" << endl;
00054 }
00055
00056 void destroyTransform2D(void)
00057 {
00058 if(g_plan_exists) {
00059 fftwnd_destroy_plan(g_plan_2D);
00060 g_plan_exists = false;
00061 }
00062 }
00063
00064 void initTransform3D(dword XSize, dword YSize, dword ZSize)
00065 {
00066 g_x3d = XSize;
00067 g_y3d = YSize;
00068 g_z3d = ZSize;
00069 }
00070
00071 void transform3D(float* x)
00072 {
00073 fftwnd_plan p;
00074 p = fftw3d_create_plan(g_z3d, g_y3d, g_x3d, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00075 fftwnd_one(p, (fftw_complex*) x, NULL);
00076 fftwnd_destroy_plan(p);
00077 dword size_i = g_z3d * g_y3d * g_x3d * 2;
00078 float size_f = (float) g_z3d * g_y3d * g_x3d;
00079 for (dword i = 0; i < size_i; ++i)
00080 x[i] /= size_f;
00081 }
00082
00083 void destroyTransform3D(void)
00084 {
00085 g_x3d = 0;
00086 g_y3d = 0;
00087 g_z3d = 0;
00088 }
00089
00090 void shift3D(float* x, dword XSize, dword YSize, dword ZSize)
00091 {
00092 dword xstep;
00093 dword ystep;
00094 dword i;
00095 dword j;
00096 dword k;
00097 dword index1;
00098 int start1 = 0;
00099 int val1 = -1;
00100 dword index2;
00101 int start2;
00102 int val2;
00103
00104 xstep = XSize * 2;
00105 ystep = YSize * XSize * 2;
00106
00107 index1 = 0;
00108 for(k = 0; k < ZSize; ++k)
00109 {
00110 index2 = index1;
00111 start2 = start1;
00112 val2 = val1;
00113 for(j = 0; j < YSize; ++j)
00114 {
00115 for(i = start2; i < xstep; i += 4)
00116 {
00117 x[index2+i] *= -1.0f;
00118 x[index2+i+1] *= -1.0f;
00119 }
00120 val2 *= -1;
00121 start2 += val2*2;
00122 index2 += xstep;
00123 }
00124 val1 *= -1;
00125 start1 += val1*2;
00126 index1 += ystep;
00127 }
00128 }
00129
00130 void shift2D(float* x, dword XSize, dword YSize)
00131 {
00132 dword xstep;
00133 dword i;
00134 dword j;
00135 dword index;
00136 int start = 0;
00137 int val = -1;
00138
00139 xstep = XSize * 2;
00140
00141 index = 0;
00142 for (j = 0; j < YSize; ++j) {
00143 for (i = start; i < xstep; i += 4) {
00144 x[index+i] *= -1.0f;
00145 x[index+i+1] *= -1.0f;
00146 }
00147 val *= -1;
00148 start += val * 2;
00149 index += xstep;
00150 }
00151 }
00152
00153 void shift_copy2D(float* dest, float* src, dword XSize, dword YSize)
00154 {
00155 dword i, j;
00156 dword xlast = XSize / 2;
00157 dword ylast = YSize / 2;
00158
00159 float* src_ptr = src;
00160 float* dest_ptr = dest;
00161
00162 for (j = 0; j < ylast; ++j) {
00163 for (i = 0; i < xlast; ++i) {
00164 *(dest_ptr++) = *(src_ptr++) * -1.0f;
00165 *(dest_ptr++) = *(src_ptr++) * -1.0f;
00166 *(dest_ptr++) = *(src_ptr++);
00167 *(dest_ptr++) = *(src_ptr++);
00168 }
00169 for (i = 0; i < xlast; ++i) {
00170 *(dest_ptr++) = *(src_ptr++);
00171 *(dest_ptr++) = *(src_ptr++);
00172 *(dest_ptr++) = *(src_ptr++) * -1.0f;
00173 *(dest_ptr++) = *(src_ptr++) * -1.0f;
00174 }
00175 }
00176 }
00177
00178 }