00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00028 #include <stdio.h>
00029 #include <math.h>
00030 #include <string.h>
00031 #include <stdlib.h>
00032 #include <complex.h>
00033
00034 #include "nfft3util.h"
00035 #include "nfft3.h"
00036
00037 void ndft_horner_trafo(nfft_plan *ths)
00038 {
00039 int j,k;
00040 double _Complex *f_hat_k, *f_j;
00041 double _Complex exp_omega_0;
00042
00043 for(j=0, f_j=ths->f; j<ths->M_total; j++, f_j++)
00044 (*f_j) =0;
00045
00046 for(j=0, f_j=ths->f; j<ths->M_total; j++, f_j++)
00047 {
00048 exp_omega_0 = cexp(+2*PI*_Complex_I*ths->x[j]);
00049 for(k=0, f_hat_k= ths->f_hat; k<ths->N[0]; k++, f_hat_k++)
00050 {
00051 (*f_j)+=(*f_hat_k);
00052 (*f_j)*=exp_omega_0;
00053 }
00054 (*f_j)*=cexp(-PI*_Complex_I*ths->N[0]*ths->x[j]);
00055 }
00056 }
00057
00058 void ndft_pre_full_trafo(nfft_plan *ths, double _Complex *A)
00059 {
00060 int j,k;
00061 double _Complex *f_hat_k, *f_j;
00062 double _Complex *A_local;
00063
00064 for(j=0, f_j=ths->f; j<ths->M_total; j++, f_j++)
00065 (*f_j) =0;
00066
00067 for(j=0, f_j=ths->f, A_local=A; j<ths->M_total; j++, f_j++)
00068 for(k=0, f_hat_k= ths->f_hat; k<ths->N[0]; k++, f_hat_k++, A_local++)
00069 (*f_j) += (*f_hat_k)*(*A_local);
00070 }
00071
00072 void ndft_pre_full_init(nfft_plan *ths, double _Complex *A)
00073 {
00074 int j,k;
00075 double _Complex *f_hat_k, *f_j, *A_local;
00076
00077 for(j=0, f_j=ths->f, A_local=A; j<ths->M_total; j++, f_j++)
00078 for(k=0, f_hat_k= ths->f_hat; k<ths->N[0]; k++, f_hat_k++, A_local++)
00079 (*A_local) = cexp(-2*PI*_Complex_I*(k-ths->N[0]/2)*ths->x[j]);
00080
00081 }
00082
00083 void ndft_time(int N, int M, unsigned test_ndft, unsigned test_pre_full)
00084 {
00085 int r;
00086 double t, t_ndft, t_horner, t_pre_full, t_nfft;
00087 double _Complex *A;
00088
00089 nfft_plan np;
00090
00091 printf("%d\t%d\t",N, M);
00092
00093 nfft_init_1d(&np, N, M);
00094
00096 nfft_vrand_shifted_unit_double(np.x, np.M_total);
00097
00098 if(test_pre_full)
00099 {
00100 A=(double _Complex*)nfft_malloc(N*M*sizeof(double _Complex));
00101 ndft_pre_full_init(&np, A);
00102 }
00103
00105 nfft_vrand_unit_complex(np.f_hat, np.N_total);
00106
00108 if(test_ndft)
00109 {
00110 t_ndft=0;
00111 r=0;
00112 while(t_ndft<0.1)
00113 {
00114 r++;
00115 t=nfft_second();
00116 ndft_trafo(&np);
00117 t=nfft_second()-t;
00118 t_ndft+=t;
00119 }
00120 t_ndft/=r;
00121
00122 printf("%.2e\t",t_ndft);
00123 }
00124 else
00125 printf("nan\t\t");
00126
00128 t_horner=0;
00129 r=0;
00130 while(t_horner<0.1)
00131 {
00132 r++;
00133 t=nfft_second();
00134 ndft_horner_trafo(&np);
00135 t=nfft_second()-t;
00136 t_horner+=t;
00137 }
00138 t_horner/=r;
00139
00140 printf("%.2e\t", t_horner);
00141
00143 if(test_pre_full)
00144 {
00145 t_pre_full=0;
00146 r=0;
00147 while(t_pre_full<0.1)
00148 {
00149 r++;
00150 t=nfft_second();
00151 ndft_pre_full_trafo(&np,A);
00152 t=nfft_second()-t;
00153 t_pre_full+=t;
00154 }
00155 t_pre_full/=r;
00156
00157 printf("%.2e\t", t_pre_full);
00158 }
00159 else
00160 printf("nan\t\t");
00161
00162 t_nfft=0;
00163 r=0;
00164 while(t_nfft<0.1)
00165 {
00166 r++;
00167 t=nfft_second();
00168 nfft_trafo(&np);
00169 t=nfft_second()-t;
00170 t_nfft+=t;
00171 }
00172 t_nfft/=r;
00173
00174 printf("%.2e\n", t_nfft);
00175
00176 fflush(stdout);
00177
00178 if(test_pre_full)
00179 nfft_free(A);
00180
00181 nfft_finalize(&np);
00182 }
00183
00184 int main(int argc,char **argv)
00185 {
00186 int l,trial;
00187
00188 if(argc<=2)
00189 {
00190 fprintf(stderr,"ndft_fast type first last trials\n");
00191 return -1;
00192 }
00193
00194 fprintf(stderr,"Testing ndft, Horner-like ndft, fully precomputed ndft.\n");
00195 fprintf(stderr,"Columns: N, M, t_ndft, t_horner, t_pre_full, t_nfft\n\n");
00196
00197
00198 if(atoi(argv[1])==0)
00199 {
00200 for(l=atoi(argv[2]); l<=atoi(argv[3]); l++)
00201 for(trial=0; trial<atoi(argv[4]); trial++)
00202 if(l<=13)
00203 ndft_time((1U<< l), (1U<< l), 1, 1);
00204 else
00205 if(l<=15)
00206 ndft_time((1U<< l), (1U<< l), 1, 0);
00207 else
00208 ndft_time((1U<< l), (1U<< l), 0, 0);
00209 }
00210
00211 return 1;
00212 }