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