ese 519
Revision 0:2076b4d80327, committed 2015-04-07
- Comitter:
- niv17
- Date:
- Tue Apr 07 21:09:22 2015 +0000
- Commit message:
- sonic initial april 7
Changed in this revision
diff -r 000000000000 -r 2076b4d80327 AudioObj.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/AudioObj.cpp Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,105 @@ +#include "AudioObj.h" +//#include <iostream> + +Location AudioObj::getLocation () const { + return this->location; +} + +void AudioObj::setLocation (const Location& loc) { + this->location = loc; +} + +void AudioObj::setLocation (float x, float y, float z) { + this->location = Location(x,y,z); +} + +Velocity AudioObj::getVelocity () const { + return this->velocity; +} + +void AudioObj::setVelocity (const Velocity& vel) { + this->velocity = vel; +} + +void AudioObj::setVelocity (float dx, float dy, float dz) { + this->velocity = Velocity(dx,dy,dz); +} + +float AudioObj::getVolume() const { + return this->volume; +} + +void AudioObj::setVolume(float vol) { + if (vol > 1 || vol < 0){ + std::cout<<"Volume not in range (0-1)"<<endl; + //throw invalid_argument("Volume not in range (0-1)"); + } + this->volume = vol; +} + +void AudioObj::setRandomVolume() { + float randVol = (rand() % 100 + 1) / 100.0; + this->setVolume(randVol); +} + +bool AudioObj::isActive() const { + return this->active; +} + +void AudioObj::setActive(bool active){ + this->active = active; +} + +bool AudioObj::isGpsObject() const { + return this->gpsObject; +} + +bool AudioObj::isBackgroundObject() const { + return this->backgroundObject; +} + +void AudioObj::setRepeat(bool rep) { + this->repeat = rep; +} + + +bool AudioObj::fillAudioData (Complex* target, unsigned int length) { + if(circularBuffer.readSizeRemaining() < length) { + return false; + } + circularBuffer.read(target, length); + return true; +} + +void AudioObj::loadCircularBuffer() { + unsigned int length = (unsigned int)circularBuffer.writeSizeRemaining(); + if(length>16384) { // TODO: How was this number chosen??? + //cout<<"In write circ Buff : "<<length<<endl; + if(!(wavObject.loadMoreData(length, repeat))) { + // end of file reached, no repeat + this->active = false; + this->isCompleted = true; + return; + } + circularBuffer.write(wavObject.complexTempData, length); + } +} + +void AudioObj::playOnceFromBeginning() { + this->repeat = false; + this->restart(); +} + +void AudioObj::restart() { + this->circularBuffer.clear(); + this->wavObject.seekToBeginning(); + this->loadCircularBuffer(); + this->active = true; +} + +void AudioObj::loadFromBeginning() { + this->circularBuffer.clear(); + this->wavObject.seekToBeginning(); + this->loadCircularBuffer(); + this->active = false; +}
diff -r 000000000000 -r 2076b4d80327 Mixer3D.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Mixer3D.cpp Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,298 @@ +#include <math.h> +#include "Mixer3D.h" +//#include "World.h" +#include "fft.h" +#include "mit_hrtf_lib.h" +#include <sys/time.h> + +Mixer3D::Mixer3D(int bufSize, int smpRate, int bitD) : +bufferSize(bufSize), sampleRate(smpRate), bitDepth(bitD)) +{ + cout << "BUF_SIZE: " << bufSize << endl; + cout << "SAMPLE_RATE: " << smpRate << endl; + cout << "BIT_DEPTH: " << bitD << endl; + inputAO = new Complex[2*bufferSize]; + overlapInput = new Complex[2 * bufferSize]; + fInput = new Complex[2 * bufferSize]; + fFilter = new Complex[2 * bufferSize]; + + azimuths = new int[World::MAX_OBJ]; + elevations = new int[World::MAX_OBJ]; + prevAzimuths = new int[World::MAX_OBJ]; + prevElevations = new int[World::MAX_OBJ]; + + updateAngles(); // initialize azimuths and elevations + + outputLeft = new Complex*[World::MAX_OBJ]; + outputRight = new Complex*[World::MAX_OBJ]; + overlapLeft = new Complex*[World::MAX_OBJ]; + overlapRight = new Complex*[World::MAX_OBJ]; + + leftFilter = new short[2 * bufferSize]; + rightFilter = new short[2 * bufferSize]; + + +} + +Mixer3D::~Mixer3D() +{ + delete [] inputAO; + delete [] outputLeft; + delete [] outputLeft; + delete [] leftFilter; + delete [] rightFilter; + delete [] complexLeftFilter; + delete [] complexRightFilter; + delete [] overlapLeft; + delete [] overlapRight; + delete [] overlapInput; + delete [] prevAzimuths; + delete [] prevElevations; + delete [] azimuths; + delete [] elevations; +} + +void Mixer3D::updateAngles() { + AudioObj* iAudioObj; + Location iLocation; + for (int i = 0; i < myWorld->getNumObj(); i++) { + iAudioObj = myWorld->getAudioObj(i); + iLocation = iAudioObj->getLocation(); + + if (iAudioObj->isGpsObject()) { + azimuths[i] = player.computeAzimuthFromGpsCoordinates(&iLocation); + elevations[i] = player.computeElevationFromGpsCoordinates(&iLocation); + } else { + azimuths[i] = player.computeAzimuth(&iLocation); + elevations[i] = player.computeElevation(&iLocation); + } + } +} + +bool Mixer3D::isPowerOfTwo(int x) { + return !(x == 0) && !(x & (x - 1)); +} + +int Mixer3D::loadHRTF(int* pAzimuth, int* pElevation, unsigned int samplerate, unsigned int diffused, Complex *&leftFilterIn, Complex *&rightFilterIn) +{ + int size = mit_hrtf_get(pAzimuth, pElevation, samplerate, diffused, leftFilter, rightFilter); + + if (size == 0) { + // TODO: Throw MIT HRTF filter does not exist exception + } + for (int i = 0; i < size; i++) + { + leftFilterIn[i] = (double)(leftFilter[i]); + rightFilterIn[i] = (double)(rightFilter[i]); + } + + return size; +} +void Mixer3D::convolution(Complex *input, Complex *filter, Complex *output, long nSig, long nFil, long nFFT) { + + if (input == NULL || filter == NULL) { + throw invalid_argument("Input and Filter must be non-NULL."); + } + + if (!isPowerOfTwo(nFFT) || nFFT < nSig || nFFT < nFil) { + throw invalid_argument("NFFT must be a power of two, bigger than the signal length, and bigger than the filter length."); + } + + // Perform FFT on both input and filter. + // TODO: "Streamline" CFFT class? + CFFT::Forward(input, fInput, (unsigned int)nFFT); + CFFT::Forward(filter, fFilter, (unsigned int)nFFT); + + for (int i = 0; i < nFFT; i++) { + output[i] = fInput[i] * fFilter[i]; + } + CFFT::Inverse(output, (unsigned int)nFFT); +} + +void Mixer3D::stereoConvolution(Complex *input, Complex *leftFilter, Complex *rightFilter, Complex *leftOutput, Complex *rightOutput, long nSIG, long nFIL, long nFFT) +{ + // TODO: Modify parameter name, input is a data member + convolution(input, leftFilter, leftOutput, nSIG, nFIL, nFFT); + convolution(input, rightFilter, rightOutput, nSIG, nFIL, nFFT); +} + +void Mixer3D::computeValidAudioObjects(vector<AudioObj*> &validAudioObjectsOut) { + AudioObj* iAudioObj; + float iDistance; + Location iLocation; + for (int i=0; i < myWorld->getNumObj(); i++) { + iAudioObj = myWorld->getAudioObj(i); + if (iAudioObj->isActive()) { + iLocation = iAudioObj->getLocation(); + if (iAudioObj->isGpsObject()) { + if (player.isTrackingGps()) { + iDistance = player.computeDistanceFromGpsCoordinates(&iLocation); + } else { + // not tracking gps, so skip + continue; + } + } else { + iDistance = player.computeDistanceFrom(&iLocation); + } + + if (iDistance < MAX_GPS_OBJ_DISTANCE) { + validAudioObjectsOut.push_back(iAudioObj); + } + } + } +} + +void Mixer3D::performMix(short *ioDataLeft, short *ioDataRight) +{ + double runtime; + struct timeval start, finish; + gettimeofday(&start, NULL); + //Zero the output arrays. + for(int i = 0; i < bufferSize; i++) + { + ioDataLeft[i] = 0; + ioDataRight[i] = 0; + } + + //Iterate through all audio objects, obtain input data and calculate resulting audio data for each object. + AudioObj* iAudioObj; + Location iLocation; + float iVolume, iDistance, iAmplitudeFactor; + vector<AudioObj*> validAudioObjects; + + // get valid audio objects + computeValidAudioObjects(validAudioObjects); + int nValidAudioObjects = (int)validAudioObjects.size(); + for (int i = 0; i < nValidAudioObjects; i++) { + iAudioObj = validAudioObjects[i]; + + // loading in input data for the iteration accordingly + if (!(iAudioObj->fillAudioData(inputAO, bufferSize))) { + cout << "Audio Object" << i << " not loading input quickly enough" << endl; + continue; + } + + iVolume = iAudioObj->getVolume(); + if (iVolume == 0) { + // skip 'muted' audio objects, but not after loading their input. + // that way they are still playing, but silently. + continue; + } + + // handle background objects + if (iAudioObj->isBackgroundObject()) { + // copy input directly to output buffers + for (int j=0; j < bufferSize; j++) { + // scale input from [-1, 1] to proper scale depending on bit depth + ioDataLeft[j] += inputAO[j].real()*iVolume*pow(2, bitDepth) / nValidAudioObjects; + ioDataRight[j] += inputAO[j].real()*iVolume*pow(2, bitDepth) / nValidAudioObjects; + } + continue; + } + + iLocation = iAudioObj->getLocation(); + + // handle gps objects + if (iAudioObj->isGpsObject()) { + if (player.isTrackingGps()) { + iDistance = player.computeDistanceFromGpsCoordinates(&iLocation); + if (iDistance > MAX_GPS_OBJ_DISTANCE) { + // not within range, skip this audio object + continue; + } + } else { + // gps tracking hasn't started yet + continue; + } + } else { + iDistance = player.computeDistanceFrom(&iLocation); + } + + if (iDistance == 0) { + iAmplitudeFactor = iVolume; + } else { + // ensure that the amplitude factor can never be > 1 to prevent clipping + if (iDistance > 1) { + iAmplitudeFactor = iVolume / pow(iDistance, 2); + } else { + iAmplitudeFactor = iVolume; + } + } + + // dynamically calculate the Azimuth and Elevation between every object and the player + updateAngles(); + + for(int j = 0; j < bufferSize * 2; j++) { + if ( j >= bufferSize ) { + // zero pad + inputAO[j] = 0; + } else { + inputAO[j] *= iAmplitudeFactor; + } + } + + if (azimuths[i] != prevAzimuths[i] || + elevations[i] != prevElevations[i]) { + // object location relative to player has changed, so fetch a new filter + filterLength = loadHRTF(&azimuths[i], &elevations[i], sampleRate, 1, complexLeftFilter[i], complexRightFilter[i]); + + // zero pad + for (int j = filterLength; j < 2 * bufferSize; j++) { + complexLeftFilter[i][j] = 0; + complexRightFilter[i][j] = 0; + } + + if (azimuths[i] < 0) { + azimuths[i] = -azimuths[i]; + } + + + //Since the filter changed, we perform a convolution on the old input with the new filter and pull out its tail. + stereoConvolution(overlapInput, complexLeftFilter[i], complexRightFilter[i], outputLeft[i], outputRight[i], bufferSize, filterLength, 2 * bufferSize); + + // update the overlap part for the next iteration + for (int j = 0; j < bufferSize; j++) { + overlapLeft[i][j] = outputLeft[i][j + bufferSize]; + overlapRight[i][j] = outputRight[i][j + bufferSize]; + } + } + + //Perform the convolution of the current input and current filter. + stereoConvolution(inputAO, complexLeftFilter[i], complexRightFilter[i], outputLeft[i], outputRight[i], bufferSize, filterLength, 2 * bufferSize); + + //Output the data to the output arrays in short integer format. In addition to pushing the output of the main + //convolution, we also need to add the overlapped tail of the last output and divide by 2. Finally, we need + //to divide by the number of audio objects to ensure no clipping. + for (int j = 0; j < bufferSize; j++) + { + ioDataLeft[j] += (short)( (outputLeft[i][j].real() + overlapLeft[i][j].real()) / (2*nValidAudioObjects)); + ioDataRight[j] += (short)( (outputRight[i][j].real() + overlapRight[i][j].real()) / (2*nValidAudioObjects)); + } + + //Updating the overlapInput for the next iteration for the correpsonding object + for (int j = 0; j < bufferSize * 2; j++) { + if (j >= bufferSize) { + overlapInput[j] = 0; + } else { + overlapInput[j] = inputAO[j]; + } + } + + // TODO: If the filter has been changed, didn't we already do this? + //Updating the default overlap information for the next iteration if the filter won't be changed + for (int j = 0 ; j < bufferSize; j++) { + overlapLeft[i][j] = outputLeft[i][j + bufferSize]; + overlapRight[i][j] = outputRight[i][j + bufferSize]; + } + + //storing the Azimuth value in this iteration for the comparison for the next iteration so that + //we can know that whether the filter needs to be changed in the next iteration. + prevAzimuths[i] = azimuths[i]; + prevElevations[i] = elevations[i]; + } + gettimeofday(&finish, NULL); + runtime = finish.tv_sec - start.tv_sec; + runtime += (finish.tv_usec - start.tv_usec) / 1000000.0; + + cout << "performMix: " << runtime << endl; +}
diff -r 000000000000 -r 2076b4d80327 fft.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fft.cpp Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,578 @@ +// fft.cpp - impelementation of class +// of fast Fourier transform - FFT +// +// The code is property of LIBROW +// You can use it on your own +// When utilizing credit LIBROW site + +// modified by Da Cao, Oct.1st/2014. +// Changed into a FHT method to compute DFT + +#include <fstream> +#include <iostream> +#include <math.h> +#include "fft.h" +//#include "../include/complex.h" + +#include <cmath> + +//#include "revbinpermute.h" +#include "../include/fhtloc2.h" +#include "../include/fht.h" +#include "../include/cmult.h" +#include "../include/constants.h" +#include "../include/fxttypes.h" +#include "../include/complextype.h" +#include "../include/revbinpermute.h" +#include "../Sonic.h" // getPermutedIndex needs BUF_SIZE + +using namespace std; + +//Convolution function which returns the frequency representation of the result. +//NFFT is the FFT size, nSIG is the size for the input, NFIL is the size of the filter. + +//fft_dit4_core_p1 define method +// trig data used in method + +unsigned int getPermutedIndex (unsigned int num) { + unsigned int NO_OF_BITS = sizeof(num) * 8; + unsigned int reverse_num = 0, i, temp; + + for (i = 0; i < NO_OF_BITS; i++) { + temp = (num & (1 << i)); + if(temp) { + reverse_num |= (1 << ((NO_OF_BITS - 1) - i)); + } + } + return (reverse_num >> (NO_OF_BITS - (int)log2(BUF_SIZE))); +} + +void +fft8_dit_core_p1(double *fr, double *fi) +// 8-point decimation in time FFT +// isign = +1 +// input data must be in revbin_permuted order +//. +// Cf. Nussbaumer p.148f +{ + // INPUT_RE: + double t1r = fr[0] + fr[1]; + double t2r = fr[2] + fr[3]; + double t7r = t1r + t2r; + double t3r = fr[4] - fr[5]; + double t4r = fr[4] + fr[5]; + double t5r = fr[6] + fr[7]; + double t8r = t4r + t5r; + double t6r = fr[6] - fr[7]; + + double m0r = t7r + t8r; + double m1r = t7r - t8r; + double m2r = t1r - t2r; + double m3r = fr[0] - fr[1]; + double m4r = M_SQRT1_2 * (t3r - t6r); + +#define m5i t6r +#define m6i t7r +#define m7i t8r + m7i = M_SQRT1_2 * (t3r + t6r); + m5i = t5r - t4r; + m6i = fr[3] - fr[2]; + + // INPUT_IM: + double t1i = fi[0] + fi[1]; + double t2i = fi[2] + fi[3]; + double t7i = t1i + t2i; + double t3i = fi[4] - fi[5]; + double t4i = fi[4] + fi[5]; + double t5i = fi[6] + fi[7]; + double t8i = t4i + t5i; + double t6i = fi[6] - fi[7]; + + double m0i = t7i + t8i; + double m1i = t7i - t8i; + double m2i = t1i - t2i; + double m3i = fi[0] - fi[1]; + double m4i = M_SQRT1_2 * (t3i - t6i); + +#define m5r t6i +#define m6r t7i +#define m7r t8i + m7r = M_SQRT1_2 * (t3i + t6i); + m5r = t4i - t5i; + m6r = fi[2] - fi[3]; + +#define s1r t1r +#define s2r t2r +#define s3r t3r +#define s4r t4r + s1r = m3r + m4r; + s2r = m3r - m4r; + s3r = m6r + m7r; + s4r = m6r - m7r; + + // OUTPUT_RE: + fr[0] = m0r; + fr[7] = s1r + s3r; + fr[6] = m2r + m5r; + fr[5] = s2r - s4r; + fr[4] = m1r; + fr[3] = s2r + s4r; + fr[2] = m2r - m5r; + fr[1] = s1r - s3r; + +#define s1i t1r +#define s2i t2r +#define s3i t3r +#define s4i t4r + s1i = m3i + m4i; + s2i = m3i - m4i; + s3i = m6i - m7i; + s4i = m6i + m7i; + + // OUTPUT_IM: + fi[0] = m0i; + fi[7] = s1i + s3i; + fi[6] = m2i + m5i; + fi[5] = s2i - s4i; + fi[4] = m1i; + fi[3] = s2i + s4i; + fi[2] = m2i - m5i; + fi[1] = s1i - s3i; +} +// ------------------------- + +#undef s1r +#undef s2r +#undef s3r +#undef s4r + +#undef s1i +#undef s2i +#undef s3i +#undef s4i + +#undef m5r +#undef m6r +#undef m7r + +#undef m5i +#undef m6i +#undef m7i +// ------------------------- +////////////////////////////////////////////////////////////////////////////////////////////////// + + +static const ulong LX = 2; + +double allTheC [680] = {1, 0.980785, 0.92388, 0.83147, 0.707107, 0.55557, 0.382683, 0.19509, 1, 0.998795, 0.995185, 0.989177, 0.980785, 0.970031, 0.95694, 0.941544, 0.92388, 0.903989, 0.881921, 0.857729, 0.83147, 0.803208, 0.77301, 0.740951, 0.707107, 0.671559, 0.634393, 0.595699, 0.55557, 0.514103, 0.471397, 0.427555, 0.382683, 0.33689, 0.290285, 0.24298, 0.19509, 0.14673, 0.0980171, 0.0490677, 1, 0.999925, 0.999699, 0.999322, 0.998795, 0.998118, 0.99729, 0.996313, 0.995185, 0.993907, 0.99248, 0.990903, 0.989177, 0.987301, 0.985278, 0.983105, 0.980785, 0.978317, 0.975702, 0.97294, 0.970031, 0.966976, 0.963776, 0.960431, 0.95694, 0.953306, 0.949528, 0.945607, 0.941544, 0.937339, 0.932993, 0.928506, 0.92388, 0.919114, 0.91421, 0.909168, 0.903989, 0.898674, 0.893224, 0.88764, 0.881921, 0.87607, 0.870087, 0.863973, 0.857729, 0.851355, 0.844854, 0.838225, 0.83147, 0.824589, 0.817585, 0.810457, 0.803208, 0.795837, 0.788346, 0.780737, 0.77301, 0.765167, 0.757209, 0.749136, 0.740951, 0.732654, 0.724247, 0.715731, 0.707107, 0.698376, 0.689541, 0.680601, 0.671559, 0.662416, 0.653173, 0.643832, 0.634393, 0.624859, 0.615232, 0.605511, 0.595699, 0.585798, 0.575808, 0.565732, 0.55557, 0.545325, 0.534998, 0.52459, 0.514103, 0.503538, 0.492898, 0.482184, 0.471397, 0.460539, 0.449611, 0.438616, 0.427555, 0.41643, 0.405241, 0.393992, 0.382683, 0.371317, 0.359895, 0.348419, 0.33689, 0.32531, 0.313682, 0.302006, 0.290285, 0.27852, 0.266713, 0.254866, 0.24298, 0.231058, 0.219101, 0.207111, 0.19509, 0.18304, 0.170962, 0.158858, 0.14673, 0.134581, 0.122411, 0.110222, 0.0980171, 0.0857973, 0.0735646, 0.0613207, 0.0490677, 0.0368072, 0.0245412, 0.0122715, 1, 0.999995, 0.999981, 0.999958, 0.999925, 0.999882, 0.999831, 0.999769, 0.999699, 0.999619, 0.999529, 0.999431, 0.999322, 0.999205, 0.999078, 0.998941, 0.998795, 0.99864, 0.998476, 0.998302, 0.998118, 0.997925, 0.997723, 0.997511, 0.99729, 0.99706, 0.99682, 0.996571, 0.996313, 0.996045, 0.995767, 0.995481, 0.995185, 0.994879, 0.994565, 0.99424, 0.993907, 0.993564, 0.993212, 0.99285, 0.99248, 0.992099, 0.99171, 0.991311, 0.990903, 0.990485, 0.990058, 0.989622, 0.989177, 0.988722, 0.988258, 0.987784, 0.987301, 0.986809, 0.986308, 0.985798, 0.985278, 0.984749, 0.98421, 0.983662, 0.983105, 0.982539, 0.981964, 0.981379, 0.980785, 0.980182, 0.97957, 0.978948, 0.978317, 0.977677, 0.977028, 0.97637, 0.975702, 0.975025, 0.974339, 0.973644, 0.97294, 0.972226, 0.971504, 0.970772, 0.970031, 0.969281, 0.968522, 0.967754, 0.966976, 0.96619, 0.965394, 0.96459, 0.963776, 0.962953, 0.962121, 0.96128, 0.960431, 0.959572, 0.958703, 0.957826, 0.95694, 0.956045, 0.955141, 0.954228, 0.953306, 0.952375, 0.951435, 0.950486, 0.949528, 0.948561, 0.947586, 0.946601, 0.945607, 0.944605, 0.943593, 0.942573, 0.941544, 0.940506, 0.939459, 0.938404, 0.937339, 0.936266, 0.935184, 0.934093, 0.932993, 0.931884, 0.930767, 0.929641, 0.928506, 0.927363, 0.92621, 0.925049, 0.92388, 0.922701, 0.921514, 0.920318, 0.919114, 0.917901, 0.916679, 0.915449, 0.91421, 0.912962, 0.911706, 0.910441, 0.909168, 0.907886, 0.906596, 0.905297, 0.903989, 0.902673, 0.901349, 0.900016, 0.898674, 0.897325, 0.895966, 0.894599, 0.893224, 0.891841, 0.890449, 0.889048, 0.88764, 0.886223, 0.884797, 0.883363, 0.881921, 0.880471, 0.879012, 0.877545, 0.87607, 0.874587, 0.873095, 0.871595, 0.870087, 0.868571, 0.867046, 0.865514, 0.863973, 0.862424, 0.860867, 0.859302, 0.857729, 0.856147, 0.854558, 0.852961, 0.851355, 0.849742, 0.84812, 0.846491, 0.844854, 0.843208, 0.841555, 0.839894, 0.838225, 0.836548, 0.834863, 0.83317, 0.83147, 0.829761, 0.828045, 0.826321, 0.824589, 0.82285, 0.821103, 0.819348, 0.817585, 0.815814, 0.814036, 0.812251, 0.810457, 0.808656, 0.806848, 0.805031, 0.803208, 0.801376, 0.799537, 0.797691, 0.795837, 0.793975, 0.792107, 0.79023, 0.788346, 0.786455, 0.784557, 0.782651, 0.780737, 0.778817, 0.776888, 0.774953, 0.77301, 0.771061, 0.769103, 0.767139, 0.765167, 0.763188, 0.761202, 0.759209, 0.757209, 0.755201, 0.753187, 0.751165, 0.749136, 0.747101, 0.745058, 0.743008, 0.740951, 0.738887, 0.736817, 0.734739, 0.732654, 0.730563, 0.728464, 0.726359, 0.724247, 0.722128, 0.720003, 0.71787, 0.715731, 0.713585, 0.711432, 0.709273, 0.707107, 0.704934, 0.702755, 0.700569, 0.698376, 0.696177, 0.693971, 0.691759, 0.689541, 0.687315, 0.685084, 0.682846, 0.680601, 0.67835, 0.676093, 0.673829, 0.671559, 0.669283, 0.667, 0.664711, 0.662416, 0.660114, 0.657807, 0.655493, 0.653173, 0.650847, 0.648514, 0.646176, 0.643832, 0.641481, 0.639124, 0.636762, 0.634393, 0.632019, 0.629638, 0.627252, 0.624859, 0.622461, 0.620057, 0.617647, 0.615232, 0.61281, 0.610383, 0.60795, 0.605511, 0.603067, 0.600616, 0.598161, 0.595699, 0.593232, 0.59076, 0.588282, 0.585798, 0.583309, 0.580814, 0.578314, 0.575808, 0.573297, 0.570781, 0.568259, 0.565732, 0.563199, 0.560662, 0.558119, 0.55557, 0.553017, 0.550458, 0.547894, 0.545325, 0.542751, 0.540171, 0.537587, 0.534998, 0.532403, 0.529804, 0.527199, 0.52459, 0.521975, 0.519356, 0.516732, 0.514103, 0.511469, 0.50883, 0.506187, 0.503538, 0.500885, 0.498228, 0.495565, 0.492898, 0.490226, 0.48755, 0.484869, 0.482184, 0.479494, 0.476799, 0.4741, 0.471397, 0.468689, 0.465976, 0.46326, 0.460539, 0.457813, 0.455084, 0.45235, 0.449611, 0.446869, 0.444122, 0.441371, 0.438616, 0.435857, 0.433094, 0.430326, 0.427555, 0.42478, 0.422, 0.419217, 0.41643, 0.413638, 0.410843, 0.408044, 0.405241, 0.402435, 0.399624, 0.39681, 0.393992, 0.39117, 0.388345, 0.385516, 0.382683, 0.379847, 0.377007, 0.374164, 0.371317, 0.368467, 0.365613, 0.362756, 0.359895, 0.357031, 0.354164, 0.351293, 0.348419, 0.345541, 0.342661, 0.339777, 0.33689, 0.334, 0.331106, 0.32821, 0.32531, 0.322408, 0.319502, 0.316593, 0.313682, 0.310767, 0.30785, 0.304929, 0.302006, 0.29908, 0.296151, 0.293219, 0.290285, 0.287347, 0.284408, 0.281465, 0.27852, 0.275572, 0.272621, 0.269668, 0.266713, 0.263755, 0.260794, 0.257831, 0.254866, 0.251898, 0.248928, 0.245955, 0.24298, 0.240003, 0.237024, 0.234042, 0.231058, 0.228072, 0.225084, 0.222094, 0.219101, 0.216107, 0.21311, 0.210112, 0.207111, 0.204109, 0.201105, 0.198098, 0.19509, 0.19208, 0.189069, 0.186055, 0.18304, 0.180023, 0.177004, 0.173984, 0.170962, 0.167938, 0.164913, 0.161886, 0.158858, 0.155828, 0.152797, 0.149765, 0.14673, 0.143695, 0.140658, 0.13762, 0.134581, 0.13154, 0.128498, 0.125455, 0.122411, 0.119365, 0.116319, 0.113271, 0.110222, 0.107172, 0.104122, 0.10107, 0.0980171, 0.0949635, 0.091909, 0.0888536, 0.0857973, 0.0827403, 0.0796824, 0.0766239, 0.0735646, 0.0705046, 0.0674439, 0.0643826, 0.0613207, 0.0582583, 0.0551952, 0.0521317, 0.0490677, 0.0460032, 0.0429383, 0.0398729, 0.0368072, 0.0337412, 0.0306748, 0.0276081, 0.0245412, 0.0214741, 0.0184067, 0.0153392, 0.0122715, 0.00920375, 0.00613588, 0.00306796}; + +double allTheS [680] = {0, 0.19509, 0.382683, 0.55557, 0.707107, 0.83147, 0.92388, 0.980785, 0, 0.0490677, 0.0980171, 0.14673, 0.19509, 0.24298, 0.290285, 0.33689, 0.382683, 0.427555, 0.471397, 0.514103, 0.55557, 0.595699, 0.634393, 0.671559, 0.707107, 0.740951, 0.77301, 0.803208, 0.83147, 0.857729, 0.881921, 0.903989, 0.92388, 0.941544, 0.95694, 0.970031, 0.980785, 0.989177, 0.995185, 0.998795, 0, 0.0122715, 0.0245412, 0.0368072, 0.0490677, 0.0613207, 0.0735646, 0.0857973, 0.0980171, 0.110222, 0.122411, 0.134581, 0.14673, 0.158858, 0.170962, 0.18304, 0.19509, 0.207111, 0.219101, 0.231058, 0.24298, 0.254866, 0.266713, 0.27852, 0.290285, 0.302006, 0.313682, 0.32531, 0.33689, 0.348419, 0.359895, 0.371317, 0.382683, 0.393992, 0.405241, 0.41643, 0.427555, 0.438616, 0.449611, 0.460539, 0.471397, 0.482184, 0.492898, 0.503538, 0.514103, 0.52459, 0.534998, 0.545325, 0.55557, 0.565732, 0.575808, 0.585798, 0.595699, 0.605511, 0.615232, 0.624859, 0.634393, 0.643832, 0.653173, 0.662416, 0.671559, 0.680601, 0.689541, 0.698376, 0.707107, 0.715731, 0.724247, 0.732654, 0.740951, 0.749136, 0.757209, 0.765167, 0.77301, 0.780737, 0.788346, 0.795837, 0.803208, 0.810457, 0.817585, 0.824589, 0.83147, 0.838225, 0.844854, 0.851355, 0.857729, 0.863973, 0.870087, 0.87607, 0.881921, 0.88764, 0.893224, 0.898674, 0.903989, 0.909168, 0.91421, 0.919114, 0.92388, 0.928506, 0.932993, 0.937339, 0.941544, 0.945607, 0.949528, 0.953306, 0.95694, 0.960431, 0.963776, 0.966976, 0.970031, 0.97294, 0.975702, 0.978317, 0.980785, 0.983105, 0.985278, 0.987301, 0.989177, 0.990903, 0.99248, 0.993907, 0.995185, 0.996313, 0.99729, 0.998118, 0.998795, 0.999322, 0.999699, 0.999925, 0, 0.00306796, 0.00613588, 0.00920375, 0.0122715, 0.0153392, 0.0184067, 0.0214741, 0.0245412, 0.0276081, 0.0306748, 0.0337412, 0.0368072, 0.0398729, 0.0429383, 0.0460032, 0.0490677, 0.0521317, 0.0551952, 0.0582583, 0.0613207, 0.0643826, 0.0674439, 0.0705046, 0.0735646, 0.0766239, 0.0796824, 0.0827403, 0.0857973, 0.0888536, 0.091909, 0.0949635, 0.0980171, 0.10107, 0.104122, 0.107172, 0.110222, 0.113271, 0.116319, 0.119365, 0.122411, 0.125455, 0.128498, 0.13154, 0.134581, 0.13762, 0.140658, 0.143695, 0.14673, 0.149765, 0.152797, 0.155828, 0.158858, 0.161886, 0.164913, 0.167938, 0.170962, 0.173984, 0.177004, 0.180023, 0.18304, 0.186055, 0.189069, 0.19208, 0.19509, 0.198098, 0.201105, 0.204109, 0.207111, 0.210112, 0.21311, 0.216107, 0.219101, 0.222094, 0.225084, 0.228072, 0.231058, 0.234042, 0.237024, 0.240003, 0.24298, 0.245955, 0.248928, 0.251898, 0.254866, 0.257831, 0.260794, 0.263755, 0.266713, 0.269668, 0.272621, 0.275572, 0.27852, 0.281465, 0.284408, 0.287347, 0.290285, 0.293219, 0.296151, 0.29908, 0.302006, 0.304929, 0.30785, 0.310767, 0.313682, 0.316593, 0.319502, 0.322408, 0.32531, 0.32821, 0.331106, 0.334, 0.33689, 0.339777, 0.342661, 0.345541, 0.348419, 0.351293, 0.354164, 0.357031, 0.359895, 0.362756, 0.365613, 0.368467, 0.371317, 0.374164, 0.377007, 0.379847, 0.382683, 0.385516, 0.388345, 0.39117, 0.393992, 0.39681, 0.399624, 0.402435, 0.405241, 0.408044, 0.410843, 0.413638, 0.41643, 0.419217, 0.422, 0.42478, 0.427555, 0.430326, 0.433094, 0.435857, 0.438616, 0.441371, 0.444122, 0.446869, 0.449611, 0.45235, 0.455084, 0.457813, 0.460539, 0.46326, 0.465976, 0.468689, 0.471397, 0.4741, 0.476799, 0.479494, 0.482184, 0.484869, 0.48755, 0.490226, 0.492898, 0.495565, 0.498228, 0.500885, 0.503538, 0.506187, 0.50883, 0.511469, 0.514103, 0.516732, 0.519356, 0.521975, 0.52459, 0.527199, 0.529804, 0.532403, 0.534998, 0.537587, 0.540171, 0.542751, 0.545325, 0.547894, 0.550458, 0.553017, 0.55557, 0.558119, 0.560662, 0.563199, 0.565732, 0.568259, 0.570781, 0.573297, 0.575808, 0.578314, 0.580814, 0.583309, 0.585798, 0.588282, 0.59076, 0.593232, 0.595699, 0.598161, 0.600616, 0.603067, 0.605511, 0.60795, 0.610383, 0.61281, 0.615232, 0.617647, 0.620057, 0.622461, 0.624859, 0.627252, 0.629638, 0.632019, 0.634393, 0.636762, 0.639124, 0.641481, 0.643832, 0.646176, 0.648514, 0.650847, 0.653173, 0.655493, 0.657807, 0.660114, 0.662416, 0.664711, 0.667, 0.669283, 0.671559, 0.673829, 0.676093, 0.67835, 0.680601, 0.682846, 0.685084, 0.687315, 0.689541, 0.691759, 0.693971, 0.696177, 0.698376, 0.700569, 0.702755, 0.704934, 0.707107, 0.709273, 0.711432, 0.713585, 0.715731, 0.71787, 0.720003, 0.722128, 0.724247, 0.726359, 0.728464, 0.730563, 0.732654, 0.734739, 0.736817, 0.738887, 0.740951, 0.743008, 0.745058, 0.747101, 0.749136, 0.751165, 0.753187, 0.755201, 0.757209, 0.759209, 0.761202, 0.763188, 0.765167, 0.767139, 0.769103, 0.771061, 0.77301, 0.774953, 0.776888, 0.778817, 0.780737, 0.782651, 0.784557, 0.786455, 0.788346, 0.79023, 0.792107, 0.793975, 0.795837, 0.797691, 0.799537, 0.801376, 0.803208, 0.805031, 0.806848, 0.808656, 0.810457, 0.812251, 0.814036, 0.815814, 0.817585, 0.819348, 0.821103, 0.82285, 0.824589, 0.826321, 0.828045, 0.829761, 0.83147, 0.83317, 0.834863, 0.836548, 0.838225, 0.839894, 0.841555, 0.843208, 0.844854, 0.846491, 0.84812, 0.849742, 0.851355, 0.852961, 0.854558, 0.856147, 0.857729, 0.859302, 0.860867, 0.862424, 0.863973, 0.865514, 0.867046, 0.868571, 0.870087, 0.871595, 0.873095, 0.874587, 0.87607, 0.877545, 0.879012, 0.880471, 0.881921, 0.883363, 0.884797, 0.886223, 0.88764, 0.889048, 0.890449, 0.891841, 0.893224, 0.894599, 0.895966, 0.897325, 0.898674, 0.900016, 0.901349, 0.902673, 0.903989, 0.905297, 0.906596, 0.907886, 0.909168, 0.910441, 0.911706, 0.912962, 0.91421, 0.915449, 0.916679, 0.917901, 0.919114, 0.920318, 0.921514, 0.922701, 0.92388, 0.925049, 0.92621, 0.927363, 0.928506, 0.929641, 0.930767, 0.931884, 0.932993, 0.934093, 0.935184, 0.936266, 0.937339, 0.938404, 0.939459, 0.940506, 0.941544, 0.942573, 0.943593, 0.944605, 0.945607, 0.946601, 0.947586, 0.948561, 0.949528, 0.950486, 0.951435, 0.952375, 0.953306, 0.954228, 0.955141, 0.956045, 0.95694, 0.957826, 0.958703, 0.959572, 0.960431, 0.96128, 0.962121, 0.962953, 0.963776, 0.96459, 0.965394, 0.96619, 0.966976, 0.967754, 0.968522, 0.969281, 0.970031, 0.970772, 0.971504, 0.972226, 0.97294, 0.973644, 0.974339, 0.975025, 0.975702, 0.97637, 0.977028, 0.977677, 0.978317, 0.978948, 0.97957, 0.980182, 0.980785, 0.981379, 0.981964, 0.982539, 0.983105, 0.983662, 0.98421, 0.984749, 0.985278, 0.985798, 0.986308, 0.986809, 0.987301, 0.987784, 0.988258, 0.988722, 0.989177, 0.989622, 0.990058, 0.990485, 0.990903, 0.991311, 0.99171, 0.992099, 0.99248, 0.99285, 0.993212, 0.993564, 0.993907, 0.99424, 0.994565, 0.994879, 0.995185, 0.995481, 0.995767, 0.996045, 0.996313, 0.996571, 0.99682, 0.99706, 0.99729, 0.997511, 0.997723, 0.997925, 0.998118, 0.998302, 0.998476, 0.99864, 0.998795, 0.998941, 0.999078, 0.999205, 0.999322, 0.999431, 0.999529, 0.999619, 0.999699, 0.999769, 0.999831, 0.999882, 0.999925, 0.999958, 0.999981, 0.999995}; + +double allTheC2 [680] = {1, 0.92388, 0.707107, 0.382683, 2.22045e-16, -0.382683, -0.707107, -0.92388, 1, 0.995185, 0.980785, 0.95694, 0.92388, 0.881921, 0.83147, 0.77301, 0.707107, 0.634393, 0.55557, 0.471397, 0.382683, 0.290285, 0.19509, 0.0980171, -3.33067e-16, -0.0980171, -0.19509, -0.290285, -0.382683, -0.471397, -0.55557, -0.634393, -0.707107, -0.77301, -0.83147, -0.881921, -0.92388, -0.95694, -0.980785, -0.995185, 1, 0.999699, 0.998795, 0.99729, 0.995185, 0.99248, 0.989177, 0.985278, 0.980785, 0.975702, 0.970031, 0.963776, 0.95694, 0.949528, 0.941544, 0.932993, 0.92388, 0.91421, 0.903989, 0.893224, 0.881921, 0.870087, 0.857729, 0.844854, 0.83147, 0.817585, 0.803208, 0.788346, 0.77301, 0.757209, 0.740951, 0.724247, 0.707107, 0.689541, 0.671559, 0.653173, 0.634393, 0.615232, 0.595699, 0.575808, 0.55557, 0.534998, 0.514103, 0.492898, 0.471397, 0.449611, 0.427555, 0.405241, 0.382683, 0.359895, 0.33689, 0.313682, 0.290285, 0.266713, 0.24298, 0.219101, 0.19509, 0.170962, 0.14673, 0.122411, 0.0980171, 0.0735646, 0.0490677, 0.0245412, 1.11022e-15, -0.0245412, -0.0490677, -0.0735646, -0.0980171, -0.122411, -0.14673, -0.170962, -0.19509, -0.219101, -0.24298, -0.266713, -0.290285, -0.313682, -0.33689, -0.359895, -0.382683, -0.405241, -0.427555, -0.449611, -0.471397, -0.492898, -0.514103, -0.534998, -0.55557, -0.575808, -0.595699, -0.615232, -0.634393, -0.653173, -0.671559, -0.689541, -0.707107, -0.724247, -0.740951, -0.757209, -0.77301, -0.788346, -0.803208, -0.817585, -0.83147, -0.844854, -0.857729, -0.870087, -0.881921, -0.893224, -0.903989, -0.91421, -0.92388, -0.932993, -0.941544, -0.949528, -0.95694, -0.963776, -0.970031, -0.975702, -0.980785, -0.985278, -0.989177, -0.99248, -0.995185, -0.99729, -0.998795, -0.999699, 1, 0.999981, 0.999925, 0.999831, 0.999699, 0.999529, 0.999322, 0.999078, 0.998795, 0.998476, 0.998118, 0.997723, 0.99729, 0.99682, 0.996313, 0.995767, 0.995185, 0.994565, 0.993907, 0.993212, 0.99248, 0.99171, 0.990903, 0.990058, 0.989177, 0.988258, 0.987301, 0.986308, 0.985278, 0.98421, 0.983105, 0.981964, 0.980785, 0.97957, 0.978317, 0.977028, 0.975702, 0.974339, 0.97294, 0.971504, 0.970031, 0.968522, 0.966976, 0.965394, 0.963776, 0.962121, 0.960431, 0.958703, 0.95694, 0.955141, 0.953306, 0.951435, 0.949528, 0.947586, 0.945607, 0.943593, 0.941544, 0.939459, 0.937339, 0.935184, 0.932993, 0.930767, 0.928506, 0.92621, 0.92388, 0.921514, 0.919114, 0.916679, 0.91421, 0.911706, 0.909168, 0.906596, 0.903989, 0.901349, 0.898674, 0.895966, 0.893224, 0.890449, 0.88764, 0.884797, 0.881921, 0.879012, 0.87607, 0.873095, 0.870087, 0.867046, 0.863973, 0.860867, 0.857729, 0.854558, 0.851355, 0.84812, 0.844854, 0.841555, 0.838225, 0.834863, 0.83147, 0.828045, 0.824589, 0.821103, 0.817585, 0.814036, 0.810457, 0.806848, 0.803208, 0.799537, 0.795837, 0.792107, 0.788346, 0.784557, 0.780737, 0.776888, 0.77301, 0.769103, 0.765167, 0.761202, 0.757209, 0.753187, 0.749136, 0.745058, 0.740951, 0.736817, 0.732654, 0.728464, 0.724247, 0.720003, 0.715731, 0.711432, 0.707107, 0.702755, 0.698376, 0.693971, 0.689541, 0.685084, 0.680601, 0.676093, 0.671559, 0.667, 0.662416, 0.657807, 0.653173, 0.648514, 0.643832, 0.639124, 0.634393, 0.629638, 0.624859, 0.620057, 0.615232, 0.610383, 0.605511, 0.600616, 0.595699, 0.59076, 0.585798, 0.580814, 0.575808, 0.570781, 0.565732, 0.560662, 0.55557, 0.550458, 0.545325, 0.540171, 0.534998, 0.529804, 0.52459, 0.519356, 0.514103, 0.50883, 0.503538, 0.498228, 0.492898, 0.48755, 0.482184, 0.476799, 0.471397, 0.465976, 0.460539, 0.455084, 0.449611, 0.444122, 0.438616, 0.433094, 0.427555, 0.422, 0.41643, 0.410843, 0.405241, 0.399624, 0.393992, 0.388345, 0.382683, 0.377007, 0.371317, 0.365613, 0.359895, 0.354164, 0.348419, 0.342661, 0.33689, 0.331106, 0.32531, 0.319502, 0.313682, 0.30785, 0.302006, 0.296151, 0.290285, 0.284408, 0.27852, 0.272621, 0.266713, 0.260794, 0.254866, 0.248928, 0.24298, 0.237024, 0.231058, 0.225084, 0.219101, 0.21311, 0.207111, 0.201105, 0.19509, 0.189069, 0.18304, 0.177004, 0.170962, 0.164913, 0.158858, 0.152797, 0.14673, 0.140658, 0.134581, 0.128498, 0.122411, 0.116319, 0.110222, 0.104122, 0.0980171, 0.091909, 0.0857973, 0.0796824, 0.0735646, 0.0674439, 0.0613207, 0.0551952, 0.0490677, 0.0429383, 0.0368072, 0.0306748, 0.0245412, 0.0184067, 0.0122715, 0.00613588, 4.21885e-15, -0.00613588, -0.0122715, -0.0184067, -0.0245412, -0.0306748, -0.0368072, -0.0429383, -0.0490677, -0.0551952, -0.0613207, -0.0674439, -0.0735646, -0.0796824, -0.0857973, -0.091909, -0.0980171, -0.104122, -0.110222, -0.116319, -0.122411, -0.128498, -0.134581, -0.140658, -0.14673, -0.152797, -0.158858, -0.164913, -0.170962, -0.177004, -0.18304, -0.189069, -0.19509, -0.201105, -0.207111, -0.21311, -0.219101, -0.225084, -0.231058, -0.237024, -0.24298, -0.248928, -0.254866, -0.260794, -0.266713, -0.272621, -0.27852, -0.284408, -0.290285, -0.296151, -0.302006, -0.30785, -0.313682, -0.319502, -0.32531, -0.331106, -0.33689, -0.342661, -0.348419, -0.354164, -0.359895, -0.365613, -0.371317, -0.377007, -0.382683, -0.388345, -0.393992, -0.399624, -0.405241, -0.410843, -0.41643, -0.422, -0.427555, -0.433094, -0.438616, -0.444122, -0.449611, -0.455084, -0.460539, -0.465976, -0.471397, -0.476799, -0.482184, -0.48755, -0.492898, -0.498228, -0.503538, -0.50883, -0.514103, -0.519356, -0.52459, -0.529804, -0.534998, -0.540171, -0.545325, -0.550458, -0.55557, -0.560662, -0.565732, -0.570781, -0.575808, -0.580814, -0.585798, -0.59076, -0.595699, -0.600616, -0.605511, -0.610383, -0.615232, -0.620057, -0.624859, -0.629638, -0.634393, -0.639124, -0.643832, -0.648514, -0.653173, -0.657807, -0.662416, -0.667, -0.671559, -0.676093, -0.680601, -0.685084, -0.689541, -0.693971, -0.698376, -0.702755, -0.707107, -0.711432, -0.715731, -0.720003, -0.724247, -0.728464, -0.732654, -0.736817, -0.740951, -0.745058, -0.749136, -0.753187, -0.757209, -0.761202, -0.765167, -0.769103, -0.77301, -0.776888, -0.780737, -0.784557, -0.788346, -0.792107, -0.795837, -0.799537, -0.803208, -0.806848, -0.810457, -0.814036, -0.817585, -0.821103, -0.824589, -0.828045, -0.83147, -0.834863, -0.838225, -0.841555, -0.844854, -0.84812, -0.851355, -0.854558, -0.857729, -0.860867, -0.863973, -0.867046, -0.870087, -0.873095, -0.87607, -0.879012, -0.881921, -0.884797, -0.88764, -0.890449, -0.893224, -0.895966, -0.898674, -0.901349, -0.903989, -0.906596, -0.909168, -0.911706, -0.91421, -0.916679, -0.919114, -0.921514, -0.92388, -0.92621, -0.928506, -0.930767, -0.932993, -0.935184, -0.937339, -0.939459, -0.941544, -0.943593, -0.945607, -0.947586, -0.949528, -0.951435, -0.953306, -0.955141, -0.95694, -0.958703, -0.960431, -0.962121, -0.963776, -0.965394, -0.966976, -0.968522, -0.970031, -0.971504, -0.97294, -0.974339, -0.975702, -0.977028, -0.978317, -0.97957, -0.980785, -0.981964, -0.983105, -0.98421, -0.985278, -0.986308, -0.987301, -0.988258, -0.989177, -0.990058, -0.990903, -0.99171, -0.99248, -0.993212, -0.993907, -0.994565, -0.995185, -0.995767, -0.996313, -0.99682, -0.99729, -0.997723, -0.998118, -0.998476, -0.998795, -0.999078, -0.999322, -0.999529, -0.999699, -0.999831, -0.999925, -0.999981}; + +double allTheS2 [680] = {0, 0.382683, 0.707107, 0.92388, 1, 0.92388, 0.707107, 0.382683, 0, 0.0980171, 0.19509, 0.290285, 0.382683, 0.471397, 0.55557, 0.634393, 0.707107, 0.77301, 0.83147, 0.881921, 0.92388, 0.95694, 0.980785, 0.995185, 1, 0.995185, 0.980785, 0.95694, 0.92388, 0.881921, 0.83147, 0.77301, 0.707107, 0.634393, 0.55557, 0.471397, 0.382683, 0.290285, 0.19509, 0.0980171, 0, 0.0245412, 0.0490677, 0.0735646, 0.0980171, 0.122411, 0.14673, 0.170962, 0.19509, 0.219101, 0.24298, 0.266713, 0.290285, 0.313682, 0.33689, 0.359895, 0.382683, 0.405241, 0.427555, 0.449611, 0.471397, 0.492898, 0.514103, 0.534998, 0.55557, 0.575808, 0.595699, 0.615232, 0.634393, 0.653173, 0.671559, 0.689541, 0.707107, 0.724247, 0.740951, 0.757209, 0.77301, 0.788346, 0.803208, 0.817585, 0.83147, 0.844854, 0.857729, 0.870087, 0.881921, 0.893224, 0.903989, 0.91421, 0.92388, 0.932993, 0.941544, 0.949528, 0.95694, 0.963776, 0.970031, 0.975702, 0.980785, 0.985278, 0.989177, 0.99248, 0.995185, 0.99729, 0.998795, 0.999699, 1, 0.999699, 0.998795, 0.99729, 0.995185, 0.99248, 0.989177, 0.985278, 0.980785, 0.975702, 0.970031, 0.963776, 0.95694, 0.949528, 0.941544, 0.932993, 0.92388, 0.91421, 0.903989, 0.893224, 0.881921, 0.870087, 0.857729, 0.844854, 0.83147, 0.817585, 0.803208, 0.788346, 0.77301, 0.757209, 0.740951, 0.724247, 0.707107, 0.689541, 0.671559, 0.653173, 0.634393, 0.615232, 0.595699, 0.575808, 0.55557, 0.534998, 0.514103, 0.492898, 0.471397, 0.449611, 0.427555, 0.405241, 0.382683, 0.359895, 0.33689, 0.313682, 0.290285, 0.266713, 0.24298, 0.219101, 0.19509, 0.170962, 0.14673, 0.122411, 0.0980171, 0.0735646, 0.0490677, 0.0245412, 0, 0.00613588, 0.0122715, 0.0184067, 0.0245412, 0.0306748, 0.0368072, 0.0429383, 0.0490677, 0.0551952, 0.0613207, 0.0674439, 0.0735646, 0.0796824, 0.0857973, 0.091909, 0.0980171, 0.104122, 0.110222, 0.116319, 0.122411, 0.128498, 0.134581, 0.140658, 0.14673, 0.152797, 0.158858, 0.164913, 0.170962, 0.177004, 0.18304, 0.189069, 0.19509, 0.201105, 0.207111, 0.21311, 0.219101, 0.225084, 0.231058, 0.237024, 0.24298, 0.248928, 0.254866, 0.260794, 0.266713, 0.272621, 0.27852, 0.284408, 0.290285, 0.296151, 0.302006, 0.30785, 0.313682, 0.319502, 0.32531, 0.331106, 0.33689, 0.342661, 0.348419, 0.354164, 0.359895, 0.365613, 0.371317, 0.377007, 0.382683, 0.388345, 0.393992, 0.399624, 0.405241, 0.410843, 0.41643, 0.422, 0.427555, 0.433094, 0.438616, 0.444122, 0.449611, 0.455084, 0.460539, 0.465976, 0.471397, 0.476799, 0.482184, 0.48755, 0.492898, 0.498228, 0.503538, 0.50883, 0.514103, 0.519356, 0.52459, 0.529804, 0.534998, 0.540171, 0.545325, 0.550458, 0.55557, 0.560662, 0.565732, 0.570781, 0.575808, 0.580814, 0.585798, 0.59076, 0.595699, 0.600616, 0.605511, 0.610383, 0.615232, 0.620057, 0.624859, 0.629638, 0.634393, 0.639124, 0.643832, 0.648514, 0.653173, 0.657807, 0.662416, 0.667, 0.671559, 0.676093, 0.680601, 0.685084, 0.689541, 0.693971, 0.698376, 0.702755, 0.707107, 0.711432, 0.715731, 0.720003, 0.724247, 0.728464, 0.732654, 0.736817, 0.740951, 0.745058, 0.749136, 0.753187, 0.757209, 0.761202, 0.765167, 0.769103, 0.77301, 0.776888, 0.780737, 0.784557, 0.788346, 0.792107, 0.795837, 0.799537, 0.803208, 0.806848, 0.810457, 0.814036, 0.817585, 0.821103, 0.824589, 0.828045, 0.83147, 0.834863, 0.838225, 0.841555, 0.844854, 0.84812, 0.851355, 0.854558, 0.857729, 0.860867, 0.863973, 0.867046, 0.870087, 0.873095, 0.87607, 0.879012, 0.881921, 0.884797, 0.88764, 0.890449, 0.893224, 0.895966, 0.898674, 0.901349, 0.903989, 0.906596, 0.909168, 0.911706, 0.91421, 0.916679, 0.919114, 0.921514, 0.92388, 0.92621, 0.928506, 0.930767, 0.932993, 0.935184, 0.937339, 0.939459, 0.941544, 0.943593, 0.945607, 0.947586, 0.949528, 0.951435, 0.953306, 0.955141, 0.95694, 0.958703, 0.960431, 0.962121, 0.963776, 0.965394, 0.966976, 0.968522, 0.970031, 0.971504, 0.97294, 0.974339, 0.975702, 0.977028, 0.978317, 0.97957, 0.980785, 0.981964, 0.983105, 0.98421, 0.985278, 0.986308, 0.987301, 0.988258, 0.989177, 0.990058, 0.990903, 0.99171, 0.99248, 0.993212, 0.993907, 0.994565, 0.995185, 0.995767, 0.996313, 0.99682, 0.99729, 0.997723, 0.998118, 0.998476, 0.998795, 0.999078, 0.999322, 0.999529, 0.999699, 0.999831, 0.999925, 0.999981, 1, 0.999981, 0.999925, 0.999831, 0.999699, 0.999529, 0.999322, 0.999078, 0.998795, 0.998476, 0.998118, 0.997723, 0.99729, 0.99682, 0.996313, 0.995767, 0.995185, 0.994565, 0.993907, 0.993212, 0.99248, 0.99171, 0.990903, 0.990058, 0.989177, 0.988258, 0.987301, 0.986308, 0.985278, 0.98421, 0.983105, 0.981964, 0.980785, 0.97957, 0.978317, 0.977028, 0.975702, 0.974339, 0.97294, 0.971504, 0.970031, 0.968522, 0.966976, 0.965394, 0.963776, 0.962121, 0.960431, 0.958703, 0.95694, 0.955141, 0.953306, 0.951435, 0.949528, 0.947586, 0.945607, 0.943593, 0.941544, 0.939459, 0.937339, 0.935184, 0.932993, 0.930767, 0.928506, 0.92621, 0.92388, 0.921514, 0.919114, 0.916679, 0.91421, 0.911706, 0.909168, 0.906596, 0.903989, 0.901349, 0.898674, 0.895966, 0.893224, 0.890449, 0.88764, 0.884797, 0.881921, 0.879012, 0.87607, 0.873095, 0.870087, 0.867046, 0.863973, 0.860867, 0.857729, 0.854558, 0.851355, 0.84812, 0.844854, 0.841555, 0.838225, 0.834863, 0.83147, 0.828045, 0.824589, 0.821103, 0.817585, 0.814036, 0.810457, 0.806848, 0.803208, 0.799537, 0.795837, 0.792107, 0.788346, 0.784557, 0.780737, 0.776888, 0.77301, 0.769103, 0.765167, 0.761202, 0.757209, 0.753187, 0.749136, 0.745058, 0.740951, 0.736817, 0.732654, 0.728464, 0.724247, 0.720003, 0.715731, 0.711432, 0.707107, 0.702755, 0.698376, 0.693971, 0.689541, 0.685084, 0.680601, 0.676093, 0.671559, 0.667, 0.662416, 0.657807, 0.653173, 0.648514, 0.643832, 0.639124, 0.634393, 0.629638, 0.624859, 0.620057, 0.615232, 0.610383, 0.605511, 0.600616, 0.595699, 0.59076, 0.585798, 0.580814, 0.575808, 0.570781, 0.565732, 0.560662, 0.55557, 0.550458, 0.545325, 0.540171, 0.534998, 0.529804, 0.52459, 0.519356, 0.514103, 0.50883, 0.503538, 0.498228, 0.492898, 0.48755, 0.482184, 0.476799, 0.471397, 0.465976, 0.460539, 0.455084, 0.449611, 0.444122, 0.438616, 0.433094, 0.427555, 0.422, 0.41643, 0.410843, 0.405241, 0.399624, 0.393992, 0.388345, 0.382683, 0.377007, 0.371317, 0.365613, 0.359895, 0.354164, 0.348419, 0.342661, 0.33689, 0.331106, 0.32531, 0.319502, 0.313682, 0.30785, 0.302006, 0.296151, 0.290285, 0.284408, 0.27852, 0.272621, 0.266713, 0.260794, 0.254866, 0.248928, 0.24298, 0.237024, 0.231058, 0.225084, 0.219101, 0.21311, 0.207111, 0.201105, 0.19509, 0.189069, 0.18304, 0.177004, 0.170962, 0.164913, 0.158858, 0.152797, 0.14673, 0.140658, 0.134581, 0.128498, 0.122411, 0.116319, 0.110222, 0.104122, 0.0980171, 0.091909, 0.0857973, 0.0796824, 0.0735646, 0.0674439, 0.0613207, 0.0551952, 0.0490677, 0.0429383, 0.0368072, 0.0306748, 0.0245412, 0.0184067, 0.0122715, 0.00613588}; + +double allTheC3 [680] = {1, 0.83147, 0.382683, -0.19509, -0.707107, -0.980785, -0.92388, -0.55557, 1, 0.989177, 0.95694, 0.903989, 0.83147, 0.740951, 0.634393, 0.514103, 0.382683, 0.24298, 0.0980171, -0.0490677, -0.19509, -0.33689, -0.471397, -0.595699, -0.707107, -0.803208, -0.881921, -0.941544, -0.980785, -0.998795, -0.995185, -0.970031, -0.92388, -0.857729, -0.77301, -0.671559, -0.55557, -0.427555, -0.290285, -0.14673, 1, 0.999322, 0.99729, 0.993907, 0.989177, 0.983105, 0.975702, 0.966976, 0.95694, 0.945607, 0.932993, 0.919114, 0.903989, 0.88764, 0.870087, 0.851355, 0.83147, 0.810457, 0.788346, 0.765167, 0.740951, 0.715731, 0.689541, 0.662416, 0.634393, 0.605511, 0.575808, 0.545325, 0.514103, 0.482184, 0.449611, 0.41643, 0.382683, 0.348419, 0.313682, 0.27852, 0.24298, 0.207111, 0.170962, 0.134581, 0.0980171, 0.0613207, 0.0245412, -0.0122715, -0.0490677, -0.0857973, -0.122411, -0.158858, -0.19509, -0.231058, -0.266713, -0.302006, -0.33689, -0.371317, -0.405241, -0.438616, -0.471397, -0.503538, -0.534998, -0.565732, -0.595699, -0.624859, -0.653173, -0.680601, -0.707107, -0.732654, -0.757209, -0.780737, -0.803208, -0.824589, -0.844854, -0.863973, -0.881921, -0.898674, -0.91421, -0.928506, -0.941544, -0.953306, -0.963776, -0.97294, -0.980785, -0.987301, -0.99248, -0.996313, -0.998795, -0.999925, -0.999699, -0.998118, -0.995185, -0.990903, -0.985278, -0.978317, -0.970031, -0.960431, -0.949528, -0.937339, -0.92388, -0.909168, -0.893224, -0.87607, -0.857729, -0.838225, -0.817585, -0.795837, -0.77301, -0.749136, -0.724247, -0.698376, -0.671559, -0.643832, -0.615232, -0.585798, -0.55557, -0.52459, -0.492898, -0.460539, -0.427555, -0.393992, -0.359895, -0.32531, -0.290285, -0.254866, -0.219101, -0.18304, -0.14673, -0.110222, -0.0735646, -0.0368072, 1, 0.999958, 0.999831, 0.999619, 0.999322, 0.998941, 0.998476, 0.997925, 0.99729, 0.996571, 0.995767, 0.994879, 0.993907, 0.99285, 0.99171, 0.990485, 0.989177, 0.987784, 0.986308, 0.984749, 0.983105, 0.981379, 0.97957, 0.977677, 0.975702, 0.973644, 0.971504, 0.969281, 0.966976, 0.96459, 0.962121, 0.959572, 0.95694, 0.954228, 0.951435, 0.948561, 0.945607, 0.942573, 0.939459, 0.936266, 0.932993, 0.929641, 0.92621, 0.922701, 0.919114, 0.915449, 0.911706, 0.907886, 0.903989, 0.900016, 0.895966, 0.891841, 0.88764, 0.883363, 0.879012, 0.874587, 0.870087, 0.865514, 0.860867, 0.856147, 0.851355, 0.846491, 0.841555, 0.836548, 0.83147, 0.826321, 0.821103, 0.815814, 0.810457, 0.805031, 0.799537, 0.793975, 0.788346, 0.782651, 0.776888, 0.771061, 0.765167, 0.759209, 0.753187, 0.747101, 0.740951, 0.734739, 0.728464, 0.722128, 0.715731, 0.709273, 0.702755, 0.696177, 0.689541, 0.682846, 0.676093, 0.669283, 0.662416, 0.655493, 0.648514, 0.641481, 0.634393, 0.627252, 0.620057, 0.61281, 0.605511, 0.598161, 0.59076, 0.583309, 0.575808, 0.568259, 0.560662, 0.553017, 0.545325, 0.537587, 0.529804, 0.521975, 0.514103, 0.506187, 0.498228, 0.490226, 0.482184, 0.4741, 0.465976, 0.457813, 0.449611, 0.441371, 0.433094, 0.42478, 0.41643, 0.408044, 0.399624, 0.39117, 0.382683, 0.374164, 0.365613, 0.357031, 0.348419, 0.339777, 0.331106, 0.322408, 0.313682, 0.304929, 0.296151, 0.287347, 0.27852, 0.269668, 0.260794, 0.251898, 0.24298, 0.234042, 0.225084, 0.216107, 0.207111, 0.198098, 0.189069, 0.180023, 0.170962, 0.161886, 0.152797, 0.143695, 0.134581, 0.125455, 0.116319, 0.107172, 0.0980171, 0.0888536, 0.0796824, 0.0705046, 0.0613207, 0.0521317, 0.0429383, 0.0337412, 0.0245412, 0.0153392, 0.00613588, -0.00306796, -0.0122715, -0.0214741, -0.0306748, -0.0398729, -0.0490677, -0.0582583, -0.0674439, -0.0766239, -0.0857973, -0.0949635, -0.104122, -0.113271, -0.122411, -0.13154, -0.140658, -0.149765, -0.158858, -0.167938, -0.177004, -0.186055, -0.19509, -0.204109, -0.21311, -0.222094, -0.231058, -0.240003, -0.248928, -0.257831, -0.266713, -0.275572, -0.284408, -0.293219, -0.302006, -0.310767, -0.319502, -0.32821, -0.33689, -0.345541, -0.354164, -0.362756, -0.371317, -0.379847, -0.388345, -0.39681, -0.405241, -0.413638, -0.422, -0.430326, -0.438616, -0.446869, -0.455084, -0.46326, -0.471397, -0.479494, -0.48755, -0.495565, -0.503538, -0.511469, -0.519356, -0.527199, -0.534998, -0.542751, -0.550458, -0.558119, -0.565732, -0.573297, -0.580814, -0.588282, -0.595699, -0.603067, -0.610383, -0.617647, -0.624859, -0.632019, -0.639124, -0.646176, -0.653173, -0.660114, -0.667, -0.673829, -0.680601, -0.687315, -0.693971, -0.700569, -0.707107, -0.713585, -0.720003, -0.726359, -0.732654, -0.738887, -0.745058, -0.751165, -0.757209, -0.763188, -0.769103, -0.774953, -0.780737, -0.786455, -0.792107, -0.797691, -0.803208, -0.808656, -0.814036, -0.819348, -0.824589, -0.829761, -0.834863, -0.839894, -0.844854, -0.849742, -0.854558, -0.859302, -0.863973, -0.868571, -0.873095, -0.877545, -0.881921, -0.886223, -0.890449, -0.894599, -0.898674, -0.902673, -0.906596, -0.910441, -0.91421, -0.917901, -0.921514, -0.925049, -0.928506, -0.931884, -0.935184, -0.938404, -0.941544, -0.944605, -0.947586, -0.950486, -0.953306, -0.956045, -0.958703, -0.96128, -0.963776, -0.96619, -0.968522, -0.970772, -0.97294, -0.975025, -0.977028, -0.978948, -0.980785, -0.982539, -0.98421, -0.985798, -0.987301, -0.988722, -0.990058, -0.991311, -0.99248, -0.993564, -0.994565, -0.995481, -0.996313, -0.99706, -0.997723, -0.998302, -0.998795, -0.999205, -0.999529, -0.999769, -0.999925, -0.999995, -0.999981, -0.999882, -0.999699, -0.999431, -0.999078, -0.99864, -0.998118, -0.997511, -0.99682, -0.996045, -0.995185, -0.99424, -0.993212, -0.992099, -0.990903, -0.989622, -0.988258, -0.986809, -0.985278, -0.983662, -0.981964, -0.980182, -0.978317, -0.97637, -0.974339, -0.972226, -0.970031, -0.967754, -0.965394, -0.962953, -0.960431, -0.957826, -0.955141, -0.952375, -0.949528, -0.946601, -0.943593, -0.940506, -0.937339, -0.934093, -0.930767, -0.927363, -0.92388, -0.920318, -0.916679, -0.912962, -0.909168, -0.905297, -0.901349, -0.897325, -0.893224, -0.889048, -0.884797, -0.880471, -0.87607, -0.871595, -0.867046, -0.862424, -0.857729, -0.852961, -0.84812, -0.843208, -0.838225, -0.83317, -0.828045, -0.82285, -0.817585, -0.812251, -0.806848, -0.801376, -0.795837, -0.79023, -0.784557, -0.778817, -0.77301, -0.767139, -0.761202, -0.755201, -0.749136, -0.743008, -0.736817, -0.730563, -0.724247, -0.71787, -0.711432, -0.704934, -0.698376, -0.691759, -0.685084, -0.67835, -0.671559, -0.664711, -0.657807, -0.650847, -0.643832, -0.636762, -0.629638, -0.622461, -0.615232, -0.60795, -0.600616, -0.593232, -0.585798, -0.578314, -0.570781, -0.563199, -0.55557, -0.547894, -0.540171, -0.532403, -0.52459, -0.516732, -0.50883, -0.500885, -0.492898, -0.484869, -0.476799, -0.468689, -0.460539, -0.45235, -0.444122, -0.435857, -0.427555, -0.419217, -0.410843, -0.402435, -0.393992, -0.385516, -0.377007, -0.368467, -0.359895, -0.351293, -0.342661, -0.334, -0.32531, -0.316593, -0.30785, -0.29908, -0.290285, -0.281465, -0.272621, -0.263755, -0.254866, -0.245955, -0.237024, -0.228072, -0.219101, -0.210112, -0.201105, -0.19208, -0.18304, -0.173984, -0.164913, -0.155828, -0.14673, -0.13762, -0.128498, -0.119365, -0.110222, -0.10107, -0.091909, -0.0827403, -0.0735646, -0.0643826, -0.0551952, -0.0460032, -0.0368072, -0.0276081, -0.0184067, -0.00920375}; + +double allTheS3 [680] = {0, 0.55557, 0.92388, 0.980785, 0.707107, 0.19509, -0.382683, -0.83147, 0, 0.14673, 0.290285, 0.427555, 0.55557, 0.671559, 0.77301, 0.857729, 0.92388, 0.970031, 0.995185, 0.998795, 0.980785, 0.941544, 0.881921, 0.803208, 0.707107, 0.595699, 0.471397, 0.33689, 0.19509, 0.0490677, -0.0980171, -0.24298, -0.382683, -0.514103, -0.634393, -0.740951, -0.83147, -0.903989, -0.95694, -0.989177, 0, 0.0368072, 0.0735646, 0.110222, 0.14673, 0.18304, 0.219101, 0.254866, 0.290285, 0.32531, 0.359895, 0.393992, 0.427555, 0.460539, 0.492898, 0.52459, 0.55557, 0.585798, 0.615232, 0.643832, 0.671559, 0.698376, 0.724247, 0.749136, 0.77301, 0.795837, 0.817585, 0.838225, 0.857729, 0.87607, 0.893224, 0.909168, 0.92388, 0.937339, 0.949528, 0.960431, 0.970031, 0.978317, 0.985278, 0.990903, 0.995185, 0.998118, 0.999699, 0.999925, 0.998795, 0.996313, 0.99248, 0.987301, 0.980785, 0.97294, 0.963776, 0.953306, 0.941544, 0.928506, 0.91421, 0.898674, 0.881921, 0.863973, 0.844854, 0.824589, 0.803208, 0.780737, 0.757209, 0.732654, 0.707107, 0.680601, 0.653173, 0.624859, 0.595699, 0.565732, 0.534998, 0.503538, 0.471397, 0.438616, 0.405241, 0.371317, 0.33689, 0.302006, 0.266713, 0.231058, 0.19509, 0.158858, 0.122411, 0.0857973, 0.0490677, 0.0122715, -0.0245412, -0.0613207, -0.0980171, -0.134581, -0.170962, -0.207111, -0.24298, -0.27852, -0.313682, -0.348419, -0.382683, -0.41643, -0.449611, -0.482184, -0.514103, -0.545325, -0.575808, -0.605511, -0.634393, -0.662416, -0.689541, -0.715731, -0.740951, -0.765167, -0.788346, -0.810457, -0.83147, -0.851355, -0.870087, -0.88764, -0.903989, -0.919114, -0.932993, -0.945607, -0.95694, -0.966976, -0.975702, -0.983105, -0.989177, -0.993907, -0.99729, -0.999322, 0, 0.00920375, 0.0184067, 0.0276081, 0.0368072, 0.0460032, 0.0551952, 0.0643826, 0.0735646, 0.0827403, 0.091909, 0.10107, 0.110222, 0.119365, 0.128498, 0.13762, 0.14673, 0.155828, 0.164913, 0.173984, 0.18304, 0.19208, 0.201105, 0.210112, 0.219101, 0.228072, 0.237024, 0.245955, 0.254866, 0.263755, 0.272621, 0.281465, 0.290285, 0.29908, 0.30785, 0.316593, 0.32531, 0.334, 0.342661, 0.351293, 0.359895, 0.368467, 0.377007, 0.385516, 0.393992, 0.402435, 0.410843, 0.419217, 0.427555, 0.435857, 0.444122, 0.45235, 0.460539, 0.468689, 0.476799, 0.484869, 0.492898, 0.500885, 0.50883, 0.516732, 0.52459, 0.532403, 0.540171, 0.547894, 0.55557, 0.563199, 0.570781, 0.578314, 0.585798, 0.593232, 0.600616, 0.60795, 0.615232, 0.622461, 0.629638, 0.636762, 0.643832, 0.650847, 0.657807, 0.664711, 0.671559, 0.67835, 0.685084, 0.691759, 0.698376, 0.704934, 0.711432, 0.71787, 0.724247, 0.730563, 0.736817, 0.743008, 0.749136, 0.755201, 0.761202, 0.767139, 0.77301, 0.778817, 0.784557, 0.79023, 0.795837, 0.801376, 0.806848, 0.812251, 0.817585, 0.82285, 0.828045, 0.83317, 0.838225, 0.843208, 0.84812, 0.852961, 0.857729, 0.862424, 0.867046, 0.871595, 0.87607, 0.880471, 0.884797, 0.889048, 0.893224, 0.897325, 0.901349, 0.905297, 0.909168, 0.912962, 0.916679, 0.920318, 0.92388, 0.927363, 0.930767, 0.934093, 0.937339, 0.940506, 0.943593, 0.946601, 0.949528, 0.952375, 0.955141, 0.957826, 0.960431, 0.962953, 0.965394, 0.967754, 0.970031, 0.972226, 0.974339, 0.97637, 0.978317, 0.980182, 0.981964, 0.983662, 0.985278, 0.986809, 0.988258, 0.989622, 0.990903, 0.992099, 0.993212, 0.99424, 0.995185, 0.996045, 0.99682, 0.997511, 0.998118, 0.99864, 0.999078, 0.999431, 0.999699, 0.999882, 0.999981, 0.999995, 0.999925, 0.999769, 0.999529, 0.999205, 0.998795, 0.998302, 0.997723, 0.99706, 0.996313, 0.995481, 0.994565, 0.993564, 0.99248, 0.991311, 0.990058, 0.988722, 0.987301, 0.985798, 0.98421, 0.982539, 0.980785, 0.978948, 0.977028, 0.975025, 0.97294, 0.970772, 0.968522, 0.96619, 0.963776, 0.96128, 0.958703, 0.956045, 0.953306, 0.950486, 0.947586, 0.944605, 0.941544, 0.938404, 0.935184, 0.931884, 0.928506, 0.925049, 0.921514, 0.917901, 0.91421, 0.910441, 0.906596, 0.902673, 0.898674, 0.894599, 0.890449, 0.886223, 0.881921, 0.877545, 0.873095, 0.868571, 0.863973, 0.859302, 0.854558, 0.849742, 0.844854, 0.839894, 0.834863, 0.829761, 0.824589, 0.819348, 0.814036, 0.808656, 0.803208, 0.797691, 0.792107, 0.786455, 0.780737, 0.774953, 0.769103, 0.763188, 0.757209, 0.751165, 0.745058, 0.738887, 0.732654, 0.726359, 0.720003, 0.713585, 0.707107, 0.700569, 0.693971, 0.687315, 0.680601, 0.673829, 0.667, 0.660114, 0.653173, 0.646176, 0.639124, 0.632019, 0.624859, 0.617647, 0.610383, 0.603067, 0.595699, 0.588282, 0.580814, 0.573297, 0.565732, 0.558119, 0.550458, 0.542751, 0.534998, 0.527199, 0.519356, 0.511469, 0.503538, 0.495565, 0.48755, 0.479494, 0.471397, 0.46326, 0.455084, 0.446869, 0.438616, 0.430326, 0.422, 0.413638, 0.405241, 0.39681, 0.388345, 0.379847, 0.371317, 0.362756, 0.354164, 0.345541, 0.33689, 0.32821, 0.319502, 0.310767, 0.302006, 0.293219, 0.284408, 0.275572, 0.266713, 0.257831, 0.248928, 0.240003, 0.231058, 0.222094, 0.21311, 0.204109, 0.19509, 0.186055, 0.177004, 0.167938, 0.158858, 0.149765, 0.140658, 0.13154, 0.122411, 0.113271, 0.104122, 0.0949635, 0.0857973, 0.0766239, 0.0674439, 0.0582583, 0.0490677, 0.0398729, 0.0306748, 0.0214741, 0.0122715, 0.00306796, -0.00613588, -0.0153392, -0.0245412, -0.0337412, -0.0429383, -0.0521317, -0.0613207, -0.0705046, -0.0796824, -0.0888536, -0.0980171, -0.107172, -0.116319, -0.125455, -0.134581, -0.143695, -0.152797, -0.161886, -0.170962, -0.180023, -0.189069, -0.198098, -0.207111, -0.216107, -0.225084, -0.234042, -0.24298, -0.251898, -0.260794, -0.269668, -0.27852, -0.287347, -0.296151, -0.304929, -0.313682, -0.322408, -0.331106, -0.339777, -0.348419, -0.357031, -0.365613, -0.374164, -0.382683, -0.39117, -0.399624, -0.408044, -0.41643, -0.42478, -0.433094, -0.441371, -0.449611, -0.457813, -0.465976, -0.4741, -0.482184, -0.490226, -0.498228, -0.506187, -0.514103, -0.521975, -0.529804, -0.537587, -0.545325, -0.553017, -0.560662, -0.568259, -0.575808, -0.583309, -0.59076, -0.598161, -0.605511, -0.61281, -0.620057, -0.627252, -0.634393, -0.641481, -0.648514, -0.655493, -0.662416, -0.669283, -0.676093, -0.682846, -0.689541, -0.696177, -0.702755, -0.709273, -0.715731, -0.722128, -0.728464, -0.734739, -0.740951, -0.747101, -0.753187, -0.759209, -0.765167, -0.771061, -0.776888, -0.782651, -0.788346, -0.793975, -0.799537, -0.805031, -0.810457, -0.815814, -0.821103, -0.826321, -0.83147, -0.836548, -0.841555, -0.846491, -0.851355, -0.856147, -0.860867, -0.865514, -0.870087, -0.874587, -0.879012, -0.883363, -0.88764, -0.891841, -0.895966, -0.900016, -0.903989, -0.907886, -0.911706, -0.915449, -0.919114, -0.922701, -0.92621, -0.929641, -0.932993, -0.936266, -0.939459, -0.942573, -0.945607, -0.948561, -0.951435, -0.954228, -0.95694, -0.959572, -0.962121, -0.96459, -0.966976, -0.969281, -0.971504, -0.973644, -0.975702, -0.977677, -0.97957, -0.981379, -0.983105, -0.984749, -0.986308, -0.987784, -0.989177, -0.990485, -0.99171, -0.99285, -0.993907, -0.994879, -0.995767, -0.996571, -0.99729, -0.997925, -0.998476, -0.998941, -0.999322, -0.999619, -0.999831, -0.999958}; + +void +fft_dit4_core_p1(double *fr, double *fi, ulong ldn) +// Auxiliary routine for fft_dit4() +// Decimation in time (DIT) radix-4 FFT +// Input data must be in revbin_permuted order +// ldn := base-2 logarithm of the array length +// Fixed isign = +1 +{ + int tmp = 0;// for getting the precalculated cos and sin + + const ulong n = (1UL<<ldn); + + if ( n<=2 ) + { + if ( n==2 ) + { + sumdiff(fr[0], fr[1]); + sumdiff(fi[0], fi[1]); + } + return; + } + + + ulong ldm = (ldn&1); + if ( ldm!=0 ) // n is not a power of 4, need a radix-8 step + { + for (ulong i0=0; i0<n; i0+=8) fft8_dit_core_p1(fr+i0, fi+i0); + } + else + { + for (ulong i0=0; i0<n; i0+=4) + { + double xr, yr, ur, vr, xi, yi, ui, vi; + ulong i1 = i0 + 1; + ulong i2 = i1 + 1; + ulong i3 = i2 + 1; + + sumdiff(fr[i0], fr[i1], xr, ur); + sumdiff(fr[i2], fr[i3], yr, vi); + sumdiff(fi[i0], fi[i1], xi, ui); + sumdiff(fi[i3], fi[i2], yi, vr); + + sumdiff(ui, vi, fi[i1], fi[i3]); + sumdiff(xi, yi, fi[i0], fi[i2]); + sumdiff(ur, vr, fr[i1], fr[i3]); + sumdiff(xr, yr, fr[i0], fr[i2]); + } + } + ldm += 2*LX; + + + for ( ; ldm<=ldn; ldm+=LX) + { + ulong m = (1UL<<ldm); + ulong m4 = (m>>LX); + + for (ulong j=0; j<m4; j++) + { + double c, s, c2, s2, c3, s3; + c = allTheC[tmp]; + s = allTheS[tmp]; + c2 = allTheC2[tmp]; + s2 = allTheS2[tmp]; + c3 = allTheC3[tmp]; + s3 = allTheS3[tmp]; + tmp = tmp + 1; + + for (ulong r=0; r<n; r+=m) + { + ulong i0 = j + r; + ulong i1 = i0 + m4; + ulong i2 = i1 + m4; + ulong i3 = i2 + m4; + + double xr, yr, ur, vr, xi, yi, ui, vi; + cmult(c2, s2, fr[i1], fi[i1], xr, xi); + + sumdiff3_r(xr, fr[i0], ur); + sumdiff3_r(xi, fi[i0], ui); + + cmult(c, s, fr[i2], fi[i2], yr, vr); + cmult(c3, s3, fr[i3], fi[i3], vi, yi); + + sumdiff(yr, vi); + sumdiff(yi, vr); + + sumdiff(ur, vr, fr[i1], fr[i3]); + sumdiff(ui, vi, fi[i1], fi[i3]); + sumdiff(xr, yr, fr[i0], fr[i2]); + sumdiff(xi, yi, fi[i0], fi[i2]); + } + } + } +} +// ------------------------- +////////////////////////////////////////////////////////////////////////////////////////////////// + + +//// store all the precalculated indexes in the array +int PermutedIndex [2048] = {0, 1024, 512, 1536, 256, 1280, 768, 1792, 128, 1152, 640, 1664, 384, 1408, 896, 1920, 64, 1088, 576, 1600, 320, 1344, 832, 1856, 192, 1216, 704, 1728, 448, 1472, 960, 1984, 32, 1056, 544, 1568, 288, 1312, 800, 1824, 160, 1184, 672, 1696, 416, 1440, 928, 1952, 96, 1120, 608, 1632, 352, 1376, 864, 1888, 224, 1248, 736, 1760, 480, 1504, 992, 2016, 16, 1040, 528, 1552, 272, 1296, 784, 1808, 144, 1168, 656, 1680, 400, 1424, 912, 1936, 80, 1104, 592, 1616, 336, 1360, 848, 1872, 208, 1232, 720, 1744, 464, 1488, 976, 2000, 48, 1072, 560, 1584, 304, 1328, 816, 1840, 176, 1200, 688, 1712, 432, 1456, 944, 1968, 112, 1136, 624, 1648, 368, 1392, 880, 1904, 240, 1264, 752, 1776, 496, 1520, 1008, 2032, 8, 1032, 520, 1544, 264, 1288, 776, 1800, 136, 1160, 648, 1672, 392, 1416, 904, 1928, 72, 1096, 584, 1608, 328, 1352, 840, 1864, 200, 1224, 712, 1736, 456, 1480, 968, 1992, 40, 1064, 552, 1576, 296, 1320, 808, 1832, 168, 1192, 680, 1704, 424, 1448, 936, 1960, 104, 1128, 616, 1640, 360, 1384, 872, 1896, 232, 1256, 744, 1768, 488, 1512, 1000, 2024, 24, 1048, 536, 1560, 280, 1304, 792, 1816, 152, 1176, 664, 1688, 408, 1432, 920, 1944, 88, 1112, 600, 1624, 344, 1368, 856, 1880, 216, 1240, 728, 1752, 472, 1496, 984, 2008, 56, 1080, 568, 1592, 312, 1336, 824, 1848, 184, 1208, 696, 1720, 440, 1464, 952, 1976, 120, 1144, 632, 1656, 376, 1400, 888, 1912, 248, 1272, 760, 1784, 504, 1528, 1016, 2040, 4, 1028, 516, 1540, 260, 1284, 772, 1796, 132, 1156, 644, 1668, 388, 1412, 900, 1924, 68, 1092, 580, 1604, 324, 1348, 836, 1860, 196, 1220, 708, 1732, 452, 1476, 964, 1988, 36, 1060, 548, 1572, 292, 1316, 804, 1828, 164, 1188, 676, 1700, 420, 1444, 932, 1956, 100, 1124, 612, 1636, 356, 1380, 868, 1892, 228, 1252, 740, 1764, 484, 1508, 996, 2020, 20, 1044, 532, 1556, 276, 1300, 788, 1812, 148, 1172, 660, 1684, 404, 1428, 916, 1940, 84, 1108, 596, 1620, 340, 1364, 852, 1876, 212, 1236, 724, 1748, 468, 1492, 980, 2004, 52, 1076, 564, 1588, 308, 1332, 820, 1844, 180, 1204, 692, 1716, 436, 1460, 948, 1972, 116, 1140, 628, 1652, 372, 1396, 884, 1908, 244, 1268, 756, 1780, 500, 1524, 1012, 2036, 12, 1036, 524, 1548, 268, 1292, 780, 1804, 140, 1164, 652, 1676, 396, 1420, 908, 1932, 76, 1100, 588, 1612, 332, 1356, 844, 1868, 204, 1228, 716, 1740, 460, 1484, 972, 1996, 44, 1068, 556, 1580, 300, 1324, 812, 1836, 172, 1196, 684, 1708, 428, 1452, 940, 1964, 108, 1132, 620, 1644, 364, 1388, 876, 1900, 236, 1260, 748, 1772, 492, 1516, 1004, 2028, 28, 1052, 540, 1564, 284, 1308, 796, 1820, 156, 1180, 668, 1692, 412, 1436, 924, 1948, 92, 1116, 604, 1628, 348, 1372, 860, 1884, 220, 1244, 732, 1756, 476, 1500, 988, 2012, 60, 1084, 572, 1596, 316, 1340, 828, 1852, 188, 1212, 700, 1724, 444, 1468, 956, 1980, 124, 1148, 636, 1660, 380, 1404, 892, 1916, 252, 1276, 764, 1788, 508, 1532, 1020, 2044, 2, 1026, 514, 1538, 258, 1282, 770, 1794, 130, 1154, 642, 1666, 386, 1410, 898, 1922, 66, 1090, 578, 1602, 322, 1346, 834, 1858, 194, 1218, 706, 1730, 450, 1474, 962, 1986, 34, 1058, 546, 1570, 290, 1314, 802, 1826, 162, 1186, 674, 1698, 418, 1442, 930, 1954, 98, 1122, 610, 1634, 354, 1378, 866, 1890, 226, 1250, 738, 1762, 482, 1506, 994, 2018, 18, 1042, 530, 1554, 274, 1298, 786, 1810, 146, 1170, 658, 1682, 402, 1426, 914, 1938, 82, 1106, 594, 1618, 338, 1362, 850, 1874, 210, 1234, 722, 1746, 466, 1490, 978, 2002, 50, 1074, 562, 1586, 306, 1330, 818, 1842, 178, 1202, 690, 1714, 434, 1458, 946, 1970, 114, 1138, 626, 1650, 370, 1394, 882, 1906, 242, 1266, 754, 1778, 498, 1522, 1010, 2034, 10, 1034, 522, 1546, 266, 1290, 778, 1802, 138, 1162, 650, 1674, 394, 1418, 906, 1930, 74, 1098, 586, 1610, 330, 1354, 842, 1866, 202, 1226, 714, 1738, 458, 1482, 970, 1994, 42, 1066, 554, 1578, 298, 1322, 810, 1834, 170, 1194, 682, 1706, 426, 1450, 938, 1962, 106, 1130, 618, 1642, 362, 1386, 874, 1898, 234, 1258, 746, 1770, 490, 1514, 1002, 2026, 26, 1050, 538, 1562, 282, 1306, 794, 1818, 154, 1178, 666, 1690, 410, 1434, 922, 1946, 90, 1114, 602, 1626, 346, 1370, 858, 1882, 218, 1242, 730, 1754, 474, 1498, 986, 2010, 58, 1082, 570, 1594, 314, 1338, 826, 1850, 186, 1210, 698, 1722, 442, 1466, 954, 1978, 122, 1146, 634, 1658, 378, 1402, 890, 1914, 250, 1274, 762, 1786, 506, 1530, 1018, 2042, 6, 1030, 518, 1542, 262, 1286, 774, 1798, 134, 1158, 646, 1670, 390, 1414, 902, 1926, 70, 1094, 582, 1606, 326, 1350, 838, 1862, 198, 1222, 710, 1734, 454, 1478, 966, 1990, 38, 1062, 550, 1574, 294, 1318, 806, 1830, 166, 1190, 678, 1702, 422, 1446, 934, 1958, 102, 1126, 614, 1638, 358, 1382, 870, 1894, 230, 1254, 742, 1766, 486, 1510, 998, 2022, 22, 1046, 534, 1558, 278, 1302, 790, 1814, 150, 1174, 662, 1686, 406, 1430, 918, 1942, 86, 1110, 598, 1622, 342, 1366, 854, 1878, 214, 1238, 726, 1750, 470, 1494, 982, 2006, 54, 1078, 566, 1590, 310, 1334, 822, 1846, 182, 1206, 694, 1718, 438, 1462, 950, 1974, 118, 1142, 630, 1654, 374, 1398, 886, 1910, 246, 1270, 758, 1782, 502, 1526, 1014, 2038, 14, 1038, 526, 1550, 270, 1294, 782, 1806, 142, 1166, 654, 1678, 398, 1422, 910, 1934, 78, 1102, 590, 1614, 334, 1358, 846, 1870, 206, 1230, 718, 1742, 462, 1486, 974, 1998, 46, 1070, 558, 1582, 302, 1326, 814, 1838, 174, 1198, 686, 1710, 430, 1454, 942, 1966, 110, 1134, 622, 1646, 366, 1390, 878, 1902, 238, 1262, 750, 1774, 494, 1518, 1006, 2030, 30, 1054, 542, 1566, 286, 1310, 798, 1822, 158, 1182, 670, 1694, 414, 1438, 926, 1950, 94, 1118, 606, 1630, 350, 1374, 862, 1886, 222, 1246, 734, 1758, 478, 1502, 990, 2014, 62, 1086, 574, 1598, 318, 1342, 830, 1854, 190, 1214, 702, 1726, 446, 1470, 958, 1982, 126, 1150, 638, 1662, 382, 1406, 894, 1918, 254, 1278, 766, 1790, 510, 1534, 1022, 2046, 1, 1025, 513, 1537, 257, 1281, 769, 1793, 129, 1153, 641, 1665, 385, 1409, 897, 1921, 65, 1089, 577, 1601, 321, 1345, 833, 1857, 193, 1217, 705, 1729, 449, 1473, 961, 1985, 33, 1057, 545, 1569, 289, 1313, 801, 1825, 161, 1185, 673, 1697, 417, 1441, 929, 1953, 97, 1121, 609, 1633, 353, 1377, 865, 1889, 225, 1249, 737, 1761, 481, 1505, 993, 2017, 17, 1041, 529, 1553, 273, 1297, 785, 1809, 145, 1169, 657, 1681, 401, 1425, 913, 1937, 81, 1105, 593, 1617, 337, 1361, 849, 1873, 209, 1233, 721, 1745, 465, 1489, 977, 2001, 49, 1073, 561, 1585, 305, 1329, 817, 1841, 177, 1201, 689, 1713, 433, 1457, 945, 1969, 113, 1137, 625, 1649, 369, 1393, 881, 1905, 241, 1265, 753, 1777, 497, 1521, 1009, 2033, 9, 1033, 521, 1545, 265, 1289, 777, 1801, 137, 1161, 649, 1673, 393, 1417, 905, 1929, 73, 1097, 585, 1609, 329, 1353, 841, 1865, 201, 1225, 713, 1737, 457, 1481, 969, 1993, 41, 1065, 553, 1577, 297, 1321, 809, 1833, 169, 1193, 681, 1705, 425, 1449, 937, 1961, 105, 1129, 617, 1641, 361, 1385, 873, 1897, 233, 1257, 745, 1769, 489, 1513, 1001, 2025, 25, 1049, 537, 1561, 281, 1305, 793, 1817, 153, 1177, 665, 1689, 409, 1433, 921, 1945, 89, 1113, 601, 1625, 345, 1369, 857, 1881, 217, 1241, 729, 1753, 473, 1497, 985, 2009, 57, 1081, 569, 1593, 313, 1337, 825, 1849, 185, 1209, 697, 1721, 441, 1465, 953, 1977, 121, 1145, 633, 1657, 377, 1401, 889, 1913, 249, 1273, 761, 1785, 505, 1529, 1017, 2041, 5, 1029, 517, 1541, 261, 1285, 773, 1797, 133, 1157, 645, 1669, 389, 1413, 901, 1925, 69, 1093, 581, 1605, 325, 1349, 837, 1861, 197, 1221, 709, 1733, 453, 1477, 965, 1989, 37, 1061, 549, 1573, 293, 1317, 805, 1829, 165, 1189, 677, 1701, 421, 1445, 933, 1957, 101, 1125, 613, 1637, 357, 1381, 869, 1893, 229, 1253, 741, 1765, 485, 1509, 997, 2021, 21, 1045, 533, 1557, 277, 1301, 789, 1813, 149, 1173, 661, 1685, 405, 1429, 917, 1941, 85, 1109, 597, 1621, 341, 1365, 853, 1877, 213, 1237, 725, 1749, 469, 1493, 981, 2005, 53, 1077, 565, 1589, 309, 1333, 821, 1845, 181, 1205, 693, 1717, 437, 1461, 949, 1973, 117, 1141, 629, 1653, 373, 1397, 885, 1909, 245, 1269, 757, 1781, 501, 1525, 1013, 2037, 13, 1037, 525, 1549, 269, 1293, 781, 1805, 141, 1165, 653, 1677, 397, 1421, 909, 1933, 77, 1101, 589, 1613, 333, 1357, 845, 1869, 205, 1229, 717, 1741, 461, 1485, 973, 1997, 45, 1069, 557, 1581, 301, 1325, 813, 1837, 173, 1197, 685, 1709, 429, 1453, 941, 1965, 109, 1133, 621, 1645, 365, 1389, 877, 1901, 237, 1261, 749, 1773, 493, 1517, 1005, 2029, 29, 1053, 541, 1565, 285, 1309, 797, 1821, 157, 1181, 669, 1693, 413, 1437, 925, 1949, 93, 1117, 605, 1629, 349, 1373, 861, 1885, 221, 1245, 733, 1757, 477, 1501, 989, 2013, 61, 1085, 573, 1597, 317, 1341, 829, 1853, 189, 1213, 701, 1725, 445, 1469, 957, 1981, 125, 1149, 637, 1661, 381, 1405, 893, 1917, 253, 1277, 765, 1789, 509, 1533, 1021, 2045, 3, 1027, 515, 1539, 259, 1283, 771, 1795, 131, 1155, 643, 1667, 387, 1411, 899, 1923, 67, 1091, 579, 1603, 323, 1347, 835, 1859, 195, 1219, 707, 1731, 451, 1475, 963, 1987, 35, 1059, 547, 1571, 291, 1315, 803, 1827, 163, 1187, 675, 1699, 419, 1443, 931, 1955, 99, 1123, 611, 1635, 355, 1379, 867, 1891, 227, 1251, 739, 1763, 483, 1507, 995, 2019, 19, 1043, 531, 1555, 275, 1299, 787, 1811, 147, 1171, 659, 1683, 403, 1427, 915, 1939, 83, 1107, 595, 1619, 339, 1363, 851, 1875, 211, 1235, 723, 1747, 467, 1491, 979, 2003, 51, 1075, 563, 1587, 307, 1331, 819, 1843, 179, 1203, 691, 1715, 435, 1459, 947, 1971, 115, 1139, 627, 1651, 371, 1395, 883, 1907, 243, 1267, 755, 1779, 499, 1523, 1011, 2035, 11, 1035, 523, 1547, 267, 1291, 779, 1803, 139, 1163, 651, 1675, 395, 1419, 907, 1931, 75, 1099, 587, 1611, 331, 1355, 843, 1867, 203, 1227, 715, 1739, 459, 1483, 971, 1995, 43, 1067, 555, 1579, 299, 1323, 811, 1835, 171, 1195, 683, 1707, 427, 1451, 939, 1963, 107, 1131, 619, 1643, 363, 1387, 875, 1899, 235, 1259, 747, 1771, 491, 1515, 1003, 2027, 27, 1051, 539, 1563, 283, 1307, 795, 1819, 155, 1179, 667, 1691, 411, 1435, 923, 1947, 91, 1115, 603, 1627, 347, 1371, 859, 1883, 219, 1243, 731, 1755, 475, 1499, 987, 2011, 59, 1083, 571, 1595, 315, 1339, 827, 1851, 187, 1211, 699, 1723, 443, 1467, 955, 1979, 123, 1147, 635, 1659, 379, 1403, 891, 1915, 251, 1275, 763, 1787, 507, 1531, 1019, 2043, 7, 1031, 519, 1543, 263, 1287, 775, 1799, 135, 1159, 647, 1671, 391, 1415, 903, 1927, 71, 1095, 583, 1607, 327, 1351, 839, 1863, 199, 1223, 711, 1735, 455, 1479, 967, 1991, 39, 1063, 551, 1575, 295, 1319, 807, 1831, 167, 1191, 679, 1703, 423, 1447, 935, 1959, 103, 1127, 615, 1639, 359, 1383, 871, 1895, 231, 1255, 743, 1767, 487, 1511, 999, 2023, 23, 1047, 535, 1559, 279, 1303, 791, 1815, 151, 1175, 663, 1687, 407, 1431, 919, 1943, 87, 1111, 599, 1623, 343, 1367, 855, 1879, 215, 1239, 727, 1751, 471, 1495, 983, 2007, 55, 1079, 567, 1591, 311, 1335, 823, 1847, 183, 1207, 695, 1719, 439, 1463, 951, 1975, 119, 1143, 631, 1655, 375, 1399, 887, 1911, 247, 1271, 759, 1783, 503, 1527, 1015, 2039, 15, 1039, 527, 1551, 271, 1295, 783, 1807, 143, 1167, 655, 1679, 399, 1423, 911, 1935, 79, 1103, 591, 1615, 335, 1359, 847, 1871, 207, 1231, 719, 1743, 463, 1487, 975, 1999, 47, 1071, 559, 1583, 303, 1327, 815, 1839, 175, 1199, 687, 1711, 431, 1455, 943, 1967, 111, 1135, 623, 1647, 367, 1391, 879, 1903, 239, 1263, 751, 1775, 495, 1519, 1007, 2031, 31, 1055, 543, 1567, 287, 1311, 799, 1823, 159, 1183, 671, 1695, 415, 1439, 927, 1951, 95, 1119, 607, 1631, 351, 1375, 863, 1887, 223, 1247, 735, 1759, 479, 1503, 991, 2015, 63, 1087, 575, 1599, 319, 1343, 831, 1855, 191, 1215, 703, 1727, 447, 1471, 959, 1983, 127, 1151, 639, 1663, 383, 1407, 895, 1919, 255, 1279, 767, 1791, 511, 1535, 1023, 2047}; + +ulong ldn = 11; + + +Complex* CFFT::convolutionF(const Complex *input, const Complex *filter, long nSIG, long NFIL, long &NFFT) +{ + //Check for invalid inputs. + if(input == NULL || filter == NULL) + { + cout << "Could not perform convolution on empty arrays!" << endl; + return NULL; + } + + bool NFFTChanged = false; + //If NFFT not a power of 2, or it is smaller than signal or filter, prompt for new. + while (log(NFFT) / log(2) != (int)(log(NFFT) / log(2)) || NFFT < nSIG || NFFT < NFIL) + { + cout << "Please input a valid NFFT, which is >= nSIG(" << nSIG << ") and >= NFIL(" << NFIL <<") : "; + cin >> NFFT; + NFFTChanged = true; + } + + //Perform zero padding. + Complex *fInput, *fFilter; + + fInput = new Complex[NFFT]; + for(int i = 0; i < nSIG; i++) + fInput[i] = input[i]; + + fFilter = new Complex[NFFT]; + for(int i = 0; i < NFIL; i++) + fFilter[i] = filter[i]; + + //Store the output data. + Complex *output = new Complex[NFFT]; + + //Perform FFT on both input and filter. + CFFT::Forward(fInput, (unsigned int)NFFT); + CFFT::Forward(fFilter, (unsigned int)NFFT); + + for(int i = 0; i < NFFT; i++) + output[i] = fInput[i] * fFilter[i];// SIMD? + + return output; +} + +//Convolution function which returns the time representation of the result. +//NFFT is the FFT size, nSIG is the size for the input, NFIL is the size of the filter. +Complex* CFFT::convolutionT(const Complex *input,const Complex *filter, long nSIG, long NFIL, long &NFFT) +{ + //Store the output data. + Complex *output = convolutionF(input, filter, nSIG, NFIL, NFFT); + + //Perform IFFT on the ouput. + CFFT::Inverse(output, (unsigned int)NFFT); + + return output; +} + +Complex* CFFT::stereoConvMonoInputF(const Complex *input,const Complex *filterLeft,const Complex *filterRight, long nSIG, long NFILL, long NFILR, long &NFFT) +{ + Complex *result = stereoConvMonoInputT(input, filterLeft, filterRight, nSIG, NFILL, NFILR, NFFT); + CFFT::Forward(result, (unsigned int)NFFT); + return result; +} + +//The size of the ouput will be 2 times of the size of FFT for the input signal and the NFFT value +//will be doubled after running the function +Complex* CFFT::stereoConvMonoInputT(const Complex *input,const Complex *filterLeft,const Complex *filterRight, long nSIG, long NFILL, long NFilR, long &NFFT) +{ + Complex *tempLeft = new Complex[NFFT]; + Complex *tempRight = new Complex[NFFT]; + Complex *result = new Complex[2*NFFT]; + + tempLeft = CFFT::convolutionT(input, filterLeft, nSIG, NFILL, NFFT); + tempRight = CFFT::convolutionT(input, filterRight, nSIG, NFILL, NFFT); + NFFT = NFFT * 2; + for (int i = 0; i < NFFT / 2; i++) + { + result[2 * i] = tempLeft[i]; + result[2 * i + 1] = tempRight[i]; + } + + delete tempLeft; + delete tempRight; + + return result; +} + +/* + Complex* CFFT::stereoConvStereoInputT(const Complex *input, const Complex *filterLeft, const Complex *filterRight, long nSIG, long NFILL, long NFILR, long &NFFT) + { + NFFT = NFFT / 2; + Complex *leftTemp = new Complex[NFFT]; + Complex *rightTemp = new Complex[NFFT]; + + + + + + return 0; + + } + + */ + + +//storing the an array into a text file +//filename is the file name you want to store the data into +//datatype represents the data you wanna store: real/real+imag/amplitude +void CFFT::storingData(Complex *data, int NFFT,string temp ,char datatype) +{ + //string temp= filename; + ofstream outputFile(temp.c_str()); + if (outputFile.is_open()) + { + switch (datatype) + { + case 'r': + //ofstream outputFile("real.txt"); + for (int i = 0; i < NFFT; i++) + outputFile << data[i].real() << endl; + break; + + case 'c': + for (int i = 0; i < NFFT; i++) + outputFile << data[i].real() << " " << data[i].imag() << endl; + break; + case 'a': + for (int i = 0; i < NFFT; i++) + outputFile << sqrt(data[i].real()*data[i].real()+data[i].imag()*data[i].imag()) << endl; + break; + } + outputFile.close(); + } +} + + + +// FORWARD FOURIER TRANSFORM +// Input - input data +// Output - transform result +// N - length of both input data and result +bool CFFT::Forward(const Complex *const Input, Complex *const Output, const unsigned int N) +{ + // Check input parameters + if (!Input || !Output || N < 1 || N & (N - 1)) + return false; + // Initialize data + Rearrange(Input, Output, N); + // Call FFT implementation + Perform(Output, N); + // Succeeded + return true; +} + +// FORWARD FOURIER TRANSFORM, INPLACE VERSION +// Data - both input data and output +// N - length of input data +bool CFFT::Forward(Complex *const Data, const unsigned int N) +{ + // Check input parameters + if (!Data || N < 1 || N & (N - 1)) + return false; + // Rearrange + Rearrange(Data, N); + // Call FFT implementation + Perform(Data, N); + // Succeeded + return true; +} + +// INVERSE FOURIER TRANSFORM +// Input - input data +// Output - transform result +// N - length of both input data and result +// Scale - if to scale result +bool CFFT::Inverse(const Complex *const Input, Complex *const Output, const unsigned int N, const bool Scale /* = true */) +{ + // Check input parameters + if (!Input || !Output || N < 1 || N & (N - 1)) + return false; + // Initialize data + Rearrange(Input, Output, N); + // Call FFT implementation + Perform(Output, N, true); + // Scale if necessary + if (Scale) + CFFT::Scale(Output, N); + // Succeeded + return true; +} + +// INVERSE FOURIER TRANSFORM, INPLACE VERSION +// Data - both input data and output +// N - length of both input data and result +// Scale - if to scale result +bool CFFT::Inverse(Complex *const Data, const unsigned int N, const bool Scale /* = true */) +{ + // Check input parameters + if (!Data || N < 1 || N & (N - 1)) + return false; + // Rearrange + Rearrange(Data, N); + // Call FFT implementation + Perform(Data, N, true); + // Scale if necessary + if (Scale) + CFFT::Scale(Data, N); + // Succeeded + return true; +} + +// Rearrange function +void CFFT::Rearrange(const Complex *const Input, Complex *const Output, const unsigned int N) +{ + // // Data entry position + // unsigned int Target = 0; + // // Process all positions of input signal + for (unsigned int Position = 0; Position < N; ++Position) + { + Output[Position] = Input[Position]; + }// SIMD? +} + +// Inplace version of rearrange function +void CFFT::Rearrange(Complex *const Data, const unsigned int N){} + + +// FFT implementation +void CFFT::Perform(Complex *const Data, const unsigned int N, const bool Inverse /* = false */) +{ + + if (Inverse) { + // Inverse transform from Complex number to real number(Complex number with real part defined and imaginary part set to zero) radix-4 decimation in Time method + double dataReal [N]; + double dataImag [N]; + + for (int i = 0; i < N; i++) { + //cout << "i: " << i << endl; + //cout << "getPermutedIndex(i): " << getPermutedIndex(i) << endl; + dataReal[i] = Data[PermutedIndex[i]].real(); + dataImag[i] = Data[PermutedIndex[i]].imag(); + } + double *Kr = dataReal; + double *Ki = dataImag; + + + //revbin_permute(Kr, 1UL<<ldn); + //revbin_permute(Ki, 1UL<<ldn); + + fft_dit4_core_p1(Kr, Ki, ldn); + + for (unsigned int Position = 0; Position < N; ++Position){ + Data[Position] = std::complex <double> (Kr[Position], 0); + } + + + } + else{ + // forward transform from real number to Complex number using FHT method + // first permute Data + double *f = new double [N]; + for (int i = 0; i < N; i++) { + f[i] = Data[i].real(); + } + revbin_permute(f, 1UL<<ldn); + // now tmpdata is the permuted version of Data + fht_loc_dit2_core(f, ldn); + ///////////// FHT to FFT conversion start /////////////// + const ulong n = (1UL<<ldn); + //int is = 1; + //if ( is>0 ) for (ulong i=1, j=n-1; i<j; i++, j--) sumdiff05(f[i], f[j]); + //else for (ulong i=1, j=n-1; i<j; i++, j--) sumdiff05_r(f[i], f[j]); + for (ulong i=1, j=n-1; i<j; i++, j--) sumdiff05(f[i], f[j]); + + for (int i = 1; i<512; i++) { + Complex tmp(f[i],-f[N-i]); + Data[i] = tmp; + } + Data[512] = Complex(f[512],0); + for (int i = 513; i<N; i++) { + Complex tmp(f[N-i],f[i]); + Data[i] = tmp; + } + Data[0] = Complex(f[0],0); + delete f; + ///////////// FHT to FFT conversion finished /////////////// + } + +} + +// Scaling of inverse FFT result +void CFFT::Scale(Complex *const Data, const unsigned int N) +{ + const double Factor = 1. / double(N); + // Scale all data entries + for (unsigned int Position = 0; Position < N; ++Position) + Data[Position] *= Factor; +}// SIMD?
diff -r 000000000000 -r 2076b4d80327 fhtdit.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fhtdit.cc Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,262 @@ +// This file is part of the FXT library. +// Copyright (C) 2010, 2012 Joerg Arndt +// License: GNU General Public License version 3 or later, +// see the file COPYING.txt in the main directory. + +#include "fxttypes.h" +#include "cmult.h" +#include "sumdiff.h" +#include "constants.h" // COS_1_PI_8, SIN_1_PI_8 +#include "sincos.h" +#include "revbinpermute.h" + +#include <cmath> // M_PI, M_SQRT1_2 + +// tuning parameter: +// define to use trig recurrence: +// (and possibly lose some precision, see below) +//#define USE_TRIG_REC +// with type 'long double' slight speed loss on my machine, +// with type 'double' little speed gain. + +#if defined USE_TRIG_REC +//#warning 'FYI: fht(double *, ulong) uses trig recursion' +#endif + +// tuning parameter: +#define INITIAL_RADIX_16 1 // 0 or 1 (default) +// +#if ( INITIAL_RADIX_16==1 ) +//#warning 'FYI: INITIAL_RADIX_16 set in fht_dit(double *, ulong)' +#else +#warning 'INITIAL_RADIX_16 is NOT SET in fht_dit(double *, ulong)' +#endif + + +void +fht_dit_core(double *f, ulong ldn) +// Decimation in time (DIT) FHT. +// Input data must be in revbin_permuted order. +// ldn := base-2 logarithm of the array length. +{ + if ( ldn<=2 ) + { + if ( ldn==1 ) // two point fht + { + sumdiff(f[0], f[1]); + } + else if ( ldn==2 ) // four point fht + { + double f0, f1, f2, f3; + sumdiff(f[0], f[1], f0, f1); + sumdiff(f[2], f[3], f2, f3); + sumdiff(f0, f2, f[0], f[2]); + sumdiff(f1, f3, f[1], f[3]); + } + return; + } + + const ulong n = (1UL<<ldn); + const double *fn = f + n; + ulong ldk = ldn & 1; + + if ( ldk==0 ) // ldn is multiple of 2 => n is a power of 4 + { +#if ( INITIAL_RADIX_16==1 ) + for (double *fi=f; fi<fn; fi+=16) // radix-16 step + { + double f0, f1, f2, f3; + sumdiff(fi[0], fi[1], f0, f1); + sumdiff(fi[2], fi[3], f2, f3); + sumdiff(f0, f2, fi[0], fi[2]); + sumdiff(f1, f3, fi[1], fi[3]); + + sumdiff(fi[4], fi[5], f0, f1); + sumdiff(fi[6], fi[7], f2, f3); + sumdiff(f0, f2, fi[4], fi[6]); + sumdiff(f1, f3, fi[5], fi[7]); + + sumdiff(fi[8], fi[9], f0, f1); + sumdiff(fi[10], fi[11], f2, f3); + sumdiff(f0, f2, fi[8], fi[10]); + sumdiff(f1, f3, fi[9], fi[11]); + + sumdiff(fi[12], fi[13], f0, f1); + sumdiff(fi[14], fi[15], f2, f3); + sumdiff(f0, f2, fi[12], fi[14]); + sumdiff(f1, f3, fi[13], fi[15]); + + sumdiff(fi[0], fi[4], f0, f1); + sumdiff(fi[8], fi[12], f2, f3); + sumdiff(f0, f2, fi[0], fi[8]); + sumdiff(f1, f3, fi[4], fi[12]); + + sumdiff(fi[2], fi[6], f0, f1); + f3 = M_SQRT2 * fi[14]; + f2 = M_SQRT2 * fi[10]; + sumdiff(f0, f2, fi[2], fi[10]); + sumdiff(f1, f3, fi[6], fi[14]); + + double a, b, g0, g1, g2, g3; + sumdiff(fi[5], fi[7], a, b); + a *= M_SQRT1_2; + b *= M_SQRT1_2; + sumdiff(fi[1], a, f0, f1); + sumdiff(fi[3], b, g0, g1); + sumdiff(fi[13], fi[15], a, b); + a *= M_SQRT1_2; + b *= M_SQRT1_2; + sumdiff(fi[9], a, f2, f3); + sumdiff(fi[11], b, g2, g3); + double c1 = COS_1_PI_8; // jjkeep + double s1 = SIN_1_PI_8; // jjkeep + cmult(s1, c1, f2, g3, b, a); + sumdiff(f0, a, fi[1], fi[9]); + sumdiff(g1, b, fi[7], fi[15]); + cmult(c1, s1, g2, f3, b, a); + sumdiff(g0, a, fi[3], fi[11]); + sumdiff(f1, b, fi[5], fi[13]); + } + ldk = 4; +#else // INITIAL_RADIX_16 + for (double *fi=f; fi<fn; fi+=4) // radix-4 step + { + double f0, f1, f2, f3; + sumdiff(fi[0], fi[1], f0, f1); + sumdiff(fi[2], fi[3], f2, f3); + sumdiff(f0, f2, fi[0], fi[2]); + sumdiff(f1, f3, fi[1], fi[3]); + } + ldk = 2; +#endif // INITIAL_RADIX_16 + } + else // ldk==1, n is no power of 4 + { + for (double *fi=f; fi<fn; fi+=8) // radix-8 step + { + double g0, f0, f1, g1; + sumdiff(fi[0], fi[1], f0, g0); + sumdiff(fi[2], fi[3], f1, g1); + + sumdiff(f0, f1); + sumdiff(g0, g1); + + double s1, c1, s2, c2; + sumdiff(fi[4], fi[5], s1, c1); + sumdiff(fi[6], fi[7], s2, c2); + + sumdiff(s1, s2); + + sumdiff(f0, s1, fi[0], fi[4]); + sumdiff(f1, s2, fi[2], fi[6]); + + c1 *= M_SQRT2; + c2 *= M_SQRT2; + + sumdiff(g0, c1, fi[1], fi[5]); + sumdiff(g1, c2, fi[3], fi[7]); + } + ldk = 3; + } + + + for ( ; ldk<ldn; ldk+=2) + { + ulong k = 1UL<<ldk; + ulong kh = k >> 1; + ulong k2 = k << 1; + ulong k3 = k2 + k; + ulong k4 = k2 << 1; + + for (double *fi=f, *gi=f+kh; fi<fn; fi+=k4, gi+=k4) + { + double f0, f1, f2, f3; + sumdiff(fi[0], fi[k], f0, f1); + sumdiff(fi[k2], fi[k3], f2, f3); + sumdiff(f0, f2, fi[0], fi[k2]); + sumdiff(f1, f3, fi[k], fi[k3]); + + sumdiff(gi[0], gi[k], f0, f1); + f3 = M_SQRT2 * gi[k3]; + f2 = M_SQRT2 * gi[k2]; + sumdiff(f0, f2, gi[0], gi[k2]); + sumdiff(f1, f3, gi[k], gi[k3]); + } + + double tt = M_PI_4/(double)kh; // jjkeep +#if defined USE_TRIG_REC + double s1 = 0.0, c1 = 1.0; // jjkeep + double al = sin(0.5*tt); // jjkeep + al *= (2.0*al); + double be = sin(tt); // jjkeep +#endif // USE_TRIG_REC + + for (ulong i=1; i<kh; i++) + { +#if defined USE_TRIG_REC + { + double t1 = c1; // jjkeep + c1 -= (al*t1+be*s1); + s1 -= (al*s1-be*t1); + } +#else + double s1, c1; // jjkeep + SinCos(tt*(double)i, &s1, &c1); // jjkeep +#endif // USE_TRIG_REC + + double c2, s2; // jjkeep + csqr(c1, s1, c2, s2); + + for (double *fi=f+i, *gi=f+k-i; fi<fn; fi+=k4, gi+=k4) + { + double a, b, g0, f0, f1, g1, f2, g2, f3, g3; + cmult(s2, c2, fi[k], gi[k], b, a); + sumdiff(fi[0], a, f0, f1); + sumdiff(gi[0], b, g0, g1); + + cmult(s2, c2, fi[k3], gi[k3], b, a); + sumdiff(fi[k2], a, f2, f3); + sumdiff(gi[k2], b, g2, g3); + + cmult(s1, c1, f2, g3, b, a); + sumdiff(f0, a, fi[0], fi[k2]); + sumdiff(g1, b, gi[k], gi[k3]); + + cmult(c1, s1, g2, f3, b, a); + sumdiff(g0, a, gi[0], gi[k2]); + sumdiff(f1, b, fi[k], fi[k3]); + } + } + } +} +// ------------------------- + + +void +fht_dit(double *f, ulong ldn) +// Fast Hartley Transform. +// Split-radix decimation in time (DIT) algorithm. +// ldn := base-2 logarithm of the array length. +{ + if ( ldn<=2 ) + { + if ( ldn==1 ) // two point fht + { + sumdiff(f[0], f[1]); + } + else if ( ldn==2 ) // four point fht + { + double f0, f1, f2, f3; + sumdiff(f[0], f[2], f0, f1); + sumdiff(f[1], f[3], f2, f3); + sumdiff(f0, f2, f[0], f[2]); + sumdiff(f1, f3, f[1], f[3]); + } + + return; + } + + revbin_permute(f, 1UL<<ldn); + fht_dit_core(f, ldn); +} +// -------------------------
diff -r 000000000000 -r 2076b4d80327 fhtdit2.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/fhtdit2.cc Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,130 @@ +// This file is part of the FXT library. +// Copyright (C) 2010, 2012 Joerg Arndt +// License: GNU General Public License version 3 or later, +// see the file COPYING.txt in the main directory. + +#include "fxttypes.h" +#include "complextype.h" +#include "sincos.h" +#include "revbinpermute.h" +#include "sumdiff.h" + +#include <cmath> // M_PI + + +void +fht_depth_first_dit2(double *f, ulong ldn) +// Radix-2 decimation in time (DIT) FHT. +// Depth-first version. +// Compared to usual fht this one +// - does more trig computations +// - is (far) better memory local +{ + const ulong n = 1UL<<ldn; + revbin_permute(f, n); + + for (ulong ldm=1; ldm<=ldn; ++ldm) + { + const ulong m = (1UL<<ldm); + const ulong mh = (m>>1); + const ulong m4 = (mh>>1); + const double phi0 = M_PI/(double)mh; + + for (ulong r=0; r<n; r+=m) + { + { // j == 0: + ulong t1 = r; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + + if ( m4 ) + { + ulong t1 = r + m4; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + + for (ulong j=1, k=mh-1; j<k; ++j, --k) + { + double s, c; + SinCos(phi0*(double)j, &s, &c); + + ulong tj = r + mh + j; + ulong tk = r + mh + k; + double fj = f[tj]; + double fk = f[tk]; + f[tj] = fj * c + fk * s; + f[tk] = fj * s - fk * c; + + ulong t1 = r + j; + ulong t2 = tj; // == t1 + mh; + sumdiff(f[t1], f[t2]); + + t1 = r + k; + t2 = tk; // == t1 + mh; + sumdiff(f[t1], f[t2]); + } + } + } +} +// ------------------------- + + +void +fht_dit2(double *f, ulong ldn) +// Radix-2 decimation in time (DIT) FHT. +{ + const ulong n = 1UL<<ldn; + revbin_permute(f, n); + + for (ulong ldm=1; ldm<=ldn; ++ldm) + { + const ulong m = (1UL<<ldm); + const ulong mh = (m>>1); + const ulong m4 = (mh>>1); + const double phi0 = M_PI/(double)mh; + + for (ulong r=0; r<n; r+=m) + { + { // j == 0: + ulong t1 = r; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + + if ( m4 ) + { + ulong t1 = r + m4; + ulong t2 = t1 + mh; + sumdiff(f[t1], f[t2]); + } + } + + for (ulong j=1, k=mh-1; j<k; ++j, --k) + { + double s, c; + SinCos(phi0*(double)j, &s, &c); + + for (ulong r=0; r<n; r+=m) + { + ulong tj = r + mh + j; + ulong tk = r + mh + k; + double fj = f[tj]; + double fk = f[tk]; + f[tj] = fj * c + fk * s; + f[tk] = fj * s - fk * c; + + ulong t1 = r + j; + ulong t2 = tj; // == t1 + mh; + sumdiff(f[t1], f[t2]); + + t1 = r + k; + t2 = tk; // == t1 + mh; + sumdiff(f[t1], f[t2]); + } + } + } +} +// ------------------------- +
diff -r 000000000000 -r 2076b4d80327 location.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/location.cpp Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,54 @@ +#include "../include/location.h" +#include <iostream> + +Location& Location::operator= (const Location& rhs) +{ + x = (rhs.x); + y = (rhs.y); + z = (rhs.z); + return *this; +} + +Location& Location::operator+= (const Location& rhs) +{ + x = x + rhs.x; + y = y + rhs.y; + z = z + rhs.z; + return *this; +} + +Location operator+ (Location lhs, const Location& rhs) +{ + lhs += rhs; + return lhs; +} + +bool Location::operator< (const Location& rhs) const +{ + if (x == rhs.x){ + if (y == rhs.y){ + return z < rhs.z; + }return y < rhs.y; + }return x < rhs.x; +} + +ostream& operator<< (ostream& os, Location obj) +{ + os << obj.getX() <<","<< obj.getY() <<","<< obj.getZ(); + return os; +} + +float Location::getX (void) const +{ + return this->x; +} + +float Location::getY (void) const +{ + return this->y; +} + +float Location::getZ (void) const +{ + return this->z; +}
diff -r 000000000000 -r 2076b4d80327 mit_hrtf_lib.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mit_hrtf_lib.c Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,467 @@ +/*############################################################################*/ +/*# #*/ +/*# MIT HRTF C Library #*/ +/*# Copyright � 2007 Aristotel Digenis #*/ +/*# #*/ +/*# Filename: mit_hrtf_lib.c #*/ +/*# Version: 0.1 #*/ +/*# Date: 04/05/2007 #*/ +/*# Author(s): Aristotel Digenis (adigenis@users.sourceforge.net) #*/ +/*# Credit: Bill Gardner and Keith Martin #*/ +/*# Licence: GNU Library or Lesser General Public License (LGPL) #*/ +/*# #*/ +/*############################################################################*/ + + +#include "../include/mit_hrtf_lib.h" +#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_44100.h" +#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_44100.h" +#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_48000.h" +#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_48000.h" +#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_88200.h" +#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_88200.h" +#include "../MIT_HRTF_Library/source/normal/mit_hrtf_normal_96000.h" +#include "../MIT_HRTF_Library/source/diffuse/mit_hrtf_diffuse_96000.h" + + +/* Internal functions for handling the indexing of the -/+40 degree elevation + data which has irregular azimuth positions. */ +int mit_hrtf_findAzimuthFor40Elev(int azimuth); +int mit_hrtf_findIndexFor40Elev(int azimuth); + + + +unsigned int mit_hrtf_availability(int azimuth, int elevation, unsigned int samplerate, unsigned int diffused) +{ + if(elevation > 90 || elevation < -40) + return 0; + + if(azimuth > 180 || azimuth < -180) + return 0; + + if(diffused > 1) + return 0; + + if(samplerate == 44100) + return MIT_HRTF_44_TAPS; + else if(samplerate == 48000) + return MIT_HRTF_48_TAPS; + else if(samplerate == 88200) + return MIT_HRTF_88_TAPS; + else if(samplerate == 96000) + return MIT_HRTF_96_TAPS; + + return 0; +} + + + +unsigned int mit_hrtf_get(int* pAzimuth, int* pElevation, unsigned int samplerate, unsigned int diffused, short* psLeft, short* psRight) +{ + int nInternalElevation = 0; + float fAzimuthIncrement = 0; + int nInternalAzimuth = 0; + int nSwitchLeftRight = 0; + int nAzimuthIndex = 0; + const mit_hrtf_filter_set_44* pFilter44 = 0; + const mit_hrtf_filter_set_48* pFilter48 = 0; + const mit_hrtf_filter_set_88* pFilter88 = 0; + const mit_hrtf_filter_set_96* pFilter96 = 0; + const short* psLeftTaps = 0; + const short* psRightTaps = 0; + const short* psTempTaps = 0; + unsigned int nTotalTaps = 0; + unsigned int niTap = 0; + + //Check if the requested HRTF exists + if(!mit_hrtf_availability(*pAzimuth, *pElevation, samplerate, diffused)) + return 0; + + //Snap elevation to the nearest available elevation in the filter set + if(*pElevation < 0) + nInternalElevation = ((*pElevation - 5) / 10) * 10; + else + nInternalElevation = ((*pElevation + 5) / 10) * 10; + + // Elevation of 50 has a maximum 176 in the azimuth plane so we need to handle that. + if(nInternalElevation == 50) + { + if(*pAzimuth < 0) + *pAzimuth = *pAzimuth < -176 ? -176 : *pAzimuth; + else + *pAzimuth = *pAzimuth > 176 ? 176 : *pAzimuth; + } + + //Snap azimuth to the nearest available azimuth in the filter set. + switch(nInternalElevation) + { + case 0: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_00 - 1); break; // 180 5 + case 10: + case -10: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_10 - 1); break; // 180 5 + case 20: + case -20: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_20 - 1); break; // 180 5 + case 30: + case -30: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_30 - 1); break; // 180 6 + case 40: + case -40: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_40 - 1); break; // 180 6.43 + case 50: fAzimuthIncrement = 176.f / (MIT_HRTF_AZI_POSITIONS_50 - 1); break; // 176 8 + case 60: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_60 - 1); break; // 180 10 + case 70: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_70 - 1); break; // 180 15 + case 80: fAzimuthIncrement = 180.f / (MIT_HRTF_AZI_POSITIONS_80 - 1); break; // 180 30 + case 90: fAzimuthIncrement = 0; break; // 0 1 + }; + + if(*pAzimuth < 0) + { + nInternalAzimuth = (int)((int)((-*pAzimuth + fAzimuthIncrement / 2.f) / fAzimuthIncrement) * fAzimuthIncrement + 0.5f); + nSwitchLeftRight = 1; + } + else + { + nInternalAzimuth = (int)((int)((*pAzimuth + fAzimuthIncrement / 2.f) / fAzimuthIncrement) * fAzimuthIncrement + 0.5f); + } + + //Determine array index for azimuth based on elevation + switch(nInternalElevation) + { + case 0: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_00 - 1)); break; + case 10: + case -10: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_10 - 1)); break; + case 20: + case -20: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_20 - 1)); break; + case 30: + case -30: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_30 - 1)); break; + case 40: + case -40: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_40 - 1)); break; + case 50: nAzimuthIndex = (int)((nInternalAzimuth / 176.f) * (MIT_HRTF_AZI_POSITIONS_50 - 1)); break; + case 60: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_60 - 1)); break; + case 70: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_70 - 1)); break; + case 80: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_80 - 1)); break; + case 90: nAzimuthIndex = (int)((nInternalAzimuth / 180.f) * (MIT_HRTF_AZI_POSITIONS_90 - 1)); break; + }; + + //The azimuths for +/- elevations need special handling + if(nInternalElevation == 40 || nInternalElevation == -40) + { + nInternalAzimuth = mit_hrtf_findAzimuthFor40Elev(nInternalAzimuth); + nAzimuthIndex = mit_hrtf_findIndexFor40Elev(nInternalAzimuth); + } + + //Assign pointer to appropriate array depending on saple rate, normal or diffused filters, elevation, and azimuth index. + switch(samplerate) + { + case 44100: + if(diffused) + pFilter44 = &diffuse_44; + else + pFilter44 = &normal_44; + + switch(nInternalElevation) + { + case -10: psLeftTaps = pFilter44->e_10[nAzimuthIndex].left; + psRightTaps = pFilter44->e_10[nAzimuthIndex].right; break; + case -20: psLeftTaps = pFilter44->e_20[nAzimuthIndex].left; + psRightTaps = pFilter44->e_20[nAzimuthIndex].right; break; + case -30: psLeftTaps = pFilter44->e_30[nAzimuthIndex].left; + psRightTaps = pFilter44->e_30[nAzimuthIndex].right; break; + case -40: psLeftTaps = pFilter44->e_40[nAzimuthIndex].left; + psRightTaps = pFilter44->e_40[nAzimuthIndex].right; break; + case 0: psLeftTaps = pFilter44->e00[nAzimuthIndex].left; + psRightTaps = pFilter44->e00[nAzimuthIndex].right; break; + case 10: psLeftTaps = pFilter44->e10[nAzimuthIndex].left; + psRightTaps = pFilter44->e10[nAzimuthIndex].right; break; + case 20: psLeftTaps = pFilter44->e20[nAzimuthIndex].left; + psRightTaps = pFilter44->e20[nAzimuthIndex].right; break; + case 30: psLeftTaps = pFilter44->e30[nAzimuthIndex].left; + psRightTaps = pFilter44->e30[nAzimuthIndex].right; break; + case 40: psLeftTaps = pFilter44->e40[nAzimuthIndex].left; + psRightTaps = pFilter44->e40[nAzimuthIndex].right; break; + case 50: psLeftTaps = pFilter44->e50[nAzimuthIndex].left; + psRightTaps = pFilter44->e50[nAzimuthIndex].right; break; + case 60: psLeftTaps = pFilter44->e60[nAzimuthIndex].left; + psRightTaps = pFilter44->e60[nAzimuthIndex].right; break; + case 70: psLeftTaps = pFilter44->e70[nAzimuthIndex].left; + psRightTaps = pFilter44->e70[nAzimuthIndex].right; break; + case 80: psLeftTaps = pFilter44->e80[nAzimuthIndex].left; + psRightTaps = pFilter44->e80[nAzimuthIndex].right; break; + case 90: psLeftTaps = pFilter44->e90[nAzimuthIndex].left; + psRightTaps = pFilter44->e90[nAzimuthIndex].right; break; + }; + + //How many taps to copy later to user's buffers + nTotalTaps = MIT_HRTF_44_TAPS; + break; + case 48000: + if(diffused) + pFilter48 = &diffuse_48; + else + pFilter48 = &normal_48; + + switch(nInternalElevation) + { + case -10: psLeftTaps = pFilter48->e_10[nAzimuthIndex].left; + psRightTaps = pFilter48->e_10[nAzimuthIndex].right; break; + case -20: psLeftTaps = pFilter48->e_20[nAzimuthIndex].left; + psRightTaps = pFilter48->e_20[nAzimuthIndex].right; break; + case -30: psLeftTaps = pFilter48->e_30[nAzimuthIndex].left; + psRightTaps = pFilter48->e_30[nAzimuthIndex].right; break; + case -40: psLeftTaps = pFilter48->e_40[nAzimuthIndex].left; + psRightTaps = pFilter48->e_40[nAzimuthIndex].right; break; + case 0: psLeftTaps = pFilter48->e00[nAzimuthIndex].left; + psRightTaps = pFilter48->e00[nAzimuthIndex].right; break; + case 10: psLeftTaps = pFilter48->e10[nAzimuthIndex].left; + psRightTaps = pFilter48->e10[nAzimuthIndex].right; break; + case 20: psLeftTaps = pFilter48->e20[nAzimuthIndex].left; + psRightTaps = pFilter48->e20[nAzimuthIndex].right; break; + case 30: psLeftTaps = pFilter48->e30[nAzimuthIndex].left; + psRightTaps = pFilter48->e30[nAzimuthIndex].right; break; + case 40: psLeftTaps = pFilter48->e40[nAzimuthIndex].left; + psRightTaps = pFilter48->e40[nAzimuthIndex].right; break; + case 50: psLeftTaps = pFilter48->e50[nAzimuthIndex].left; + psRightTaps = pFilter48->e50[nAzimuthIndex].right; break; + case 60: psLeftTaps = pFilter48->e60[nAzimuthIndex].left; + psRightTaps = pFilter48->e60[nAzimuthIndex].right; break; + case 70: psLeftTaps = pFilter48->e70[nAzimuthIndex].left; + psRightTaps = pFilter48->e70[nAzimuthIndex].right; break; + case 80: psLeftTaps = pFilter48->e80[nAzimuthIndex].left; + psRightTaps = pFilter48->e80[nAzimuthIndex].right; break; + case 90: psLeftTaps = pFilter48->e90[nAzimuthIndex].left; + psRightTaps = pFilter48->e90[nAzimuthIndex].right; break; + }; + + //How many taps to copy later to user's buffers + nTotalTaps = MIT_HRTF_48_TAPS; + break; + case 88200: + if(diffused) + pFilter88 = &diffuse_88; + else + pFilter88 = &normal_88; + + switch(nInternalElevation) + { + case -10: psLeftTaps = pFilter88->e_10[nAzimuthIndex].left; + psRightTaps = pFilter88->e_10[nAzimuthIndex].right; break; + case -20: psLeftTaps = pFilter88->e_20[nAzimuthIndex].left; + psRightTaps = pFilter88->e_20[nAzimuthIndex].right; break; + case -30: psLeftTaps = pFilter88->e_30[nAzimuthIndex].left; + psRightTaps = pFilter88->e_30[nAzimuthIndex].right; break; + case -40: psLeftTaps = pFilter88->e_40[nAzimuthIndex].left; + psRightTaps = pFilter88->e_40[nAzimuthIndex].right; break; + case 0: psLeftTaps = pFilter88->e00[nAzimuthIndex].left; + psRightTaps = pFilter88->e00[nAzimuthIndex].right; break; + case 10: psLeftTaps = pFilter88->e10[nAzimuthIndex].left; + psRightTaps = pFilter88->e10[nAzimuthIndex].right; break; + case 20: psLeftTaps = pFilter88->e20[nAzimuthIndex].left; + psRightTaps = pFilter88->e20[nAzimuthIndex].right; break; + case 30: psLeftTaps = pFilter88->e30[nAzimuthIndex].left; + psRightTaps = pFilter88->e30[nAzimuthIndex].right; break; + case 40: psLeftTaps = pFilter88->e40[nAzimuthIndex].left; + psRightTaps = pFilter88->e40[nAzimuthIndex].right; break; + case 50: psLeftTaps = pFilter88->e50[nAzimuthIndex].left; + psRightTaps = pFilter88->e50[nAzimuthIndex].right; break; + case 60: psLeftTaps = pFilter88->e60[nAzimuthIndex].left; + psRightTaps = pFilter88->e60[nAzimuthIndex].right; break; + case 70: psLeftTaps = pFilter88->e70[nAzimuthIndex].left; + psRightTaps = pFilter88->e70[nAzimuthIndex].right; break; + case 80: psLeftTaps = pFilter88->e80[nAzimuthIndex].left; + psRightTaps = pFilter88->e80[nAzimuthIndex].right; break; + case 90: psLeftTaps = pFilter88->e90[nAzimuthIndex].left; + psRightTaps = pFilter88->e90[nAzimuthIndex].right; break; + }; + + //How many taps to copy later to user's buffers + nTotalTaps = MIT_HRTF_88_TAPS; + break; + case 96000: + if(diffused) + pFilter96 = &diffuse_96; + else + pFilter96 = &normal_96; + + switch(nInternalElevation) + { + case -10: psLeftTaps = pFilter96->e_10[nAzimuthIndex].left; + psRightTaps = pFilter96->e_10[nAzimuthIndex].right; break; + case -20: psLeftTaps = pFilter96->e_20[nAzimuthIndex].left; + psRightTaps = pFilter96->e_20[nAzimuthIndex].right; break; + case -30: psLeftTaps = pFilter96->e_30[nAzimuthIndex].left; + psRightTaps = pFilter96->e_30[nAzimuthIndex].right; break; + case -40: psLeftTaps = pFilter96->e_40[nAzimuthIndex].left; + psRightTaps = pFilter96->e_40[nAzimuthIndex].right; break; + case 0: psLeftTaps = pFilter96->e00[nAzimuthIndex].left; + psRightTaps = pFilter96->e00[nAzimuthIndex].right; break; + case 10: psLeftTaps = pFilter96->e10[nAzimuthIndex].left; + psRightTaps = pFilter96->e10[nAzimuthIndex].right; break; + case 20: psLeftTaps = pFilter96->e20[nAzimuthIndex].left; + psRightTaps = pFilter96->e20[nAzimuthIndex].right; break; + case 30: psLeftTaps = pFilter96->e30[nAzimuthIndex].left; + psRightTaps = pFilter96->e30[nAzimuthIndex].right; break; + case 40: psLeftTaps = pFilter96->e40[nAzimuthIndex].left; + psRightTaps = pFilter96->e40[nAzimuthIndex].right; break; + case 50: psLeftTaps = pFilter96->e50[nAzimuthIndex].left; + psRightTaps = pFilter96->e50[nAzimuthIndex].right; break; + case 60: psLeftTaps = pFilter96->e60[nAzimuthIndex].left; + psRightTaps = pFilter96->e60[nAzimuthIndex].right; break; + case 70: psLeftTaps = pFilter96->e70[nAzimuthIndex].left; + psRightTaps = pFilter96->e70[nAzimuthIndex].right; break; + case 80: psLeftTaps = pFilter96->e80[nAzimuthIndex].left; + psRightTaps = pFilter96->e80[nAzimuthIndex].right; break; + case 90: psLeftTaps = pFilter96->e90[nAzimuthIndex].left; + psRightTaps = pFilter96->e90[nAzimuthIndex].right; break; + }; + + //How many taps to copy later to user's buffers + nTotalTaps = MIT_HRTF_96_TAPS; + break; + }; + + //Switch left and right ear if the azimuth is to the left of front centre (azimuth < 0) + if(nSwitchLeftRight) + { + psTempTaps = psRightTaps; + psRightTaps = psLeftTaps; + psLeftTaps = psTempTaps; + } + + //Copy taps to user's arrays + for(niTap = 0; niTap < nTotalTaps; niTap++) + { + psLeft[niTap] = psLeftTaps[niTap]; + psRight[niTap] = psRightTaps[niTap]; + } + + //Assign the real azimuth and elevation used + *pAzimuth = nInternalAzimuth; + *pElevation = nInternalElevation; + + return nTotalTaps; +} + + + +int mit_hrtf_findAzimuthFor40Elev(int azimuth) +{ + if(azimuth >= 0 && azimuth < 4) + return 0; + else if(azimuth >= 4 && azimuth < 10) + return 6; + else if(azimuth >= 10 && azimuth < 17) + return 13; + else if(azimuth >= 17 && azimuth < 23) + return 19; + else if(azimuth >= 23 && azimuth < 30) + return 26; + else if(azimuth >= 30 && azimuth < 36) + return 32; + else if(azimuth >= 36 && azimuth < 43) + return 39; + else if(azimuth >= 43 && azimuth < 49) + return 45; + else if(azimuth >= 49 && azimuth < 55) + return 51; + else if(azimuth >= 55 && azimuth < 62) + return 58; + else if(azimuth >= 62 && azimuth < 68) + return 64; + else if(azimuth >= 68 && azimuth < 75) + return 71; + else if(azimuth >= 75 && azimuth < 81) + return 77; + else if(azimuth >= 81 && azimuth < 88) + return 84; + else if(azimuth >= 88 && azimuth < 94) + return 90; + else if(azimuth >= 94 && azimuth < 100) + return 96; + else if(azimuth >= 100 && azimuth < 107) + return 103; + else if(azimuth >= 107 && azimuth < 113) + return 109; + else if(azimuth >= 113 && azimuth < 120) + return 116; + else if(azimuth >= 120 && azimuth < 126) + return 122; + else if(azimuth >= 126 && azimuth < 133) + return 129; + else if(azimuth >= 133 && azimuth < 139) + return 135; + else if(azimuth >= 139 && azimuth < 145) + return 141; + else if(azimuth >= 145 && azimuth < 152) + return 148; + else if(azimuth >= 152 && azimuth < 158) + return 154; + else if(azimuth >= 158 && azimuth < 165) + return 161; + else if(azimuth >= 165 && azimuth < 171) + return 167; + else if(azimuth >= 171 && azimuth < 178) + return 174; + else + return 180; +}; + + + +int mit_hrtf_findIndexFor40Elev(int azimuth) +{ + if(azimuth >= 0 && azimuth < 4) + return 0; + else if(azimuth >= 4 && azimuth < 10) + return 1; + else if(azimuth >= 10 && azimuth < 17) + return 2; + else if(azimuth >= 17 && azimuth < 23) + return 3; + else if(azimuth >= 23 && azimuth < 30) + return 4; + else if(azimuth >= 30 && azimuth < 36) + return 5; + else if(azimuth >= 36 && azimuth < 43) + return 6; + else if(azimuth >= 43 && azimuth < 49) + return 7; + else if(azimuth >= 49 && azimuth < 55) + return 8; + else if(azimuth >= 55 && azimuth < 62) + return 9; + else if(azimuth >= 62 && azimuth < 68) + return 10; + else if(azimuth >= 68 && azimuth < 75) + return 11; + else if(azimuth >= 75 && azimuth < 81) + return 12; + else if(azimuth >= 81 && azimuth < 88) + return 13; + else if(azimuth >= 88 && azimuth < 94) + return 14; + else if(azimuth >= 94 && azimuth < 100) + return 15; + else if(azimuth >= 100 && azimuth < 107) + return 16; + else if(azimuth >= 107 && azimuth < 113) + return 17; + else if(azimuth >= 113 && azimuth < 120) + return 18; + else if(azimuth >= 120 && azimuth < 126) + return 19; + else if(azimuth >= 126 && azimuth < 133) + return 20; + else if(azimuth >= 133 && azimuth < 139) + return 21; + else if(azimuth >= 139 && azimuth < 145) + return 22; + else if(azimuth >= 145 && azimuth < 152) + return 23; + else if(azimuth >= 152 && azimuth < 158) + return 24; + else if(azimuth >= 158 && azimuth < 165) + return 25; + else if(azimuth >= 165 && azimuth < 171) + return 26; + else if(azimuth >= 171 && azimuth < 178) + return 27; + else + return 28; +} \ No newline at end of file
diff -r 000000000000 -r 2076b4d80327 velocity.cpp --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/velocity.cpp Tue Apr 07 21:09:22 2015 +0000 @@ -0,0 +1,54 @@ +#include "../include/velocity.h" +#include <iostream> + +Velocity& Velocity::operator= (const Velocity& rhs) +{ + dx = (rhs.dx); + dy = (rhs.dy); + dz = (rhs.dz); + return *this; +} + +Velocity& Velocity::operator+= (const Velocity& rhs) +{ + dx = dx + rhs.dx; + dy = dy + rhs.dy; + dz = dz + rhs.dz; + return *this; +} + +Velocity operator+ (Velocity lhs, const Velocity& rhs) +{ + lhs += rhs; + return lhs; +} + +bool Velocity::operator< (const Velocity& rhs) const +{ + if (dx == rhs.dx){ + if (dy == rhs.dy){ + return dz < rhs.dz; + }return dy < rhs.dy; + }return dx < rhs.dx; +} + +std::ostream& operator<< (std::ostream& os, const Velocity& obj) +{ + os << obj.getdX() <<","<< obj.getdY() <<","<< obj.getdZ(); + return os; +} + +float Velocity::getdX (void) const +{ + return this->dx; +} + +float Velocity::getdY (void) const +{ + return this->dy; +} + +float Velocity::getdZ (void) const +{ + return this->dz; +}