《数字图像处理技术》课程教学资源（实例）实验 - C语言附录实例

附录快速付里哀变换与反变换程序实例目 #include <math.h> #include <malloc.h> #define pi (double)3.14159265359 /*复数定义*/ typedef struct doublere double im; COMPLEX: /朱复数加运算/ COMPLEX Add(COMPLEX c1,COMPLEX c2) COMPLEX c: c.re=cl.re+c2.re: c.im-c1.im+c2.im; return c: /朱复数减运算*/ COMPLEX Sub (COMPLEX c1,COMPLEX c2) COMPLEX c: c.re=cl.re-c2.re; c.im=c1.im-c2.im: return c: /*复数乘运算* COMPLEX Mul(COMPLEX c1,COMPLEX c2) COMPLEX c: return c: /余快速付里哀变换

附录一、快速付里哀变换与反变换程序实例 #include <math.h> #include <malloc.h> #define pi (double)3.14159265359 /*复数定义*/ typedef struct { double re; double im; }COMPLEX; /*复数加运算*/ COMPLEX Add(COMPLEX c1, COMPLEX c2) { COMPLEX c; c.re=c1.re+c2.re; c.im=c1.im+c2.im; return c; } /*复数减运算*/ COMPLEX Sub(COMPLEX c1, COMPLEX c2) { COMPLEX c; c.re=c1.re-c2.re; c.im=c1.im-c2.im; return c; } /*复数乘运算*/ COMPLEX Mul(COMPLEX c1, COMPLEX c2) { COMPLEX c; c.re=c1.re*c2.re-c1.im*c2.im; c.im=c1.re*c2.im+c2.re*c1.im; return c; } /*快速付里哀变换

TD为时域值，FD为频域值，power为2的幂数*/ void FFT(COMPLEX*TD,COMPLEX FD,int power) int count; int i,j,k,bfsize,p: double angle: COMPLEXW*X1,*X2,*X /计算付里哀变换点数*/ count=1《<power: /*分配运算所需存储器* W=(COMPLEX *)malloc (sizeof (COMPLEX)*count/2) x1-(COMPLEX)mal(izof(COMPLEX)t) OMPLEX)malloc(sizeof(COMPLEX)*count) /计算加权系数* for(i=0:i<count/2:i++) W[i].im=sin(angle): /将时域点写入存储器*/ of(C)) for(k-0:k<power:k++) for(j=0:j<1<<k;j++) bfsize=1<<powe -k; for(i-0:i<bfsize/2:i++) D=i#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++)

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++) {

n=0 for(i-0:i<power:i++) if (j&(1<<i))p+=1<<power-i-1; FD[]=X1[p]: /释放存储器* free(W) free(X1): free(X2): /*快速付里哀反变换，利用快速付里哀变换 FD为频域值，TD为时域值，power为2的幂数*/ void IFFT(COMPLEX *FD,COMPLEX *TD,int power) int i,count; /计算付里哀反变换点数*/ count=1<<power: *分配云算所需存储婴* x=(COMPLEX)malloc(siz eof(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<cot:i+) TD[i].re/=count; TD[i],in-TD[i].im/count /*释放存储器* free(x): 附录二快速余弦变换与反变换程序实似

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++) { TD[i].re/=count; TD[i].im=-TD[i].im/count; } /*释放存储器*/ free(x); } 附录二、快速余弦变换与反变换程序实例

(利用2N点付里哀变换实现快速余弦变换) /快速离散余弦变换，利用快速付里哀变换 f为时域值，F为变换域值，power为2的幂数*/ void DCT(double *f,double #F,int power) int i count: COMPLEX *X: double s /计算离散余弦变换点数*/ count=1<<power: /分配运算所需存储器* X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2) /*延拓*/ memset (X,0,sizeof (COMPLEX)*count*2) /将时域点写入存储器* for(i=0:i<count:i++) x[i】.re=f[i] /调用快速付里哀变换* FFT(X,X,power+1): /*调整系数*/ s=1/sqrt count) F[O]=X[O]. re*s s*=sqrt(2) for(i=1:i<count:i++) F[i]=(仅[i].re*cos(i和i/(count*2)+X[i.i*sin(ipi/(counta*2)*s /体释放存储器*/ free): /快速离散余弦反变换，利用快速付里哀反变换 F为变换域值，f为时域值，poer为2的幂数* void IDCT(double *F,double *f,int power) int i count. COMPLEX *X: /计算离散余弦反变换点数*/ count=1<<power:

（利用 2N 点付里哀变换实现快速余弦变换） /*快速离散余弦变换,利用快速付里哀变换 f 为时域值,F 为变换域值,power 为 2 的幂数*/ void DCT(double *f, double *F, int power) { int i,count; COMPLEX *X; double s; /*计算离散余弦变换点数*/ count=1<<power; /*分配运算所需存储器*/ X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2); /*延拓*/ memset(X,0,sizeof(COMPLEX)*count*2); /*将时域点写入存储器*/ for(i=0;i<count;i++) { X[i].re=f[i]; } /*调用快速付里哀变换*/ FFT(X,X,power+1); /*调整系数*/ s=1/sqrt(count); F[0]=X[0].re*s; s*=sqrt(2); for(i=1;i<count;i++) { F[i]=(X[i].re*cos(i*pi/(count*2))+X[i].im*sin(i*pi/(count*2)))*s; } /*释放存储器*/ free(X); } /*快速离散余弦反变换,利用快速付里哀反变换 F 为变换域值,f 为时域值,power 为 2 的幂数*/ void IDCT(double *F, double *f, int power) { int i,count; COMPLEX *X; double s; /*计算离散余弦反变换点数*/ count=1<<power;

/*分配运算所需存储器*/ X=(COMPLEX*)malloc(sizeof(COMPLEX)*count*2) /*延拓* memset(化，0，sizeof(C0 MPLEX)*counta*2) /*调整频域点并写入存储器*/ for(i=0:i<count:i++) X[i].re=F[i]*cos(i*pi/(count*2)) X[i】.im=F[i]*sin(iw柳i/(count*2) /*调用快速付里哀反变换*/ IFFT(X,X.power+1): /调整系数* s=1/sqrt(count); for(i=0:i<count:i+】 f[i]=(1-s0rt(2))*s*F[0]+s0rt(2)*s*X[i].re*count*2: /释放存储器* free(X): 附录三.快速Walsh-dm变换与反变换程序实例 /*快速沃尔什-哈达玛变换 f为时域值，F为变换域值，power为2的幂数*/ void WALh(double*f,double*，int power) int count int i,j,k,bfsize,p: doub1e1,*2,*: /计算快速沃尔什变换点数* count=1<<power /*分配运算所需存储器* X1=(double *)malloc(sizeof(double)*count): X2=(double *)malloc(sizeof(double)*count) *将时域点写入存储器*/ for(k=0;k<power:k++) for(=0:i(1(k:i++)】

/*分配运算所需存储器*/ X=(COMPLEX *)malloc(sizeof(COMPLEX)*count*2); /*延拓*/ memset(X,0,sizeof(COMPLEX)*count*2); /*调整频域点并写入存储器*/ for(i=0;i<count;i++) { X[i].re=F[i]*cos(i*pi/(count*2)); X[i].im=F[i]*sin(i*pi/(count*2)); } /*调用快速付里哀反变换*/ IFFT(X,X,power+1); /*调整系数*/ s=1/sqrt(count); for(i=0;i<count;i++) { f[i]=(1-sqrt(2))*s*F[0]+sqrt(2)*s*X[i].re*count*2; } /*释放存储器*/ free(X); } 附录三、快速 Walsh-Hadamard 变换与反变换程序实例 /*快速沃尔什-哈达玛变换 f 为时域值,F 为变换域值,power 为 2 的幂数*/ void WALh(double *f, double *W, int power) { int count; int i,j,k,bfsize,p; double *X1,*X2,*X; /*计算快速沃尔什变换点数*/ count=1<<power; /*分配运算所需存储器*/ X1=(double *)malloc(sizeof(double)*count); X2=(double *)malloc(sizeof(double)*count); /*将时域点写入存储器*/ memcpy(X1,f,sizeof(double)*count); /*蝶形运算*/ for(k=0;k<power;k++) { for(j=0;j<1<<k;j++)