1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
|
// Test1.cpp
// Test Declaration of Variable
// Test Basic Math
// Test using Libraries
// Test rand and FFT
// Test funtion input and output, calling other functions
//
//
// Copyright 2010
// Sebastian Dietze
// University of California, San Diego
//Date: 2013/09/05
#include <stdio.h> /* printf, scanf, puts, NULL */
#include <stdlib.h> /* srand, rand */
#include "cycle.h"
#include <ctime>
#include <time.h> /* clock_t, clock, CLOCKS_PER_SEC */
#include "mtrand.h" /* Mersenne Twister PRNG must be installed*/
#include <cstdio>
#include <math.h> /* pow, exp, etc */
#include <complex> /* Complex Numbers */
#include <exception>
#include <new>
extern "C" { /* FFTW3.3 must be installed */
#include <fftw3.h> /* Link using -lfftw3f -lm for floating point operations */
} /* define this after <complex> */
using namespace std;
//Define decleration variable
typedef std::complex<float> csingle;
typedef std::complex<double> cdouble;
//function declerations
void doLoop(int, csingle[], int[]);
void makeArray(csingle[], int[], float, float);
const double pi = std::acos(-1);
int main()
{
int theSize[3] = {5,3,4}, n = 1000000;
long Npts = *(theSize+0)* *(theSize+1)* *(theSize+2);
float theAmp = 20.0, thePhase = 2*pi;
csingle* theArray;
try{
theArray = new csingle[Npts]; //dynamic creation of array referenced with pointer
}
catch (std::bad_alloc){
printf("Error allocating memory.");
return 0;
}
//pass the array using pointer
makeArray(theArray,theSize,theAmp,thePhase);
for(int ii = 0; ii<10; ii++){
printf("Before %+08.4f %+08.4f i\n",std::real(*(theArray+ii)),std::imag(*(theArray+ii)));
}
printf("\n");
printf("performing 2*%d FFTs\n",n);
std::clock_t t1 = std::clock();
doLoop(n, theArray, theSize);
std::clock_t t2 = std::clock();
for(int ii = 0; ii<10; ii++){
printf("After %+08.4f %+08.4f i\n",std::real(*(theArray+ii)),std::imag(*(theArray+ii)));
}
printf("\n");
printf("Elapse Time: %f\n",(t2 - t1 ) / (double) CLOCKS_PER_SEC);
delete[] theArray;
return 0;
}
void doLoop(int theCount, csingle theArray[], int theSize[])
{
long Npts = *(theSize+0)* *(theSize+1)* *(theSize+2);
fftwf_plan pf, pr;
pf = fftwf_plan_dft(3,theSize, reinterpret_cast<fftwf_complex*>(theArray), reinterpret_cast<fftwf_complex*>(theArray), FFTW_FORWARD, FFTW_ESTIMATE);
pr = fftwf_plan_dft(3,theSize, reinterpret_cast<fftwf_complex*>(theArray), reinterpret_cast<fftwf_complex*>(theArray), FFTW_BACKWARD, FFTW_ESTIMATE);
for(int ii = 0; ii<theCount; ii++)
{
//FFT
fftwf_execute(pf);
//iFFT
fftwf_execute(pr);
for(int jj = 0; jj<Npts; jj++){
*(theArray + jj) = *(theArray + jj)/float(Npts);
}
}
fftwf_destroy_plan(pf);
fftwf_destroy_plan(pr);
//return 0;
}
void makeArray(csingle theArray[], int theSize[], float theAmp, float thePhase)
{
unsigned long t = time(NULL);
int ind;
/* initialize random seed: */
MTRand drand(&t, 1);//Mersenne Twister
for(int kk = 0; kk<*(theSize+2); kk++){
for(int jj = 0; jj<*(theSize+1); jj++){
for(int ii = 0; ii<*(theSize+0); ii++){
ind = (kk* *(theSize+2) + jj)* *(theSize+0) + ii;
*(theArray+ind) = (csingle)std::polar(theAmp*drand(),thePhase*drand());
}}}
return;
}
|