DQFT.cpp0100644000076400010010000006556207661116645012645 0ustar AdministratorNone#include "stdafx.h" //#include //#include //#include "basic.hpp" #include "qft.h" DQFT dqft; #define sec1 0.5411961001461970 #define sec2 0.70710678118654752 #define sec3 1.3065629648763763 static const double Sec1=sec1; static const double Sec2=sec2; static const double Sec3=sec3; DSecantTable DQFT::Secants(MaxQFT); #define _USE_ASM void DQFT::Dct(double *x) { if (lnh==2) { #ifndef _USE_ASM const double t1=x[0]-x[4]; const double t2=0.70710678118654*(x[1]-x[3]); const double t3=x[1]+x[3]; const double t4=x[0]+x[4]; const double t5=t4+x[2]; x[0]=(double)(t5+t3); x[1]=(double)(t1+t2); x[2]=(double)(t4-x[2]); x[3]=(double)(t1-t2); x[4]=(double)(t5-t3); #else _asm{ mov eax,x; #define X(a) qword ptr [eax + (a*8)] // {0 1 2 3 4 5 6 7 8} fld X(0); //ld r1,X(0) { r1} fld X(1); //ld r2,X(1) { r2, r1} fld st(1); //ld r4,r1 { r4, r2, r1} fxch st(2); // { r1, r2, r4} fsub X(4); //r1-X(4) { r1, r2, r4} fxch st(1); // { r2, r1, r4} fsub X(3); //r2-X(3) { r2, r1, r4} fxch st(2); // { r4, r1, r2} fadd X(4); //r4+X(4) { r4, r1, r2} fld X(1); //ld r3,X(1) { r3, r4, r1, r2} fxch st(3); // { r2, r4, r1, r3} fmul Sec2; //r2*1/sqrt2 { r2, r4, r1, r3} fld X(2); //ld r5,X(2) { r5, r2, r4, r1, r3} fxch st(4); // { r3, r2, r4, r1, r5} fadd X(3); //r3+X(3) { r3, r2, r4, r1, r5} fxch st(4); // { r5, r2, r4, r1, r3} fadd st(0),st(2); //r5+r4 { r5, r2, r4, r1, r3} fxch st(2); // { r4, r2, r5, r1, r3} fsub X(2); //r4-X[2] { r4, r2, r5, r1, r3} fld st(1); //ld r6,r2 { r6, r4, r2, r5, r1, r3} fld st(5); //ld r7,r3 { r7, r6, r4, r2, r5, r1, r3} fxch st(1); // { r6, r7, r4, r2, r5, r1, r3} fadd st(0),st(5); //r6+r1 { r6, r7, r4, r2, r5, r1, r3} fxch st(2); // { r4, r7, r6, r2, r5, r1, r3} fstp X(2); //from r4 { r7, r6, r2, r5, r1, r3} fadd st(0),st(3); //r7+r5 { r7, r6, r2, r5, r1, r3} fxch st(2); // { r2, r6, r7, r5, r1, r3} fsubp st(4),st(0); //r1-r2 { r6, r7, r5, r1, r3} fxch st(4); // { r3, r7, r5, r1, r6} fsubp st(2),st(0); //r5-r3 { r7, r5, r1, r6} fxch st(3); // { r6, r5, r1, r7} fstp X(1); //from r6 { r5, r1, r7} fxch st(2); // { r7, r1, r5} fstp X(0); //from r7 { r1, r5} fstp X(3); //from r1 { r5} fstp X(4); //from r5 {} } #undef X #endif return; } { const int l=(1<>1; const int hp1=h+1; double *secants = Secants[lnh]; double *x1=currentScratch; double *x2=x1+hp1; int i; currentScratch=x2+hp1; #ifndef _USE_ASM for (i=0;i>1; const int hm1=h-1; double *x1=currentScratch; double *x2=x1+hm1; int i; double *secants = Secants[lnh]+1; currentScratch=x2+hm1; #ifndef _USE_ASM for (i=0;i>1; R[0]=x[0]; R[h]=x[h]; #ifndef _USE_ASM for (int i=1;i>1]=0.0; } void DQFT::Forward(double *R, double *I, float *x, int _lnh) { // Profile(QFT_Foward); int h; // int i; int n=1<<_lnh; CriticalSectionEntry ce(cs); assert(_lnh<=maxln); lnh = _lnh-1; h=n>>1; R[0]=x[0]; R[h]=x[h]; #ifndef _USE_ASM for (int i=1;i>1]=0.0; } void DQFT::ReverseNoUnpack(double *R, double *I, int _lnh) { // Profile(QFT_Reverse); int h; int n=1<<_lnh; CriticalSectionEntry ce(cs); assert(_lnh<=maxln); lnh = _lnh-1; h=n>>1; R[0]*=0.5; R[h]*=0.5; Dct(R); Dst(I+1); I[0]=0.0f; } void DQFT::Reverse(double *x, double *R, double *I, int _lnh) { // Profile(QFT_Reverse); int h; int n=1<<_lnh; CriticalSectionEntry ce(cs); assert(_lnh<=maxln); lnh = _lnh-1; h=n>>1; R[0]*=0.5; R[h]*=0.5; Dct(R); Dst(I+1); const double nr=2.0/n; x[0]=R[0]*nr; #ifndef _USE_ASM for (int i=1;i>1; R[0]*=0.5; R[h]*=0.5; Dct(R); Dst(I+1); const double nr=2.0/n; x[0]=float(R[0]*nr); #ifndef _USE_ASM for (int i=1;i p(cos(-2*PI/n),sin(-2*PI/n)); vector fe(hh+1,0.0); vector fo(hh+1,0.0); fe[0]=(f[0]+f[h])*0.5; fo[0]=(f[0]-f[h])*0.5; int j; for (j=1;j<=hh;++j) { fe[j]=(f[j]+conjugate(f[h-j]))*0.5; fo[j]=(f[j]-conjugate(f[h-j]))*p^(-j)*0.5; }; */ double *rotations = Secants.Rotations(_lnh)+2; const long n = 1<<_lnh; const int h=n>>1; const int hh=h>>1; re[0]=(R[0]+R[h])*0.5; ro[0]=(R[0]-R[h])*0.5; ie[0]= io[0]= 0.0; for (int j=1;j r(h+1,0.0); complex p(cos(-2*PI/n),sin(-2*PI/n)); r[0]=fe[0]+fo[0]; r[h]=fe[0]-fo[0]; for(j=1;j>1; const int hh=h>>1; double sr=*re - *ro; *R=re[0]+ro[0]; *I=0.0; for(int i=1;i12){ const long n = 1<>1; const int hh=h>>1; const int ansLen = hh+1; int i; double *rotations = Secants.Rotations(lnh); double *fo = currentScratch; double *er = fo+h; double *ei = er+ansLen; double *fe = ei+ansLen; double *or, *oi; currentScratch = fe+h; for(i=0;i #include #include #include "qft.h" using namespace std; union int64Union { struct { int low; int high; } ints; __int64 int64; }; inline __int64 GetTimeStamp() { int64Union retval; __asm{ rdtsc; mov retval.ints.low,eax; mov retval.ints.high,edx; } return retval.int64; } void testsingle(int fftlenln2) { const int fftlen = 1< destI(freqlen); vector destR(freqlen); vector source(fftlen); //if we leave zeros in we'd be testing the slowness of denormalized //numbers... So we fill with noise int i,j; srand(100); for (i=0;i destI(freqlen); vector destR(freqlen); vector source(fftlen); //if we leave zeros in we'd be testing the slowness of denormalized //numbers... So we fill with noise int i,j; srand(100); for (i=0;i //#include //#include "basic.hpp" #include "qft.h" QFT qft; #define sec1 0.5411961001461970 #define sec2 0.70710678118654752 #define sec3 1.3065629648763763 static const double Sec1=sec1; static const double Sec2=sec2; static const double Sec3=sec3; SecantTable QFT::Secants(MaxQFT); #define _USE_ASM void QFT::Dct(float *x) { if (lnh==2) { #ifndef _USE_ASM const double t1=x[0]-x[4]; const double t2=0.70710678118654*(x[1]-x[3]); const double t3=x[1]+x[3]; const double t4=x[0]+x[4]; const double t5=t4+x[2]; x[0]=(float)(t5+t3); x[1]=(float)(t1+t2); x[2]=(float)(t4-x[2]); x[3]=(float)(t1-t2); x[4]=(float)(t5-t3); #else _asm{ mov eax,x; #define X(a) dword ptr [eax + (a*4)] // {0 1 2 3 4 5 6 7 8} fld X(0); //ld r1,X(0) { r1} fld X(1); //ld r2,X(1) { r2, r1} fld st(1); //ld r4,r1 { r4, r2, r1} fxch st(2); // { r1, r2, r4} fsub X(4); //r1-X(4) { r1, r2, r4} fxch st(1); // { r2, r1, r4} fsub X(3); //r2-X(3) { r2, r1, r4} fxch st(2); // { r4, r1, r2} fadd X(4); //r4+X(4) { r4, r1, r2} fld X(1); //ld r3,X(1) { r3, r4, r1, r2} fxch st(3); // { r2, r4, r1, r3} fmul Sec2; //r2*1/sqrt2 { r2, r4, r1, r3} fld X(2); //ld r5,X(2) { r5, r2, r4, r1, r3} fxch st(4); // { r3, r2, r4, r1, r5} fadd X(3); //r3+X(3) { r3, r2, r4, r1, r5} fxch st(4); // { r5, r2, r4, r1, r3} fadd st(0),st(2); //r5+r4 { r5, r2, r4, r1, r3} fxch st(2); // { r4, r2, r5, r1, r3} fsub X(2); //r4-X[2] { r4, r2, r5, r1, r3} fld st(1); //ld r6,r2 { r6, r4, r2, r5, r1, r3} fld st(5); //ld r7,r3 { r7, r6, r4, r2, r5, r1, r3} fxch st(1); // { r6, r7, r4, r2, r5, r1, r3} fadd st(0),st(5); //r6+r1 { r6, r7, r4, r2, r5, r1, r3} fxch st(2); // { r4, r7, r6, r2, r5, r1, r3} fstp X(2); //from r4 { r7, r6, r2, r5, r1, r3} fadd st(0),st(3); //r7+r5 { r7, r6, r2, r5, r1, r3} fxch st(2); // { r2, r6, r7, r5, r1, r3} fsubp st(4),st(0); //r1-r2 { r6, r7, r5, r1, r3} fxch st(4); // { r3, r7, r5, r1, r6} fsubp st(2),st(0); //r5-r3 { r7, r5, r1, r6} fxch st(3); // { r6, r5, r1, r7} fstp X(1); //from r6 { r5, r1, r7} fxch st(2); // { r7, r1, r5} fstp X(0); //from r7 { r1, r5} fstp X(3); //from r1 { r5} fstp X(4); //from r5 {} } #undef X #endif return; } { const int l=(1<>1; const int hp1=h+1; float *secants = Secants[lnh]; float *x1=currentScratch; float *x2=x1+hp1; currentScratch=x2+hp1; // */ #ifndef _USE_ASM int i; for (i=0;i>1; const int hm1=h-1; float *x1=currentScratch; float *x2=x1+hm1; float *secants = Secants[lnh]+1; currentScratch=x2+hm1; #ifndef _USE_ASM int i; for (i=0;i>1; R[0]=x[0]; R[h]=x[h]; #ifndef _USE_ASM int i; for (i=1;i>1]=0.0; } void QFT::ReverseNoUnpack(float *R, float *I, int _lnh) { // Profile(QFT_Reverse); int h; int n=1<<_lnh; CriticalSectionEntry ce(cs); assert(_lnh<=maxln); lnh = _lnh-1; h=n>>1; R[0]*=0.5; R[h]*=0.5; Dct(R); Dst(I+1); I[0]=0.0f; } void QFT::ReverseHalf(float *x, float *R, float *I, int _lnh) { // Profile(QFT_Reverse); int h; int n=1<<_lnh; CriticalSectionEntry ce(cs); assert(_lnh<=maxln); lnh = _lnh-1; h=n>>1; R[0]*=0.5; R[h]*=0.5; Dct(R); Dst(I+1); // const double nr=2.0/n; We assume that the normalization was in the convolution already // multiply the cooeficients by 2.0/_lnh beforehand x[0]=float(R[0]); for (int i=1;i>1; R[0]*=0.5; R[h]*=0.5; Dct(R); Dst(I+1); const double nr=2.0/n; x[0]=float(R[0]*nr); #ifndef _USE_ASM for (int i=1;i p(cos(-2*PI/n),sin(-2*PI/n)); vector fe(hh+1,0.0); vector fo(hh+1,0.0); fe[0]=(f[0]+f[h])*0.5; fo[0]=(f[0]-f[h])*0.5; int j; for (j=1;j<=hh;++j) { fe[j]=(f[j]+conjugate(f[h-j]))*0.5; fo[j]=(f[j]-conjugate(f[h-j]))*p^(-j)*0.5; }; */ float *rotations = Secants.Rotations(_lnh)+2; const long n = 1<<_lnh; const int h=n>>1; const int hh=h>>1; re[0]=float((R[0]+R[h])*0.5); ro[0]=float((R[0]-R[h])*0.5); ie[0]= io[0]= 0.0; for (int j=1;j r(h+1,0.0); complex p(cos(-2*PI/n),sin(-2*PI/n)); r[0]=fe[0]+fo[0]; r[h]=fe[0]-fo[0]; for(j=1;j>1; const int hh=h>>1; float sr=*re - *ro; *R=re[0]+ro[0]; *I=0.0; for(int i=1;i12){ const long n = 1<>1; const int hh=h>>1; const int ansLen = hh+1; int i; float *rotations = Secants.Rotations(lnh); float *fo = currentScratch; float *er = fo+h; float *ei = er+ansLen; float *fe = ei+ansLen; float *or, *oi; currentScratch = fe+h; #ifndef _USE_ASM for(i=0;i #include //#include "basic.hpp" #include #include "qft.h" using namespace std; SecantTable::~SecantTable() { for (int i=2;i<=maxln;++i){ delete [] secants[i]; delete [] rotations[i]; } delete [] secants; delete [] rotations; } SecantTable::SecantTable(int _maxln) { int i,j; const double pi = 3.141592653589793238462643383; maxln = _maxln; secants=new float *[maxln+1]; rotations=new float *[maxln+1]; for (i=2;i<=maxln;++i){ complex a(2.0,0.0); double theta=pi/(1< p(cos(theta),sin(theta)); int entries = (1<<(i-1)); float *table= new float [entries]; secants[i]=table; for (j=0;j(1.0,0.0); p= complex(cos(theta),sin(theta)); for (j=0;j<=entries;++j) { *table++ = a.real(); *table++ = a.imag(); a*=p; } } } DSecantTable::~DSecantTable() { for (int i=2;i<=maxln;++i){ delete [] secants[i]; delete [] rotations[i]; } delete [] secants; delete [] rotations; } //Note these tables are built using complex multiplication to //rotate the previous entry. It would be more accurate to //use sin and cos for every entry. DSecantTable::DSecantTable(int _maxln) { int i,j; const double pi = 3.141592653589793238462643383; maxln = _maxln; secants=new double *[maxln+1]; rotations=new double *[maxln+1]; for (i=2;i<=maxln;++i){ complex a(2.0,0.0); double theta=pi/(1< p(cos(theta),sin(theta)); int entries = (1<<(i-1)); double *table= new double [entries]; secants[i]=table; for (j=0;j(1.0,0.0); p= complex(cos(theta),sin(theta)); for (j=0;j<=entries;++j) { *table++ = a.real(); *table++ = a.imag(); a*=p; } } } CriticalSection.h0100644000076400010010000000210106672323672014610 0ustar AdministratorNone#ifndef INCLUDED_CRITICALSECTION_H #define INCLUDED_CRITICALSECTION_H #ifndef _WINDOWS_ #include #endif class CriticalSectionEntry; class CriticalSection { public: CriticalSection() { InitializeCriticalSection(&cs); } ~CriticalSection() { DeleteCriticalSection(&cs); } protected: friend class CriticalSectionEntry; void Enter() { EnterCriticalSection(&cs); } void Leave() { LeaveCriticalSection(&cs); } CRITICAL_SECTION cs; private: CriticalSection(const CriticalSection&){} void operator =(const CriticalSection&){} }; class AdvancedCriticalSection :public CriticalSection { public: //just changing the access to Enter and Leave: using CriticalSection::Enter; using CriticalSection::Leave; }; class CriticalSectionEntry { public: explicit CriticalSectionEntry(const CriticalSection& c) : cs(const_cast(&c)) { cs->Enter(); } ~CriticalSectionEntry() { cs->Leave(); } private: CriticalSection *cs; }; #endif StdAfx.h0100644000076400010010000000123307660353650012724 0ustar AdministratorNone// stdafx.h : include file for standard system include files, // or project specific include files that are used frequently, but // are changed infrequently // #if !defined(AFX_STDAFX_H__C9DE0B39_8BC9_4C7C_95A3_C4E91E0E3CC2__INCLUDED_) #define AFX_STDAFX_H__C9DE0B39_8BC9_4C7C_95A3_C4E91E0E3CC2__INCLUDED_ #if _MSC_VER > 1000 #pragma once #endif // _MSC_VER > 1000 // TODO: reference additional headers your program requires here //{{AFX_INSERT_LOCATION}} // Microsoft Visual C++ will insert additional declarations immediately before the previous line. #endif // !defined(AFX_STDAFX_H__C9DE0B39_8BC9_4C7C_95A3_C4E91E0E3CC2__INCLUDED_) qft.h0100644000076400010010000003374707661120050012327 0ustar AdministratorNone#ifndef _QFT_H #define _QFT_H #include #include "criticalsection.h" //This class initializes the tables used by the //QFT class and the DQFT class class SecantTable { float **secants; float **rotations; int maxln; public: float *operator [](int i) { return secants[i]; } float *Rotations(int i) { return rotations[i]; } SecantTable(int i); ~SecantTable(); }; //Since the QFT routines are not entirely in-place //they need to preallocate scratch memory. // //MaxQFT sets the size of the scratch pad // //maximum block length (2^17) const int MaxQFT = 17; //Other limitations: //The DCT routine can't handle a block size less than 4 //The DST routine can't handle a block size less than 8 //The FFT routines (both forward and reverse) can't handle //block sizes less than 16 //Block sizes must be powers of two //Warning there's some weirdness where at some block sizes //(128 & 256 I think) //allocating the data and real and imaginary blocks //consecutively slows the program down by 10000% //??? I don't know why //Usually I cludge around this by using STL or new to //allocate data. //Note that since scratch pad memory is set aside //it would not ordinarily be safe to use QFT from //multiple threads simultaniously. //There are two provisions made for this: //1. you can create more than one QFT object // and use different ones in different threads //2. As currently written, a QFT object has a critical // section that makes it block when called called // simultanious from different threads. //Alternately you might consider making the scratch //area thread local. /* The QFT and DQFT (double precision) classes supply the following functions: 1. Real valued FFT and inverse FFT functions. Note that separate arrays are used for real and imaginary component of the resulting spectrum. 2. Decomposition of a spectrum into a separate spectrum of the even samples and a spectrum of the odd samples. This can be useful for building filter banks. 3. Reconstituting a spectrum from separate spectrums of the even samples and odd samples. This can be useful for building filter banks. 4. A discrete Sin transform (a QFT decomposes an FFT into a DST and DCT). 5. A discrete Cos transfrom. 6. Since a QFT does it's last stage calculating from the outside in the last part can be left unpacked and only calculated as needed in the case where the entire spectrum isn't needed (I used this for calculating correlations and convolutions where I only needed half of the results). ReverseNoUnpack() UnpackStep() and NegUnpackStep() implement this functionality NOTE Reverse() normalizes its results (divides by one half the block length), but ReverseNoUnpack() does not. 7. Also if you only want the first half of the results you can call ReverseHalf() NOTE Reverse() normalizes its results (divides by one half the block length), but ReverseHalf() does not. 8. QFT is less numerically stable than regular FFTs. With single precision calculations, a block length of 2^15 brings the accuracy down to being barely accurate enough. At that size, single precision calculations tested sound files would occasionally have a sample off by 2, and a couple off by 1 per block. Full volume white noise would generate a few samples off by as much as 6 per block at the end, beginning and middle. No matter what the inputs the errors are always at the same positions in the block. There some sort of cancelation that gets more delicate as the block size gets bigger. For the sake of doing convolutions and the like where the forward transform is done only once for one of the inputs, I created a AccurateForward() function. It uses a regular FFT algorithm for blocks larger than 2^12, and decomposes into even and odd FFTs recursively. In any case you can always use the double precision routines to get more accuracy. DQFT even has routines that take floats as inputs and return double precision spectrum outputs. */ class QFT { static SecantTable Secants; float *scratch; float *currentScratch; int lnh; int maxln; //Since QFT uses a permanently allocated scratch space //you can only be running one QFT per thread at a time. //The reason I made this a class is so you CAN run //simultanious transforms (by making multiple QFT objects). CriticalSection cs; void Dct(float *x); void Dst(float *x); void jhfftu(float *R, float *I, float *x); public: QFT(int _max=MaxQFT) :scratch(new float[(1<<(_max+3))+2*(_max+1)]) ,currentScratch(scratch) ,maxln(_max) { assert(maxln<=MaxQFT); } ~QFT() { delete [] scratch; } //Discrete Sine transform. _lnh is the base two log of the block length. //Note that block sizes less than 8 (lnh=3) don't work. void Dst(float *x,int _lnh) { CriticalSectionEntry ce(cs); lnh = _lnh; Dst(x); } //Discrete Cosine transform. _lnh is the base two log of the block length. //Note that block sizes less than 4 (lnh=2) don't work. void Dct(float *x,int _lnh) { CriticalSectionEntry ce(cs); lnh = _lnh; x[0]*=.5; x[1< void ForwardFromTemplate(float *R, float *I, T& x, int pos, int _lnh) { int h; int i; int n=1<<_lnh; CriticalSectionEntry ce(cs); assert(_lnh<=maxln); h=n>>1; R[0]=x[pos]; R[h]=x[pos+h]; for (i=1;i void ForwardFromTemplateSkipOdd(float *R, float *I, T& x, int pos, int _lnh) { int h; int i; int n=pos+(2<<_lnh); CriticalSectionEntry ce(cs); assert(_lnh<=maxln); h=1<<(_lnh-1); R[0]=x[pos]; R[h]=x[pos+h+h]; for (i=1;i12 void AccurateForward(float *R, float *I, float *x, int _lnh) { // Profile(QFT_AccurateForward); CriticalSectionEntry ce(cs); lnh=_lnh; jhfftu(R, I, x); } //note, the reverse transform changes R and I in place //make copies if you want to preserve them void Reverse(float *x, float *R, float *I, int _lnh); //for convolutions there's some edges you don't need //this is a special case where it only calculates the //first have of the reverse transform void ReverseHalf(float *x, float *R, float *I, int _lnh); void QFT::ReverseNoUnpack(float *R, float *I, int _lnh); //this only works for the first half static double UnpackStep(int pos, float *R, float *I) { return R[pos]-I[pos]; } //gets you the last half of the spectrum, indexed backwards from the end static double NegUnpackStep(int pos, float *R, float *I) { return R[pos]+I[pos]; } //turns a spectrum into a spectrum of even samples and a spectrum of odd samples static void DecomposeEvenOdd(float *re,float *ie,float *ro,float *io, float *R, float *I, int _lnh); //combines a spectrum of even samples with a spectrum of odd samples static void CombineEvenOdd(float *R, float *I, float *re,float *ie,float *ro,float *io, int _lnh); }; //Note these tables are built using complex multiplication to //rotate the previous entry. It would be more accurate to //use sin and cos for every entry. class DSecantTable { double **secants; double **rotations; int maxln; public: double *operator [](int i) { return secants[i]; } double *Rotations(int i) { return rotations[i]; } DSecantTable(int i); ~DSecantTable(); }; //maximum block length (2^17) const int MaxDQFT = 17; class DQFT { static DSecantTable Secants; double *scratch; double *currentScratch; int lnh; int maxln; //Since DQFT uses a permanently allocated scratch space //you can only be running one DQFT per thread at a time. //The reason I made this a class is so you CAN run //simultanious transforms (by making multiple DQFT objects). CriticalSection cs; void Dct(double *x); void Dst(double *x); void jhfftu(double *R, double *I, double *x); public: DQFT(int _max=MaxDQFT) :scratch(new double[(1<<(_max+3))+2*(_max+1)+(1<<_max)]) ,currentScratch(scratch) ,maxln(_max) { assert(maxln<=MaxQFT); } ~DQFT() { delete [] scratch; } //Discrete Sine transform. _lnh is the base two log of the block length. //Note that block sizes less than 8 (lnh=3) don't work. void Dst(double *x,int _lnh) { CriticalSectionEntry ce(cs); lnh = _lnh; Dst(x); } //Discrete Cosine transform. _lnh is the base two log of the block length. //Note that block sizes less than 4 (lnh=2) don't work. void Dct(double *x,int _lnh) { CriticalSectionEntry ce(cs); lnh = _lnh; x[0]*=.5; x[1< void ForwardFromTemplate(double *R, double *I, T& x, int pos, int _lnh) { int h; int i; int n=1<<_lnh; CriticalSectionEntry ce(cs); assert(_lnh<=maxln); h=n>>1; R[0]=x[pos]; R[h]=x[pos+h]; for (i=1;i void ForwardFromTemplateSkipOdd(double *R, double *I, T& x, int pos, int _lnh) { int h; int i; int n=pos+(2<<_lnh); CriticalSectionEntry ce(cs); assert(_lnh<=maxln); h=1<<(_lnh-1); R[0]=x[pos]; R[h]=x[pos+h+h]; for (i=1;i12 void AccurateForward(double *R, double *I, double *x, int _lnh) { // Profile(QFT_AccurateForward); CriticalSectionEntry ce(cs); lnh=_lnh; jhfftu(R, I, x); } //Starts with an array of floats instead of doubles, but does all //of it's calculations to double accuracy and puts the results in //arrays of doubles void AccurateForward(double *R, double *I, float *x, int _lnh) { int n=1<<_lnh; double *temp = currentScratch; currentScratch += n; for (int i=0;i Package=<5> {{{ }}} Package=<4> {{{ }}} ############################################################################### Global: Package=<5> {{{ }}} Package=<3> {{{ }}} ############################################################################### qft.dsp0100644000076400010010000001122407661120073012655 0ustar AdministratorNone# Microsoft Developer Studio Project File - Name="qft" - Package Owner=<4> # Microsoft Developer Studio Generated Build File, Format Version 6.00 # ** DO NOT EDIT ** # TARGTYPE "Win32 (x86) Console Application" 0x0103 CFG=qft - Win32 Debug !MESSAGE This is not a valid makefile. To build this project using NMAKE, !MESSAGE use the Export Makefile command and run !MESSAGE !MESSAGE NMAKE /f "qft.mak". !MESSAGE !MESSAGE You can specify a configuration when running NMAKE !MESSAGE by defining the macro CFG on the command line. For example: !MESSAGE !MESSAGE NMAKE /f "qft.mak" CFG="qft - Win32 Debug" !MESSAGE !MESSAGE Possible choices for configuration are: !MESSAGE !MESSAGE "qft - Win32 Release" (based on "Win32 (x86) Console Application") !MESSAGE "qft - Win32 Debug" (based on "Win32 (x86) Console Application") !MESSAGE # Begin Project # PROP AllowPerConfigDependencies 0 # PROP Scc_ProjName "" # PROP Scc_LocalPath "" CPP=cl.exe RSC=rc.exe !IF "$(CFG)" == "qft - Win32 Release" # PROP BASE Use_MFC 0 # PROP BASE Use_Debug_Libraries 0 # PROP BASE Output_Dir "Release" # PROP BASE Intermediate_Dir "Release" # PROP BASE Target_Dir "" # PROP Use_MFC 0 # PROP Use_Debug_Libraries 0 # PROP Output_Dir "Release" # PROP Intermediate_Dir "Release" # PROP Target_Dir "" # ADD BASE CPP /nologo /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /Yu"stdafx.h" /FD /c # ADD CPP /nologo /G6 /Gd /MT /W3 /GX /O2 /D "WIN32" /D "NDEBUG" /D "_CONSOLE" /D "_MBCS" /Yu"stdafx.h" /FD /c # ADD BASE RSC /l 0x409 /d "NDEBUG" # ADD RSC /l 0x409 /d "NDEBUG" BSC32=bscmake.exe # ADD BASE BSC32 /nologo # ADD BSC32 /nologo LINK32=link.exe # ADD BASE LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /machine:I386 # ADD LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /machine:I386 !ELSEIF "$(CFG)" == "qft - Win32 Debug" # PROP BASE Use_MFC 0 # PROP BASE Use_Debug_Libraries 1 # PROP BASE Output_Dir "Debug" # PROP BASE Intermediate_Dir "Debug" # PROP BASE Target_Dir "" # PROP Use_MFC 0 # PROP Use_Debug_Libraries 1 # PROP Output_Dir "Debug" # PROP Intermediate_Dir "Debug" # PROP Target_Dir "" # ADD BASE CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /Yu"stdafx.h" /FD /GZ /c # ADD CPP /nologo /W3 /Gm /GX /ZI /Od /D "WIN32" /D "_DEBUG" /D "_CONSOLE" /D "_MBCS" /Yu"stdafx.h" /FD /GZ /c # ADD BASE RSC /l 0x409 /d "_DEBUG" # ADD RSC /l 0x409 /d "_DEBUG" BSC32=bscmake.exe # ADD BASE BSC32 /nologo # ADD BSC32 /nologo LINK32=link.exe # ADD BASE LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept # ADD LINK32 kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib /nologo /subsystem:console /debug /machine:I386 /pdbtype:sept !ENDIF # Begin Target # Name "qft - Win32 Release" # Name "qft - Win32 Debug" # Begin Group "Source Files" # PROP Default_Filter "cpp;c;cxx;rc;def;r;odl;idl;hpj;bat" # Begin Source File SOURCE=.\DQFT.cpp # End Source File # Begin Source File SOURCE=.\main.cpp # End Source File # Begin Source File SOURCE=.\qft.cpp # End Source File # Begin Source File SOURCE=.\secant.cpp # End Source File # Begin Source File SOURCE=.\StdAfx.cpp # ADD CPP /Yc"stdafx.h" # End Source File # End Group # Begin Group "Header Files" # PROP Default_Filter "h;hpp;hxx;hm;inl" # Begin Source File SOURCE=.\CriticalSection.h # End Source File # Begin Source File SOURCE=.\qft.h # End Source File # Begin Source File SOURCE=.\StdAfx.h # End Source File # End Group # Begin Group "Resource Files" # PROP Default_Filter "ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe" # End Group # End Target # End Project