信号处理 FFT与IFFT的实际应用程序设计与调试

实验4 FFT与IFFT的实际应用程序设计与调试
一、实验目的
掌握二维FFT算法原理，掌握利用二维FFT与IFFT实现去除燥声的方法，理解FFT变换在信号处理中的作用。

二、实验原理
1．二维FFT（IFFT）算法原理
在图像处理的广泛领域中，傅立叶变换起着非常重要的作用，包括图像的增强、图像分析、图像复原和图像压缩等。

在图像数据的数字处理中常用的是二维离散傅立叶变换，它能把空间域的图像转变到频域上进行研究，从而很容易地了解到图像各空间频域成分，进行相应处理。

二维FFT（IFFT）实际上就是对二维数组的行和列分别进行一维FFT（IFFT）。

其流程如下：
①先逐行进行一维FFT（IFFT），计算结果存入矩阵中，
②再逐列进行一维FFT（IFFT），即先取得第j列的数据col[i] = x[i][j]；对
第j列数据进行一维FFT（IFFT）计算结果在数组col中，将结果放回矩
阵第j列中。

2．文本数据转化成数值
用．和O组成的“国”字图形存储在文本文件中，需要转换成数值才能进行FFT变换。

其处理方法如下：
首先设置点(.)对应的数值POINT_V AL=0，圈(o)对应的数值
CIRCLE_V AL=1.0；
然后，初始化Data[][]用于存储图像数据的复数矩阵；
然后打开文件，读取文件中的文本数据，每次读取一个字符，如果为'o'，则Data[i][j].real=CIRCLE_V AL；如果如果为'．'，则Data[i][j].real=POINT_V AL；
最后，关闭文件。

3．数值转化成文本数据
经过处理后的Data[][]中的数值数据，需要转换成用．和O组成的“国”字图形，才知其是否已经去除噪声。

其处理方法如下：
首先设置数值转化为字符的控制阈值TH=0.5；
然后，打开文件，依次读取Data[][]中数值，如果大于TH，向文件输出'o'，如果小于TH，向文件输出'．'；
最后，关闭文件
4．去除燥声
用"．"和"o"组成的“国”字图形，个别地方有燥声，调用二维FFT与IFFT 实现对它进行二维FFT变换与反变换还得到原图形，适当调整参数可在反变换后去除燥声。

处理流程如下：
①打开数据文件（Data.txt）读取数据，并将文本数据转化成数值，存入矩
阵中
②进行二维FFT变换
③对频谱进行处理
④进行二维FFT的逆变换
⑤将变换结果转换成文本，存人结果文件（RES.txt）中
三、实验内容
编程实现去除图形中的燥声。

文件Data.txt中是用“．”和“o”组成的“国”字图形，个别地方有燥声，编成实现去除燥声。

用文件res.txt存放去除燥声后的图形。

四、实验代码
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>
#define PI (4.0*atan(1.0))
#define MAX_N 64
/* 复数的数据类型*/
typedef struct {
float real, imag;
} COMPLEX;
/*倒序的实现*/
void ReverseOrder(COMPLEX A[],int N)
{
int NV2=N/2;
int NM1=N-1;
int I,J,K=0;
COMPLEX T;
I=J=1;
while(I<=NM1)
{
if(I< J)
{
/*借助于中间变量T，将A[J－1]的内容和A[I－1]的内容互换*/
T=A[J-1];
A[J-1]=A[I-1];
A[I-1]=T;
}
/*求倒码*/
K=NV2;
while(K< J)
{
J-=K;
K/=2;
}
J+=K;
I++;
}
}
/*傅立叶正变换*/
void fft(COMPLEX A[],int M)
{
COMPLEX U,W,T;
int LE,LE1,I,J,IP;
int N=(int)pow(2,M);
int L;
float temp;
/*由于采用时间抽选奇偶分解方式，所以在参加运算前首先要对时间序列进行倒序*/
ReverseOrder(A,N);
L=1;
while(L<=M)
{
LE=(int)pow(2,L);
LE1=LE/2;
U.real=1.0f;
U.imag=0.0f;
/*计算W算子的值*/
W.real=(float)cos(PI/(1.0*LE1));
W.imag=(float)-1.0*sin(PI/(1.0*LE1));
/*if(abs(W.real)<1.0e-12)
W.real=0.0f;
if(abs(W.imag)< 1.0e-12)
W.imag=0.0f;*/
J=1;
while(J<=LE1)
{
I=J;
while(I<=N)
{
IP=I+LE1;
/*复数运算A×U*/
T.real=(float)A[IP-1].real*U.real-A[IP-1].imag*U.imag;
T.imag=(float)A[IP-1].real*U.imag+A[IP-1].imag*U.real;
/*复数运算A－T*/
A[IP-1].real=(float)A[I-1].real-T.real;
A[IP-1].imag=(float)A[I-1].imag-T.imag;
/*复数运算A＋T*/
A[I-1].real+=T.real;
A[I-1].imag+=T.imag;
I+=LE;
}
temp=U.real;
/*复数运算U×W*/
U.real=(float)U.real*W.real-U.imag*W.imag;
U.imag=(float)temp*W.imag+U.imag*W.real;
J++;
}
L++;
}
}
/*傅立叶反变换*/
void ifft(COMPLEX A[],int M)
{
COMPLEX U,W,T;
int LE,LE1,I,J,IP;
int N=(int)pow(2,M);
int L;
float temp,scale;
/*由于采用时间抽选奇偶分解方式，所以在参加运算前首先要对时间序列进行倒序*/
ReverseOrder(A,N);
L=1;
while(L<=M)
{
LE=(int)pow(2,L);
LE1=LE/2;
U.real=1.0f;
U.imag=0.0f;
/*计算W算子的值*/
W.real=(float)cos(PI/(1.0*LE1));
W.imag=(float)1.0*sin(PI/(1.0*LE1));
/*if(abs(W.real)<1.0e-12)
W.real=0.0f;
if(abs(W.imag)< 1.0e-12)
W.imag=0.0f;*/
J=1;
while(J<=LE1)
{
I=J;
while(I<=N)
{
IP=I+LE1;
/*复数运算A×U*/
T.real=(float)A[IP-1].real*U.real-A[IP-1].imag*U.imag;
T.imag=(float)A[IP-1].real*U.imag+A[IP-1].imag*U.real;
/*复数运算A－T*/
A[IP-1].real=(float)A[I-1].real-T.real;
A[IP-1].imag=(float)A[I-1].imag-T.imag;
/*复数运算A＋T*/
A[I-1].real+=T.real;
A[I-1].imag+=T.imag;
I+=LE;
}
temp=U.real;
/*复数运算U×W*/
U.real=(float)U.real*W.real-U.imag*W.imag;
U.imag=(float)temp*W.imag+U.imag*W.real;
J++;
}
L++;
}
/* 反变换乘1/N */
scale = (float)(1.0/N);
for(I = 0 ; I < N ; I++)
{
A[I].real = scale*A[I].real;
A[I].imag = scale*A[I].imag;
}
}
/*二维傅立叶正变换*/
void fft2(COMPLEX A[MAX_N][MAX_N],int M)
{
int N=(int)pow(2,M);
COMPLEX col[MAX_N];
int i,j;
/*先逐行进行一维-FFT*/
for (i=0; i<N; i++)
fft(A[i], M); /* 计算结果再存入矩阵A中*/
/*再逐列进行一维-FFT*/
for (j=0; j<N; j++)
{
/*取得第j列的数据*/
for ( i=0; i<N; i++)
col[i] = A[i][j];
/*对第j列数据进行一维-FFT*/
fft(col, M); /* 计算结果在数组col中*/
/* 将结果放回矩阵第j列中*/
for (i=0; i<N; i++)
A[i][j] = col[i];
}
}
/*二维傅立叶反变换*/
void ifft2(COMPLEX A[MAX_N][MAX_N],int M)
{
int N=(int)pow(2,M);
COMPLEX col[MAX_N];
int i,j;
/*先逐行进行一维-FFT*/
for (i=0; i<N; i++)
ifft(A[i], M); /* 计算结果再存入矩阵A中*/ /*再逐列进行一维-FFT*/
for (j=0; j<N; j++)
{
/*取得第j列的数据*/
for ( i=0; i<N; i++)
col[i] = A[i][j];
/*对第j列数据进行一维-FFT*/
ifft(col, M); /* 计算结果在数组col中*/
/* 将结果放回矩阵第j列中*/
for (i=0; i<N; i++)
A[i][j] = col[i];
}
}
void main()
{
/*设定"文本--> 数值" 的参数*/
float POINT_V AL=0.0; /*点(.)对应的数值*/
float CIRCLE_V AL=1.0; /*圈(o)对应的数值*/
float TH=0.5; /*数值转化为字符的控制阈值*/
double tempValue;
COMPLEX Data[MAX_N][MAX_N];/* 存储图像数据的复数矩阵*/ int i,j,N,M;
char inputchar;
FILE *fp;
if((fp=fopen("e:\\FFT\\Data.c","r"))==NULL)
{
printf("Can not Open File");
exit(0);
}
/*将文本数据转化成数值，存入矩阵*/
N=64;
for (i=0; i<N; i++)
{
for (j=0; j<N; j++)
{
inputchar=fgetc(fp);
//printf("%c",inputchar);
if (inputchar=='\n')
inputchar=fgetc(fp);
if (inputchar=='o')
{
Data[i][j].real=CIRCLE_V AL;
Data[i][j].imag=0.0;
}
else if (inputchar=='.')
{
Data[i][j].real=POINT_V AL;
Data[i][j].imag=0.0;
}
}
}
fclose(fp);
/*完成对数据的2D-FFT的正变换、频谱处理和逆变换*/
M=log(N)/log(2)+1;
fft2(Data, M); //正变换
printf("M=%d\n",M);
for (i=0; i<N; i++)
{
for ( j=0; j<N; j++)
{
tempValue=sqrt(Data[i][j].real*Data[i][j].real+Data[i][j].imag*Data[i][j].imag);
if(tempValue<10.0)
{
Data[i][j].real=0;
Data[i][j].imag=0;
}
}
printf("\n");
}
ifft2(Data, M); //逆变换
if((fp=fopen("e:\\FFT\\res.txt","w"))==NULL)
{
printf("Can not Open File");
exit(0);
}
for (i=0; i<N; i++)
{
for ( j=0; j<N; j++)
{
if ((Data[i][j].real)>TH)
fputc('o',fp);
else{
fputc('.',fp);
fputc(' ',fp);
} //机器编码问题，"o"占两个字符，"."占一个字符，需要补个空格
}
fputc('\n',fp);
}
fclose(fp);
}
五、实验截图。