aboutsummaryrefslogtreecommitdiff
path: root/externals/grill/vasp/source/rvfft.cpp
diff options
context:
space:
mode:
authorThomas Grill <xovo@users.sourceforge.net>2009-04-01 21:13:09 +0000
committerThomas Grill <xovo@users.sourceforge.net>2009-04-01 21:13:09 +0000
commit0ed7a8b68dd73e2b0473b8127aeca99f3bac9061 (patch)
tree5c67818b38a5cc2f9caa5ca7f8640ca356adf02b /externals/grill/vasp/source/rvfft.cpp
parentbb4c7f6a245394d09dac9adfb2efb093d3d98452 (diff)
cleaned up grill externals - replaced with svn:externals to svn.grrrr.org/ext/trunk/
svn path=/trunk/; revision=10951
Diffstat (limited to 'externals/grill/vasp/source/rvfft.cpp')
-rw-r--r--externals/grill/vasp/source/rvfft.cpp357
1 files changed, 0 insertions, 357 deletions
diff --git a/externals/grill/vasp/source/rvfft.cpp b/externals/grill/vasp/source/rvfft.cpp
deleted file mode 100644
index b7d81c93..00000000
--- a/externals/grill/vasp/source/rvfft.cpp
+++ /dev/null
@@ -1,357 +0,0 @@
-#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