diff options
author | Thomas Grill <xovo@users.sourceforge.net> | 2002-12-02 19:21:08 +0000 |
---|---|---|
committer | Thomas Grill <xovo@users.sourceforge.net> | 2002-12-02 19:21:08 +0000 |
commit | 9815096db22c73cacdbb65512d1b61d633db7fa8 (patch) | |
tree | 4a6582ead85b8efd031f68e717fbc8a5b3a3df3f /externals/grill/vasp/source/rvfft.cpp | |
parent | 0a109da279e9df66fb5ea7d6bdaeffed16592f02 (diff) |
"version 0.1.1"
svn path=/trunk/; revision=267
Diffstat (limited to 'externals/grill/vasp/source/rvfft.cpp')
-rw-r--r-- | externals/grill/vasp/source/rvfft.cpp | 357 |
1 files changed, 357 insertions, 0 deletions
diff --git a/externals/grill/vasp/source/rvfft.cpp b/externals/grill/vasp/source/rvfft.cpp new file mode 100644 index 00000000..b7d81c93 --- /dev/null +++ b/externals/grill/vasp/source/rvfft.cpp @@ -0,0 +1,357 @@ +#include <math.h> + +#ifdef _MSC_VER +#pragma warning(disable: 4244) +#endif + +#define PI 3.14159265358979 + +///////////////////////////////////////////////////////// +// Sorensen in-place split-radix FFT for real values +// data: array of floats: +// re(0),re(1),re(2),...,re(size-1) +// +// output: +// re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1) +// normalized by array length +// +// Source: +// Sorensen et al: Real-Valued Fast Fourier Transform Algorithms, +// IEEE Trans. ASSP, ASSP-35, No. 6, June 1987 + +void realfft_split(float *data,int n) +{ + + int i,j,k,i5,i6,i7,i8,i0,id,i1,i2,i3,i4,n2,n4,n8; + float t1,t2,t3,t4,t5,t6,a3,ss1,ss3,cc1,cc3,a,e,sqrt2; + + sqrt2=sqrt(2.0); + n4=n-1; + + //data shuffling + for (i=0,j=0,n2=n/2; i<n4 ; i++){ + if (i<j){ + t1=data[j]; + data[j]=data[i]; + data[i]=t1; + } + k=n2; + while (k<=j){ + j-=k; + k>>=1; + } + j+=k; + } + +/*----------------------*/ + + //length two butterflies + i0=0; + id=4; + do{ + for (; i0<n4; i0+=id){ + i1=i0+1; + t1=data[i0]; + data[i0]=t1+data[i1]; + data[i1]=t1-data[i1]; + } + id<<=1; + i0=id-2; + id<<=1; + } while ( i0<n4 ); + + /*----------------------*/ + //L shaped butterflies +n2=2; +for(k=n;k>2;k>>=1){ + n2<<=1; + n4=n2>>2; + n8=n2>>3; + e = 2*PI/(n2); + i1=0; + id=n2<<1; + do{ + for (; i1<n; i1+=id){ + i2=i1+n4; + i3=i2+n4; + i4=i3+n4; + t1=data[i4]+data[i3]; + data[i4]-=data[i3]; + data[i3]=data[i1]-t1; + data[i1]+=t1; + if (n4!=1){ + i0=i1+n8; + i2+=n8; + i3+=n8; + i4+=n8; + t1=(data[i3]+data[i4])/sqrt2; + t2=(data[i3]-data[i4])/sqrt2; + data[i4]=data[i2]-t1; + data[i3]=-data[i2]-t1; + data[i2]=data[i0]-t2; + data[i0]+=t2; + } + } + id<<=1; + i1=id-n2; + id<<=1; + } while ( i1<n ); + a=e; + for (j=2; j<=n8; j++){ + a3=3*a; + cc1=cos(a); + ss1=sin(a); + cc3=cos(a3); + ss3=sin(a3); + a=j*e; + i=0; + id=n2<<1; + do{ + for (; i<n; i+=id){ + i1=i+j-1; + i2=i1+n4; + i3=i2+n4; + i4=i3+n4; + i5=i+n4-j+1; + i6=i5+n4; + i7=i6+n4; + i8=i7+n4; + t1=data[i3]*cc1+data[i7]*ss1; + t2=data[i7]*cc1-data[i3]*ss1; + t3=data[i4]*cc3+data[i8]*ss3; + t4=data[i8]*cc3-data[i4]*ss3; + t5=t1+t3; + t6=t2+t4; + t3=t1-t3; + t4=t2-t4; + t2=data[i6]+t6; + data[i3]=t6-data[i6]; + data[i8]=t2; + t2=data[i2]-t3; + data[i7]=-data[i2]-t3; + data[i4]=t2; + t1=data[i1]+t5; + data[i6]=data[i1]-t5; + data[i1]=t1; + t1=data[i5]+t4; + data[i5]-=t4; + data[i2]=t1; + } + id<<=1; + i=id-n2; + id<<=1; + } while(i<n); + } + } +} + + +///////////////////////////////////////////////////////// +// Sorensen in-place inverse split-radix FFT for real values +// data: array of doubles: +// re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1) +// +// output: +// re(0),re(1),re(2),...,re(size-1) +// NOT normalized by array length +// +// Source: +// Sorensen et al: Real-Valued Fast Fourier Transform Algorithms, +// IEEE Trans. ASSP, ASSP-35, No. 6, June 1987 + +void irealfft_split(float *data,int n){ + + int i,j,k,i5,i6,i7,i8,i0,id,i1,i2,i3,i4,n2,n4,n8,n1; + float t1,t2,t3,t4,t5,a3,ss1,ss3,cc1,cc3,a,e,sqrt2; + + sqrt2=sqrt(2.0); + +n1=n-1; +n2=n<<1; +for(k=n;k>2;k>>=1){ + id=n2; + n2>>=1; + n4=n2>>2; + n8=n2>>3; + e = 2*PI/(n2); + i1=0; + do{ + for (; i1<n; i1+=id){ + i2=i1+n4; + i3=i2+n4; + i4=i3+n4; + t1=data[i1]-data[i3]; + data[i1]+=data[i3]; + data[i2]*=2; + data[i3]=t1-2*data[i4]; + data[i4]=t1+2*data[i4]; + if (n4!=1){ + i0=i1+n8; + i2+=n8; + i3+=n8; + i4+=n8; + t1=(data[i2]-data[i0])/sqrt2; + t2=(data[i4]+data[i3])/sqrt2; + data[i0]+=data[i2]; + data[i2]=data[i4]-data[i3]; + data[i3]=2*(-t2-t1); + data[i4]=2*(-t2+t1); + } + } + id<<=1; + i1=id-n2; + id<<=1; + } while ( i1<n1 ); + a=e; + for (j=2; j<=n8; j++){ + a3=3*a; + cc1=cos(a); + ss1=sin(a); + cc3=cos(a3); + ss3=sin(a3); + a=j*e; + i=0; + id=n2<<1; + do{ + for (; i<n; i+=id){ + i1=i+j-1; + i2=i1+n4; + i3=i2+n4; + i4=i3+n4; + i5=i+n4-j+1; + i6=i5+n4; + i7=i6+n4; + i8=i7+n4; + t1=data[i1]-data[i6]; + data[i1]+=data[i6]; + t2=data[i5]-data[i2]; + data[i5]+=data[i2]; + t3=data[i8]+data[i3]; + data[i6]=data[i8]-data[i3]; + t4=data[i4]+data[i7]; + data[i2]=data[i4]-data[i7]; + t5=t1-t4; + t1+=t4; + t4=t2-t3; + t2+=t3; + data[i3]=t5*cc1+t4*ss1; + data[i7]=-t4*cc1+t5*ss1; + data[i4]=t1*cc3-t2*ss3; + data[i8]=t2*cc3+t1*ss3; + } + id<<=1; + i=id-n2; + id<<=1; + } while(i<n1); + } + } + + /*----------------------*/ + i0=0; + id=4; + do{ + for (; i0<n1; i0+=id){ + i1=i0+1; + t1=data[i0]; + data[i0]=t1+data[i1]; + data[i1]=t1-data[i1]; + } + id<<=1; + i0=id-2; + id<<=1; + } while ( i0<n1 ); + +/*----------------------*/ + +//data shuffling + for (i=0,j=0,n2=n/2; i<n1 ; i++){ + if (i<j){ + t1=data[j]; + data[j]=data[i]; + data[i]=t1; + } + k=n2; + while (k<=j){ + j-=k; + k>>=1; + } + j+=k; + } +} + + +#if 0 +///////////////////////////////////////////////////////// +// Sorensen in-place radix-2 FFT for real values +// data: array of floats: +// re(0),re(1),re(2),...,re(size-1) +// +// output: +// re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1) +// normalized by array length +// +// Source: +// Sorensen et al: Real-Valued Fast Fourier Transform Algorithms, +// IEEE Trans. ASSP, ASSP-35, No. 6, June 1987 + +void realfft_radix2(float *data,int n){ + + float xt,a,e, t1, t2, cc, ss; + int i, j, k, n1, n2, n3, n4, i1, i2, i3, i4; + + n4=n-1; + //data shuffling + for (i=0,j=0,n2=n/2; i<n4 ; i++){ + if (i<j){ + xt=data[j]; + data[j]=data[i]; + data[i]=xt; + } + k=n2; + while (k<=j){ + j-=k; + k>>=1; + } + j+=k; + } + +/* -------------------- */ + for (i=0; i<n; i += 2) + { + xt = data[i]; + data[i] = xt + data[i+1]; + data[i+1] = xt - data[i+1]; + } +/* ------------------------ */ + n2 = 1; + for (k=n;k>2;k>>=1){ + n4 = n2; + n2 = n4 << 1; + n1 = n2 << 1; + e = 2*PI/(n1); + for (i=0; i<n; i+=n1){ + xt = data[i]; + data[i] = xt + data[i+n2]; + data[i+n2] = xt-data[i+n2]; + data[i+n4+n2] = -data[i+n4+n2]; + a = e; + n3=n4-1; + for (j = 1; j <=n3; j++){ + i1 = i+j; + i2 = i - j + n2; + i3 = i1 + n2; + i4 = i - j + n1; + cc = cos(a); + ss = sin(a); + a += e; + t1 = data[i3] * cc + data[i4] * ss; + t2 = data[i3] * ss - data[i4] * cc; + data[i4] = data[i2] - t2; + data[i3] = -data[i2] - t2; + data[i2] = data[i1] - t1; + data[i1] += t1; + } + } + } +} +#endif |