| 1 | #include "sopnamsp.h"
|
|---|
| 2 | #include "machdefs.h"
|
|---|
| 3 |
|
|---|
| 4 | #include <math.h>
|
|---|
| 5 | #include <ctype.h>
|
|---|
| 6 | #include <iostream>
|
|---|
| 7 | #include <typeinfo>
|
|---|
| 8 |
|
|---|
| 9 | #include "srandgen.h"
|
|---|
| 10 | #include "matharr.h"
|
|---|
| 11 | #include "fftpserver.h"
|
|---|
| 12 | #include "fftmserver.h"
|
|---|
| 13 | #include "fftmayer.h"
|
|---|
| 14 | #include "fftwserver.h"
|
|---|
| 15 | #include "ntoolsinit.h"
|
|---|
| 16 |
|
|---|
| 17 | #include "timing.h"
|
|---|
| 18 |
|
|---|
| 19 | /* --------------------------------------------------------- */
|
|---|
| 20 | /* ---- Programme de test de calcul de FFT --- */
|
|---|
| 21 | /* ---- par les FFTServer de SOPHYA --- */
|
|---|
| 22 | /* 2000-2005 - C. Magneville, R. Ansari */
|
|---|
| 23 | /* Test FFT-1D real<>complex complex<>complex */
|
|---|
| 24 | /* Aide/Liste d'arguments : */
|
|---|
| 25 | /* csh> tfft */
|
|---|
| 26 | /* Tests de base ( ==> Rc=0 ) */
|
|---|
| 27 | /* csh> tfft 15 P */
|
|---|
| 28 | /* csh> tfft 16 P */
|
|---|
| 29 | /* csh> tfft 15 W */
|
|---|
| 30 | /* csh> tfft 16 W */
|
|---|
| 31 | /* csh> tfft 1531 P D 0 50 0.0002 */
|
|---|
| 32 | /* csh> tfft 1220 P D 0 50 0.0002 */
|
|---|
| 33 | /* csh> tfft 1531 W D 0 50 0.0002 */
|
|---|
| 34 | /* csh> tfft 1220 W D 0 50 0.0002 */
|
|---|
| 35 | /* --------------------------------------------------------- */
|
|---|
| 36 |
|
|---|
| 37 | static bool inp_typ_random = false ; // true -> random input
|
|---|
| 38 |
|
|---|
| 39 | template <class T>
|
|---|
| 40 | inline T module(complex<T> c)
|
|---|
| 41 | {
|
|---|
| 42 | return (sqrt(c.real()*c.real()+c.imag()*c.imag()));
|
|---|
| 43 | }
|
|---|
| 44 |
|
|---|
| 45 | // Max Matrix print elts
|
|---|
| 46 | static int nprt = 2;
|
|---|
| 47 | static int nprtfc = 8;
|
|---|
| 48 |
|
|---|
| 49 | static int prtlev = 0;
|
|---|
| 50 |
|
|---|
| 51 | //-------------------------------------------------
|
|---|
| 52 | template <class T>
|
|---|
| 53 | void TestFFTPack(T seuil, sa_size_t num)
|
|---|
| 54 | {
|
|---|
| 55 | int i;
|
|---|
| 56 | T fact = 1./num;
|
|---|
| 57 |
|
|---|
| 58 | TVector< complex<T> > inc(num), bkc(num), difc(num);
|
|---|
| 59 | TVector< T > in(num), ino(num), bk(num),dif(num);
|
|---|
| 60 | TVector< complex<T> > outc(num);
|
|---|
| 61 |
|
|---|
| 62 | cout << " DBG/1 outc " << outc.NElts() << endl;
|
|---|
| 63 | outc.ReSize(32);
|
|---|
| 64 | cout << " DBG/2 outc " << outc.NElts() << endl;
|
|---|
| 65 | outc.ReSize(10);
|
|---|
| 66 | cout << " DBG/3 outc " << outc.NElts() << endl;
|
|---|
| 67 | outc.ReSize(48);
|
|---|
| 68 | cout << " DBG/4 outc " << outc.NElts() << endl;
|
|---|
| 69 |
|
|---|
| 70 | if (inp_typ_random)
|
|---|
| 71 | for (i=0; i<num ; i++){
|
|---|
| 72 | ino[i] = in[i] = GaussianRand(1.,0.);
|
|---|
| 73 | inc[i] = complex<T> (in[i], 0.);
|
|---|
| 74 | }
|
|---|
| 75 | else for (i=0; i<num ; i++){
|
|---|
| 76 | ino[i] = in[i] = 0.5 + cos(2*M_PI*(double)i/(double)num)
|
|---|
| 77 | + 2*sin(4*M_PI*(double)i/(double)num);
|
|---|
| 78 | inc[i] = complex<T> (in[i], 0.);
|
|---|
| 79 | }
|
|---|
| 80 |
|
|---|
| 81 |
|
|---|
| 82 | cout << "Input / L = " << num << in;
|
|---|
| 83 | cout << endl;
|
|---|
| 84 |
|
|---|
| 85 | cout << " >>>> Testing FFTPackServer " << endl;
|
|---|
| 86 | FFTPackServer fftp;
|
|---|
| 87 | cout << " Testing FFTPackServer " << endl;
|
|---|
| 88 | fftp.fftf(in.NElts(), in.Data());
|
|---|
| 89 | // in /= (num/2.);
|
|---|
| 90 | cout << " fftp.fftf(in.NElts(), in.Data()) FORWARD: " << in << endl;
|
|---|
| 91 | cout << endl;
|
|---|
| 92 | fftp.fftb(in.NElts(), in.Data());
|
|---|
| 93 | cout << " fftp.fftb(in.NElts(), in.Data()) BACKWARD: " << in <<endl;
|
|---|
| 94 | cout << endl;
|
|---|
| 95 | dif = ino-in;
|
|---|
| 96 | cout << " dif , NElts= " << dif.NElts() << dif << endl;
|
|---|
| 97 |
|
|---|
| 98 | int ndiff = 0;
|
|---|
| 99 | T maxdif=0., vdif;
|
|---|
| 100 | for(i=0; i<num; i++) {
|
|---|
| 101 | vdif = fabs(dif(i));
|
|---|
| 102 | if (vdif > seuil) ndiff++;
|
|---|
| 103 | if (vdif > maxdif) maxdif = vdif;
|
|---|
| 104 | }
|
|---|
| 105 | cout << " Difference, Seuil= " << seuil << " NDiff= " << ndiff
|
|---|
| 106 | << " MaxDiff= " << maxdif << endl;
|
|---|
| 107 |
|
|---|
| 108 | }
|
|---|
| 109 |
|
|---|
| 110 | //-------------------------------------------------
|
|---|
| 111 | template <class T>
|
|---|
| 112 | int TestFFTS(T seuil, FFTServerInterface & ffts, sa_size_t num)
|
|---|
| 113 | {
|
|---|
| 114 |
|
|---|
| 115 | cout <<" ===> TestFFTS " << ffts.getInfo() << " ArrSz= " << num << endl;
|
|---|
| 116 | int i;
|
|---|
| 117 |
|
|---|
| 118 | T fact = 1.;
|
|---|
| 119 |
|
|---|
| 120 | TVector< complex<T> > inc(num), bkc(num), difc(num);
|
|---|
| 121 | TVector< T > in(num), ino(num), bk(num),dif(num);
|
|---|
| 122 | TVector< complex<T> > outc(num);
|
|---|
| 123 |
|
|---|
| 124 | if (inp_typ_random) {
|
|---|
| 125 | cout << " TestFFTS/Random input vector ... " << endl;
|
|---|
| 126 | ino = in = RandomSequence() ; // Nombre aleatoire gaussienne - sigma=1, mean=0
|
|---|
| 127 | ComplexMathArray< T > cma;
|
|---|
| 128 | TVector< T > im(num);
|
|---|
| 129 | im = RandomSequence(RandomSequence::Flat);
|
|---|
| 130 | inc = cma.FillFrom(in, im);
|
|---|
| 131 | }
|
|---|
| 132 | else {
|
|---|
| 133 | cout << " TestFFTS/Random input vector = 0.5+cos(2x)+2sin(4x)... " << endl;
|
|---|
| 134 | for (i=0; i<num ; i++){
|
|---|
| 135 | ino[i] = in[i] = 0.5 + cos(2*M_PI*(double)i/(double)num)
|
|---|
| 136 | + 2*sin(4*M_PI*(double)i/(double)num);
|
|---|
| 137 | inc[i] = complex<T> (in[i], 0.);
|
|---|
| 138 | }
|
|---|
| 139 | }
|
|---|
| 140 |
|
|---|
| 141 | cout << " Testing FFTServer " << ffts.getInfo() << endl;
|
|---|
| 142 |
|
|---|
| 143 | cout << "Input / Length= " << num << endl;
|
|---|
| 144 | if (prtlev > 0)
|
|---|
| 145 | cout << in << endl;
|
|---|
| 146 |
|
|---|
| 147 | int ndiff = 0;
|
|---|
| 148 | int rc = 0;
|
|---|
| 149 |
|
|---|
| 150 | cout << "\n ---- Testing FFT-1D(T, complex<T>) ---- " << endl;
|
|---|
| 151 | ffts.FFTForward(in, outc);
|
|---|
| 152 | if (prtlev > 0)
|
|---|
| 153 | cout << " FourierCoefs: outc= \n" << outc << endl;
|
|---|
| 154 |
|
|---|
| 155 | ffts.FFTBackward(outc, bk);
|
|---|
| 156 | if (prtlev > 0)
|
|---|
| 157 | cout << " Backward: bk= \n" << bk << endl;
|
|---|
| 158 |
|
|---|
| 159 | dif = bk*fact - in;
|
|---|
| 160 | if (prtlev > 0)
|
|---|
| 161 | cout << " Difference dif= \n" << dif << endl;
|
|---|
| 162 |
|
|---|
| 163 | ndiff = 0;
|
|---|
| 164 | T maxdif=0., vdif;
|
|---|
| 165 | for(i=0; i<num; i++) {
|
|---|
| 166 | vdif = fabs(dif(i));
|
|---|
| 167 | if (vdif > seuil) ndiff++;
|
|---|
| 168 | if (vdif > maxdif) maxdif = vdif;
|
|---|
| 169 | }
|
|---|
| 170 | cout << " Difference, Seuil= " << seuil << " NDiff= " << ndiff
|
|---|
| 171 | << " MaxDiff= " << maxdif << endl;
|
|---|
| 172 |
|
|---|
| 173 | if (ndiff != 0) rc += 4;
|
|---|
| 174 |
|
|---|
| 175 | cout << "\n ---- Testing FFT-1D(complex<T>, complex<T>) ---- " << endl;
|
|---|
| 176 | ffts.FFTForward(inc, outc);
|
|---|
| 177 | if (prtlev > 0)
|
|---|
| 178 | cout << " FourierCoef , outc= \n" << outc << endl;
|
|---|
| 179 |
|
|---|
| 180 | ffts.FFTBackward(outc, bkc);
|
|---|
| 181 | if (prtlev > 0)
|
|---|
| 182 | cout << " Backward , bkc= \n " << bkc << endl;
|
|---|
| 183 |
|
|---|
| 184 | difc = bkc*complex<T>(fact,0.) - inc;
|
|---|
| 185 | if (prtlev > 0)
|
|---|
| 186 | cout << " Difference , difc= \n " << difc << endl;
|
|---|
| 187 |
|
|---|
| 188 | ndiff = 0;
|
|---|
| 189 | maxdif=0.;
|
|---|
| 190 | for(i=0; i<num; i++) {
|
|---|
| 191 | vdif = fabs(module(difc(i)));
|
|---|
| 192 | if (vdif > seuil) ndiff++;
|
|---|
| 193 | if (vdif > maxdif) maxdif = vdif;
|
|---|
| 194 | }
|
|---|
| 195 | cout << " Difference, Seuil= " << seuil << " NDiff= " << ndiff
|
|---|
| 196 | << " MaxDiff= " << maxdif << endl;
|
|---|
| 197 |
|
|---|
| 198 | if (ndiff != 0) rc += 8;
|
|---|
| 199 | return rc;
|
|---|
| 200 | }
|
|---|
| 201 |
|
|---|
| 202 | inline void MayerFFTForw(r_4* d, int sz) {
|
|---|
| 203 | fht_r4(d, sz);
|
|---|
| 204 | }
|
|---|
| 205 | inline void MayerFFTForw(r_8* d, int sz) {
|
|---|
| 206 | fht_r8(d, sz);
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 | //-------------------------------------------------
|
|---|
| 210 | template <class T>
|
|---|
| 211 | int MultiFFTTest(T t, sa_size_t sz, int nloop, bool fgp)
|
|---|
| 212 | {
|
|---|
| 213 | cout <<" ===> MultiFFTTest<T=r_" << sizeof(T) << "> NLoop= " << nloop << " ArrSz= " << sz << endl;
|
|---|
| 214 | TVector< T > in(sz), incopie(sz);
|
|---|
| 215 |
|
|---|
| 216 | incopie = in = RandomSequence();
|
|---|
| 217 | FFTPackServer fftp(false);
|
|---|
| 218 | TVector< complex<T> > outc;
|
|---|
| 219 |
|
|---|
| 220 | PrtTim("MultiFFT-LoopStart");
|
|---|
| 221 |
|
|---|
| 222 | if (fgp) {
|
|---|
| 223 | cout << " --- Test effectue avec FFTPack " << endl;
|
|---|
| 224 | for(int kk=0; kk<nloop; kk++) {
|
|---|
| 225 | in = incopie;
|
|---|
| 226 | // fftp.FFTForward(in, outc);
|
|---|
| 227 | fftp.fftf(sz, in.Data());
|
|---|
| 228 | }
|
|---|
| 229 | }
|
|---|
| 230 | else {
|
|---|
| 231 | cout << " --- Test effectue avec FFTMayer " << endl;
|
|---|
| 232 | for(int kk=0; kk<nloop; kk++) {
|
|---|
| 233 | in = incopie;
|
|---|
| 234 | // in += 0.01;
|
|---|
| 235 | MayerFFTForw(in.Data(), sz);
|
|---|
| 236 | }
|
|---|
| 237 | }
|
|---|
| 238 | PrtTim("MultiFFT-LoopEnd");
|
|---|
| 239 | cout << " ----- End ---- MultiFFTTest<T> ----" << endl;
|
|---|
| 240 | return 0;
|
|---|
| 241 | }
|
|---|
| 242 |
|
|---|
| 243 |
|
|---|
| 244 |
|
|---|
| 245 | //-------------------------------------------------
|
|---|
| 246 | //-------------------------------------------------
|
|---|
| 247 | //-------------------------------------------------
|
|---|
| 248 | int main(int narg, char* arg[])
|
|---|
| 249 | {
|
|---|
| 250 |
|
|---|
| 251 | SophyaInit();
|
|---|
| 252 | InitTim(); // Initializing the CPU timer
|
|---|
| 253 |
|
|---|
| 254 | if (narg < 3) {
|
|---|
| 255 | cout << "tfft/ args error - \n Usage tfft size p/P/M/W [f/d/F/D PrtLev=0 MaxNPrt=50 diffthr] \n"
|
|---|
| 256 | << " OR tfft size P/M FL/DL/FZ/DZ [NFFT=10] \n"
|
|---|
| 257 | << " 1D real/complex-FFT test (FFTServer) \n"
|
|---|
| 258 | << " size: input vector length \n "
|
|---|
| 259 | << " p=FFTPackTest P=FFTPack, M=FFTMayer, W= FFTWServer \n "
|
|---|
| 260 | << " F/f:float, D/d:double F/D:random in_vector (default=D) \n"
|
|---|
| 261 | << " FL/DL : perform NFFT on a serie of float[size] or double[size] "
|
|---|
| 262 | << " FZ/DZ : perform NFFT on a series of complex<float>[size] or complex<double>[size] "
|
|---|
| 263 | << " diffthr : Threshold for diff checks (=10^-6/10^-4 double/float)"
|
|---|
| 264 | << endl;
|
|---|
| 265 | return(1);
|
|---|
| 266 | }
|
|---|
| 267 |
|
|---|
| 268 | FFTPackServer fftp;
|
|---|
| 269 | FFTMayerServer fftm;
|
|---|
| 270 |
|
|---|
| 271 | FFTWServer fftw;
|
|---|
| 272 |
|
|---|
| 273 | sa_size_t sz = atol(arg[1]);
|
|---|
| 274 | int nprt = 50;
|
|---|
| 275 | if (narg > 4) prtlev = atoi(arg[4]);
|
|---|
| 276 | if (narg > 5) nprt = atoi(arg[5]);
|
|---|
| 277 | BaseArray::SetMaxPrint(nprt, prtlev);
|
|---|
| 278 | float fs = 1.e-4;
|
|---|
| 279 | double ds = 1.e-6;
|
|---|
| 280 | if (narg > 6) fs = ds = atof(arg[6]);
|
|---|
| 281 |
|
|---|
| 282 | if (sz < 2) sz = 2;
|
|---|
| 283 |
|
|---|
| 284 | char dtyp = 'D';
|
|---|
| 285 | if (narg > 3) dtyp = *arg[3];
|
|---|
| 286 | inp_typ_random = true;
|
|---|
| 287 | if (islower(dtyp)) inp_typ_random = false;
|
|---|
| 288 | dtyp = toupper(dtyp);
|
|---|
| 289 |
|
|---|
| 290 |
|
|---|
| 291 | FFTServerInterface * ffts;
|
|---|
| 292 | if (*arg[2] == 'M') ffts = fftm.Clone();
|
|---|
| 293 | else if (*arg[2] == 'W') ffts = fftw.Clone();
|
|---|
| 294 | else ffts = fftp.Clone();
|
|---|
| 295 |
|
|---|
| 296 | int nloop = 0;
|
|---|
| 297 | bool fgpack = true;
|
|---|
| 298 | if (arg[3][1] == 'L') { // MultiFFT test
|
|---|
| 299 | nloop = 10;
|
|---|
| 300 | if (narg > 4) nloop = atoi(arg[4]);
|
|---|
| 301 | if (*arg[2] == 'M') fgpack = false;
|
|---|
| 302 | }
|
|---|
| 303 |
|
|---|
| 304 |
|
|---|
| 305 | cout << "\n ============================================= \n"
|
|---|
| 306 | << " ------ Testing FFTServer " << typeid(*ffts).name() << "\n"
|
|---|
| 307 | << " VecSize= " << sz << " InputType= "
|
|---|
| 308 | << ( (inp_typ_random ) ? " random " : " fixed(sin+cos) ") << "\n"
|
|---|
| 309 | << " =============================================== " << endl;
|
|---|
| 310 |
|
|---|
| 311 | int rc = 0;
|
|---|
| 312 | try {
|
|---|
| 313 | if (dtyp == 'D') {
|
|---|
| 314 | cout << " ------ Testing FFTServer for double (r_8)----- " << endl;
|
|---|
| 315 | if (*arg[2] == 'p') TestFFTPack(ds, sz);
|
|---|
| 316 | else if (nloop > 0) rc = MultiFFTTest(ds, sz, nloop, fgpack);
|
|---|
| 317 | else rc = TestFFTS(ds, *ffts, sz);
|
|---|
| 318 | }
|
|---|
| 319 | else {
|
|---|
| 320 | cout << " ------ Testing FFTServer for float (r_4)----- " << endl;
|
|---|
| 321 | if (*arg[2] == 'p') TestFFTPack(fs, sz);
|
|---|
| 322 | else if (nloop > 0) rc = MultiFFTTest(fs, sz, nloop, fgpack);
|
|---|
| 323 | else rc = TestFFTS(fs, *ffts, sz);
|
|---|
| 324 | }
|
|---|
| 325 | }
|
|---|
| 326 | catch(PThrowable& exc ) {
|
|---|
| 327 | cerr << "TestFFT-main() , Catched exception: \n" << exc.Msg() << endl;
|
|---|
| 328 | rc = 97;
|
|---|
| 329 | }
|
|---|
| 330 | catch(std::exception ex) {
|
|---|
| 331 | cerr << "TestFFT-main() , Catched exception ! " << (string)(ex.what()) << endl;
|
|---|
| 332 | rc = 98;
|
|---|
| 333 | }
|
|---|
| 334 | catch(...) {
|
|---|
| 335 | cerr << "TestFFT-main() , Catched ... exception ! " << endl;
|
|---|
| 336 | rc = 99;
|
|---|
| 337 | }
|
|---|
| 338 | PrtTim("End of tfft ");
|
|---|
| 339 | delete ffts;
|
|---|
| 340 | cout << " =========== End of tfft Rc= " << rc << " =============" << endl;
|
|---|
| 341 | }
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
