附录一、快速付里哀变换与反变换程序实例
#include <>
#include <>
#define pi (double)3.**********
/*复数定义*/
typedef struct
{
double re;
double im;
}COMPLEX;
/*复数加运算*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
COMPLEX c;
=+;
=+;
return c;
}
/*复数减运算*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
COMPLEX c;
=-;
=-;
return c;
}
/*复数乘运算*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
COMPLEX c;
=*-*;
=*+*;
return c;
}
/*快速付里哀变换
TD为时域值,FD为频域值,power为2的幂数*/
void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{
int count;
int i,j,k,bfsize,p;
double angle;
COMPLEX *W,*X1,*X2,*X;
/*计算付里哀变换点数*/
count=1<<power;
/*分配运算所需存储器*/
W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);
X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
/*计算加权系数*/
for(i=0;i<count/2;i++)
{
angle=-i*pi*2/count;
W[i].re=cos(angle);
W[i].im=sin(angle);
}
/*将时域点写入存储器*/
memcpy(X1,TD,sizeof(COMPLEX)*count);
/*蝶形运算*/
for(k=0;k<power;k++)
{
for(j=0;j<1<<k;j++)
{
bfsize=1<<power-k;
for(i=0;i<bfsize/2;i++)
{
p=j*bfsize;
X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
}
}
X=X1;
X1=X2;
X2=X;
}
/*重新排序*/
for(j=0;j<count;j++)
{
p=0;
for(i=0;i<power;i++)
{
if (j&(1<<i)) p+=1<<power-i-1;
}
FD[j]=X1[p];
}
/*释放存储器*/
free(W);
free(X1);
free(X2);
}
/*快速付里哀反变换,利用快速付里哀变换
FD为频域值,TD为时域值,power为2的幂数*/
void IFFT(COMPLEX *FD, COMPLEX *TD, int power)
{
int i,count;
COMPLEX *x;
/*计算付里哀反变换点数*/
count=1<<power;
/*分配运算所需存储器*/
x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
/*将频域点写入存储器*/
memcpy(x,FD,sizeof(COMPLEX)*count);
/*求频域点的共轭*/
for(i=0;i<count;i++)
{
x[i].im=-x[i].im;
}
/*调用快速付里哀变换*/
FFT(x,TD,power);
/*求时域点的共轭*/
for(i=0;i<count;i++)
{
数字图像处理C语言程序实例 来自淘豆网www.taodocs.com转载请标明出处.