/////////////////离散傅里叶变化//////////
#include"math.h"
void dft(x,y,a,b,n,sign)
int n,sign;
double x[],y[],a[],b[];
{
int i,k;
double c,d,q,w,s;
q=6.28318530718/n;
for(k=0;k<n;k++)
{
w=k*q;
a[k]=b[k]=0.0;
for(i=0;i<n;i++)
{
d=i*w;
c=cos(d);
s=sin(d)*sign;
a[k]+=c*x[i]+s*y[i];
b[k]+=c*y[i]-s*x[i];
}
}
if(sign==-1)
{
c=1.0/n;
for(k=0;k<n;k++)
{a[k]=c*a[k];
b[k]=c*b[k];
}
}
}
/////////////快速傅里叶变化//////
#include"math.h"
void fht4(x,n)
int n;
double x[];
{
int i,j,k,m,n1,n2,n4,l11,l12,l13,l14,l21,l22,l23,l24;
double a,e,c1,c2,c3,s1,s2,s3,t3,t1,t2,t4,t12,t13,t14,t22,t23,t24;
for(j=1,i=1;i<16;i++)
{
m=i;
j=4*j;
if(j==n)break;
}
n1=n-1;
for(j=0,i=0;i<n1;i++)
{
if(i<j)
{
a=x[j];
x[j]=x[i];
x[i]=a;
}
k=n/4;
while(3*k<(j+1))
{
j=j-3*k;
k=k/4;
}
j=j+k;
}
for(i=0;i<n;i+=4)
{
t1=x[i]+x[i+1];
t2=x[i]-x[i+1];
t3=x[i+2]+x[i+3];
t4=x[i+2]-x[i-3];
x[i]=t1+t3;
x[i+1]=t1-t3;
x[i+2]=t2+t4;
x[i+3]=t2-t4;
}
n2=1;
for (k=1;k<m;k++)
{
n4=2*n2;
n2=4*n2;
n1=4*n2;
e=6.283185307179586/n1;
for(i=0;i<n;i+=n1)
{
l12=i+n2;
l13=l12+n2;
l14=l13+n2;
t1=x[i]+x[l12];
t2=x[i]-x[l12];
t3=x[l13]+x[l14];
t4=x[l13]-x[l14];
x[i]=t1+t3;
x[l12]=t1-t3;
x[l13]=t2+t4;
x[l14]=t2-t4;
l21=i+n4;
l22=l21+n2;
l23=l22+n2;
l24=l23+n2;
t1=x[l21];
t2=x[l22]*sqrt(2.0);
t3=x[l23];
t4=x[l24]*sqrt(2.0);
x[l21]=t1+t2+t3;
x[l22]=t1-t3+t4;
x[l23]=t1-t2+t3;
x[l24]=t1-t3-t4;
a=e;
for(j=1;j<n4;j++)
{
l11=i+j;
l12=l11+n2;
l13=l12+n2;
l14=l13+n2;
l21=i+n2-j;
l22=l21+n2;
l23=l22+n2;
l24=l23+n2;
c1=cos(a);
s1=sin(a);
c2=cos(2*a);
s2=sin(2*a);
c3=cos(3*a);
s3=sin(3*a);
a=(j+1)*e;
t12=x[l12]*c1+x[l22]*s1;
t13=x[l13]*c2+x[l23]*s2;
t14=x[l14]*c3+x[l24]*s3;
t22=x[l12]*s1-x[l22]*c1;
t23=x[l13]*s2-x[l23]*c2;
t24=x[l14]*s3-x[l24]*c3;
t1=x[l21]+t23;
t2=x[l21]-t23;
t3=t22+t24;
t4=t12-t14;
x[l24]=t2-t3;
x[l23]=t1-t4;
x[l22]=t2+t3;
x[l21]=t1+t4;
t1=x[l11]+t13;
t2=x[l11]-t13;
t3=t24-t22;
t4=t12+t14;
x[l14]=t2-t3;
x[l13]=t1-t4;
x[l12]=t2+t3;
x[l11]=t1+t4;
}
}
}
}
///////希尔伯特变化/////////
#include"math.h"
#include"stdio.h"
#include "hilbth.c"
main()
{
int i,n;
double x[64];
FILE *fp;
n=64;
for(i=0;i<n;i++)
x[i]=sin(2*3.14159265*i/32);
fp=fopen("hilb1.dat","w");
for(i=0;i<n;i++)
fprintf(fp,"%d %f\n",i,x[i]);
fclose(fp);
printf("\n oringal sequence\n");
for(i=0;i<n/2;i+=4)
{
printf("%10.7f%10.7f",x[i],x[i+1]);
printf("%10.7f%10.7f",x[i+2],x[i+3]);
}
hilbth(x,n);
printf("\n discrete hiblt trasform\n");
for(i=0;i<n/2;i+=4)
{
printf("%10.7f%10.7",x[i],x[i+1]);
printf("%10.7f%10.7",x[i+2],x[i+3]);
}
fp=fopen("hilb2.dat","w");
for(i=0;i<n;i++)
fprintf(fp,"%d %f\n",i,x[i]);
fclose(fp);
}
