代码转自:https://my.oschina.net/VenusV/blog/1797459
#include <time.h>
#define WIDTH 5
#define HEIGHT 5
struct Complex_{
double real;//实部
double imagin;//虚部
};
typedef struct Complex_ DFT_Complex;
#define VALUE_MAX 255
#define M_PI 3.1415926535897932
int Initdata(double (*src)[WIDTH],int size_w,int size_h){
srand((int)time(0));
for(int i=0;i<size_w;i++){
for(int j=0;j<size_h;j++){
src[i][j]=rand()%VALUE_MAX;
printf("%lf ",src[i][j]);
}
printf(";\n");
}
return 0;
}
//2维傅里叶变换函数
int DFT2D(double (*src)[WIDTH],DFT_Complex (*dst)[WIDTH],int size_w,int size_h){
for(int u=0;u<size_w;u++){
for(int v=0;v<size_h;v++){
double real=0.0;
double imagin=0.0;
for(int i=0;i<size_w;i++){
for(int j=0;j<size_h;j++){
double I=src[i][j];
double x=M_PI*2*((double)i*u/(double)size_w+(double)j*v/(double)size_h);
real+=cos(x)*I;
imagin+=-sin(x)*I;
}
}
dst[u][v].real=real;
dst[u][v].imagin=imagin;
if(imagin>=0)
printf("%lf+%lfj ",real,imagin);
else
printf("%lf%lfj ",real,imagin);
}
printf(";\n");
}
return 0;
}
//2维逆傅里叶变换函数
int IDFT2D(DFT_Complex (*src)[WIDTH],DFT_Complex (*dst)[WIDTH],int size_w,int size_h){
for(int i=0;i<size_w;i++){
for(int j=0;j<size_h;j++){
double real=0.0;
double imagin=0.0;
for(int u=0;u<size_w;u++){
for(int v=0;v<size_h;v++){
double R=src[u][v].real;
double I=src[u][v].imagin;
double x=M_PI*2*((double)i*u/(double)size_w+(double)j*v/(double)size_h);
real+=R*cos(x)-I*sin(x);
imagin+=I*cos(x)+R*sin(x);
}
}
dst[i][j].real=(1./(size_w*size_h))*real;
dst[i][j].imagin=(1./(size_w*size_h))*imagin;
if(imagin>=0)
printf("%lf+%lfj ",dst[i][j].real,dst[i][j].imagin);
else
printf("%lf%lfj ",dst[i][j].real,dst[i][j].imagin);
}
printf(";\n");
}
return 0;
}
int _tmain(int argc, _TCHAR* argv[])
{
double src[WIDTH][HEIGHT];
DFT_Complex dst[WIDTH][HEIGHT];
DFT_Complex dst_[WIDTH][HEIGHT];
Initdata(src, WIDTH, HEIGHT);
printf("\n\n");
DFT2D(src,dst,WIDTH,HEIGHT);
printf("\n\n");
IDFT2D(dst,dst_,WIDTH,HEIGHT);
system("pause");
return 0;
}
对应DCT和IDCT公式可和代码一一对应。
也就是原始图像仅仅是实数,DCT后变为同样大小的复数。
IDCT输入复数,输出也是复数,但是实际虚部都为0
如果要滤波的话就是在DCT后的复数直接和滤波器相乘。