📄 lab2_fft.c
字号:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define PI 3.1415926535898
#define ROW 128
#define COL 128
#define true 1
#define false 0
#define LEVEL 127
#define coe 0x1000 //coefficient for integer calculating
int n_fft_int(unsigned char *inimage, int *prdata, int *pidata, int *ap,
int row, int col);
int fft_inline_asm(int n, int *ar, int *ai);
int fft_asm(int n, int*ar, int*ai);
int sqr_int(int ff);
int FileRead(char *filename, unsigned char *data, int length);
int FileWrite(char *filename, unsigned char *data, int length);
int n_fft(unsigned char *inimage, double *prdata,double *pidata,double *ap,
int row, int col);
int fft(int n, double *ar, double *ai);
void swap(double *f1, double *f2);
void swap_int(int *f1, int *f2);
double sqr(double ff);
int SIN[128]; //sin(x*PI/128) << 12 (* 4096)
int COS[128]; //cos(x*PI/128) << 12 (* 4096)
int main()
{
int i = 0, j = 0;
unsigned char image[ROW*COL];
int rdata[ROW*COL];
int idata[ROW*COL];
int amplitude[ROW*COL];
/* double rdata_double[ROW*COL];
double idata_double[ROW*COL];
double amplitude_double[ROW*COL];
*/
unsigned char output[ROW*COL];
for(i = 0;i < 128;i++)
{
SIN[i] = (int)(sin(i*PI/128)*coe);
COS[i] = (int)(cos(i*PI/128)*coe);
}
if(false == FileRead("d:/fft_data/lab2_d.raw", image, ROW*COL) )
{
printf("open file failed!\n");
}
else
{
printf("open file successful!\n");
}
#if 0
n_fft(image, rdata_double, idata_double, amplitude_double, ROW, COL);
for(i = 0;i < ROW;i++)
{
for(j = 0;j < COL;j++)
{
output[i*COL + j] = (unsigned char)amplitude_double[i*COL + j];
}
}
if(false == FileWrite("d:/fft_data/lab2_d_C_FFT.raw", output, ROW*COL) )
{
printf("Write file failed!\n");
}
else
{
printf("Write file successful!\n");
}
#endif
n_fft_int(image, rdata, idata, amplitude, ROW, COL);
for(i = 0;i < ROW;i++)
{
for(j = 0;j < COL;j++)
{
output[i*COL + j] = (unsigned char)amplitude[i*COL + j];
}
}
if(false == FileWrite("d:/fft_data/lab2_d_ASM_FFT.raw", output, ROW*COL) )
{
printf("Write file failed!\n");
}
else
{
printf("Write file successful!\n");
}
return false;
}
//whole fft
#if 0
#define INTEGER_FFT
#else
#define ASM_FFT
#endif
int n_fft_int(unsigned char *inimage, int *prdata, int *pidata, int *ap,
int row, int col)
{
int i,j;
int gmax;
int *ar;
int *ai;
ar=malloc(2*col*sizeof(int));
ai=malloc(2*col*sizeof(int));
#if 1
for(i=0;i<row;i++)
{
for(j=0;j<col;j++)
{
*(prdata+i*col+j)=(int)*(inimage+i*col+j);
*(pidata+i*col+j)=0;
}
}
#else
for(i=0;i<row;i++)
{
for(j=0;j<col;j++)
{
(prdata[i*col+j])=(int)(inimage[i*col+j]);
(pidata[i*col+j])=0;
}
}
#endif
for(i=0;i<row;i++)
{
for(j=0;j<col;j++)
{
*(ar+j)=*(prdata+i*col+j);
*(ai+j)=*(pidata+i*col+j);
}
#ifdef INTEGER_FFT
fft_int(col,ar,ai);
#endif
#ifdef ASM_FFT
fft_asm(col,ar,ai);
#endif
for(j=0;j<col;j++)
{
*(prdata+i*col+j)=*(ar+j);
*(pidata+i*col+j)=*(ai+j);
}
}
for(i=0;i<col;i++)
{
for(j=0;j<row;j++)
{
*(ar+j)=*(prdata+j*col+i);
*(ai+j)=*(pidata+j*col+i);
}
#ifdef INTEGER_FFT
fft_int(row,ar,ai);
#endif
#ifdef ASM_FFT
fft_asm(row,ar,ai);
#endif
for(j=0;j<row;j++)
{
*(prdata+j*col+i)=*(ar+j);
*(pidata+j*col+i)=*(ai+j);
}
}
//printf("the end\n");
gmax=0;
for(i=0;i<row;i++)
{
for(j=1;j<col;j++)
{
*(ap+col*i+j)=sqrt(sqr_int(*(prdata+i*col+j))+sqr_int(*(pidata+i*col+j)));
if(gmax<*(ap+col*i+j))
gmax=*(ap+col*i+j);
}
}
for(i=0;i<row;i++)
{
for(j=1;j<col;j++)
{
*(ap+i*col+j)=(unsigned char)(*(ap+col*i+j)*LEVEL/gmax);
}
}
free(ar);
free(ai);
return 1;
}
void swap_int(int *f1, int *f2)
{
int ttf;
ttf= *f1;
*f1= *f2;
*f2=ttf;
return;
}
//Interger fast fourier convert
int fft_int(int n, int *ar, int *ai)
{
int n2, a, c, d, f, g, h, j;
int b, e, k, i, n1, co, si;
int m, p, q, r, s;
n1=log10(n)/log10(2); // n1 = 7; n = 128
n2=-1; // n2 = -1;
a=n;
//b=2*PI*coe/n;
b=1;
for(c=1;c<=n1;c++)
{
d=a;
a=a/2;
b=b*2;
e=0;
for(f=0;f<a;f++)
{
co=COS[e]; //cos(e) * coe
si=SIN[e]*n2; // sin(e) * coe
e=e+b; // step 2PI/128
for(g=d;g<=n;g+=d)
{
h=g-d+f;
j=h+a;
k=*(ar+h)-*(ar+j);
i=*(ai+h)-*(ai+j);
*(ar+h)=*(ar+h)+*(ar+j);
*(ai+h)=*(ai+h)+*(ai+j);
*(ar+j)=(co*k+si*i)/coe;
*(ai+j)=(co*i-si*k)/coe;
}
}
}
m=0;
p=n/2;
q=n-1;
for(r=0;r<q;r++)
{
if((r-m)<0)
{
swap_int((int*)ar+r, (int*)ar+m);
swap_int((int*)ai+r, (int*)ai+m);
}
s=p;
while((s-m)<1)
{
m=m-s;
s=s/2;
}
m=m+s;
}
#if 0
for(s=0;s<n;s++)
{
*(ar+s)=*(ar+s)/n;
*(ai+s)=*(ai+s)/n;
}
#else
for(s=0;s<n;s++)
{
*(ar+s)=*(ar+s)/1;
*(ai+s)=*(ai+s)/1;
}
#endif
return 1;
}
//ASM fast fourier convert
int fft_asm(int n, int *ar, int *ai)
{
int n2, a, c, d, f, g, h, j;
int b, e, k, i, n1, co, si;
int m, p, q, r, s, t;
n1=log10(n)/log10(2); // n1 = 7; n = 128
n2=-1; // n2 = -1;
__asm
{
MOV a, n //a=n;
MOV b,1 //b=1;
}
for(c=1;c<=n1;c++)
{
__asm
{
MOV d,a //d=a;
MOV a,a,LSR 1 //a=a/2;
MOV b,b,LSL 1 //b=b*2;
MOV e,0 //e=0;
}
for(f=0;f<a;f++)
{
__asm
{
LDR co,[COS,e,LSL 2] //co=COS[e]; //cos(e) * coe
LDR si,[SIN,e,LSL 2] //si=SIN[e]*n2; // sin(e) * coe
ADD e,e,b //e=e+b; // step 2PI/128
}
for(g=d;g<=n;g+=d)
{
__asm
{
ADD h,g,f //h=g-d+f;
SUB h,h,d
ADD j,h,a //j=h+a;
LDR k,[ar,h,LSL 2] //k=*(ar+h)-*(ar+j);
LDR r,[ar,j,LSL 2]
SUB k,k,r
LDR i,[ai,h,LSL 2] //i=*(ai+h)-*(ai+j);
LDR r,[ai,j,LSL 2]
SUB i,i,r
LDR r,[ar,h,LSL 2] //*(ar+h)=*(ar+h)+*(ar+j);
LDR s,[ar,j,LSL 2]
ADD r,r,s
STR r,[ar,h,LSL 2]
LDR r,[ai,h,LSL 2] //*(ai+h)=*(ai+h)+*(ai+j);
LDR s,[ai,j,LSL 2]
ADD r,r,s
STR r,[ai,h,LSL 2]
MUL r,co,k //*(ar+j)=(co*k+si*i)/coe;
MUL s,si,i
ADD t,r,s
MOV t,t,ASR 12
STR t,[ar,j,LSL 2]
MUL r,co,i //*(ai+j)=(co*i-si*k)/coe;
MUL s,si,k
SUB t,r,s
MOV t,t,ASR 12
STR t,[ai,j,LSL 2]
}
}
}
}
m=0;
p=n/2;
q=n-1;
for(r=0;r<q;r++)
{
if((r-m)<0)
{
swap_int((int*)ar+r, (int*)ar+m);
swap_int((int*)ai+r, (int*)ai+m);
}
s=p;
while((s-m)<1)
{
m=m-s;
s=s/2;
}
m=m+s;
}
#if 0
for(s=0;s<n;s++)
{
*(ar+s)=*(ar+s)/n;
*(ai+s)=*(ai+s)/n;
}
#else
for(s=0;s<n;s++)
{
*(ar+s)=*(ar+s)/1;
*(ai+s)=*(ai+s)/1;
}
#endif
return 1;
}
//asm fast fourier convert
int fft_inline_asm(int n, int *ar, int *ai)
{
int a, b, d, e;
int c, f, g;
int h, j;
int k, i;
int co, si;
int m, r, s, t;
__asm
{
MOV a,128 //a=128;
MOV b,1 //b=1;
MOV c,1 //for(c=1;c<=7;c++)
LOOP_C: //{
MOV d,a // d=a;
MOV a,a,ASR 1 // a=a/2; //?LSR a,a,1
MOV b,b,LSL 1 // b=2*b; //?LSL b,b,1
MOV e,0 // e=0;
MOV f,0 // for(f=0;f<a;f++)
LOOP_F: // {
LDR co,[COS,e,LSL 2]// co=COS[e];
LDR si,[SIN,e,LSL 2]// si=SIN[e];//n2;
ADD e,e,b // e=e+b;
MOV g,d // for(g=d;g<=128;g+=d)
LOOP_G: // {
ADD h,g,f // h=g-d+f;
SUB h,h,d
ADD j,h,a // j=h+a;
LDR k,[ar,h,LSL 2] // k=*(ar+h)-*(ar+j);
LDR r,[ar,j,LSL 2]
SUB k,k,r
LDR i,[ai,h,LSL 2] // i=*(ai+h)-*(ai+j);
LDR r,[ai,j,LSL 2]
SUB i,i,r
LDR r,[ar,h,LSL 2] // *(ar+h)=*(ar+h)+*(ar+j);
LDR s,[ar,j,LSL 2]
ADD r,r,s
STR r,[ar,h,LSL 2]
LDR r,[ai,h,LSL 2] // *(ai+h)=*(ai+h)+*(ai+j);
LDR s,[ai,j,LSL 2]
ADD r,r,s
STR r,[ai,h,LSL 2]
MUL r,co,k // *(ar+j)=(co*k-si*i)/1024/4;//256/128;
MUL s,si,i
SUB t,r,s
MOV t,t,ASR 12
STR t,[ar,j,LSL 2]
MUL r,co,i // *(ai+j)=(co*i+si*k)/1024/4;//256/128;
MUL s,si,k
ADD t,r,s
MOV t,t,ASR 12
STR t,[ai,j,LSL 2]
ADD g,g,d // }
CMP g,128
BLE LOOP_G
ADD f,f,1 // }
CMP f,a
BLT LOOP_F
ADD c,c,1 //}
CMP c,7
BLE LOOP_C
MOV m,0 //m=0;
//MOV p,64 //p=128/2;
//MOV q,127 //q=128-1;
MOV r,0 //for(r=0;r<q;r++)
LOOP_R: //{
CMP r,m // if((r-m)<0)
// {
LDRLT s,[ar,r,LSL 2] // swap((int *)ar+r, (int *)ar+m);
ADDLT t,ar,m,LSL 2
SWPLT s,s,[t]
STRLT s,[ar,r,LSL 2]
LDRLT s,[ai,r,LSL 2] // swap((int *)ai+r, (int *)ai+m);
ADDLT t,ai,m,LSL 2
SWPLT s,s,[t]
STRLT s,[ai,r,LSL 2]
// }
MOV s,64 // s=p;
BEGIN_WHILE:
CMP s,m // while((s-m)<1) s<m+1 s<=m //?WHILE s < m ... WEND
BGT END_WHILE
// {
SUB m,m,s // m=m-s;
MOV s,s,ASR 1 // s=s/2;
B BEGIN_WHILE
END_WHILE: // }
ADD m,m,s // m=m+s;
ADD r,r,1 //}
CMP r,127
BLT LOOP_R
MOV s,0 //for(s=0;s<128;s++)
LOOP_S: //{
LDR t,[ar,s,LSL 2] // *(ar+s)=*(ar+s)/8;
MOV t,t,ASR 3
STR t,[ar,s,LSL 2]
LDR t,[ai,s,LSL 2] // *(ai+s)=*(ai+s)/8;
MOV t,t,ASR 3
STR t,[ai,s,LSL 2]
ADD s,s,1 //}
CMP s,128
BLT LOOP_S
}
return 1;
}
//square
int sqr_int(int ff)
{
return (ff*ff);
}
int FileRead(char *filename, unsigned char *data, int length)
{
FILE *fd;
if ((fd = fopen(filename, "rb")) == NULL)
{
return false;
}
fread(data, sizeof(unsigned char), length, fd);
fclose(fd);
return true;
}
int FileWrite(char *filename, unsigned char *data, int length)
{
FILE *fd;
if ((fd = fopen(filename, "wb")) == NULL)
{
return false;
}
fwrite(data, sizeof(unsigned char), length, fd);
fclose(fd);
return true;
}
//double fast fourier convert
int fft(int n, double *ar, double *ai)
{
int n2, a, c, d, f, g, h, j;
double b, e, k, i, n1, co, si;
int m, p, q, r, s;
n1=log10(n)/log10(2);
n2=-1;
a=n;
b=2*PI/n;
for(c=1;c<=n1;c++)
{
d=a;
a=a/2;
e=0;
for(f=0;f<a;f++)
{
co=cos(e);
si=sin(e)*n2;
e=e+b;
for(g=d;g<=n;g+=d)
{
h=g-d+f;
j=h+a;
k=*(ar+h)-*(ar+j);
i=*(ai+h)-*(ai+j);
*(ar+h)=*(ar+h)+*(ar+j);
*(ai+h)=*(ai+h)+*(ai+j);
*(ar+j)=co*k+si*i;
*(ai+j)=co*i-si*k;
}
}
b=2*b;
}
m=0;
p=n/2;
q=n-1;
for(r=0;r<q;r++)
{
if((r-m)<0)
{
swap((double *)ar+r, (double *)ar+m);
swap((double *)ai+r, (double *)ai+m);
}
s=p;
while((s-m)<1)
{
m=m-s;
s=s/2;
}
m=m+s;
}
for(s=0;s<n;s++)
{
*(ar+s)=*(ar+s)/n;
*(ai+s)=*(ai+s)/n;
}
return 1;
}
/* f1, f2 swap */
void swap(double *f1, double *f2)
{
double ttf;
ttf= *f1;
*f1= *f2;
*f2=ttf;
return;
}
//
double sqr(double ff)
{
return(ff*ff);
}
//whole fft
int n_fft(unsigned char *inimage, double *prdata,double *pidata,double *ap,
int row, int col)
{
int i,j;
double *ar;
double *ai;
double gmax;
ar=malloc(2*col*sizeof(double));
ai=malloc(2*col*sizeof(double));
for(i=0;i<row;i++)
for(j=0;j<col;j++)
{
*(prdata+i*col+j)=(double)*(inimage+i*col+j);
*(pidata+i*col+j)=0;
}
for(i=0;i<row;i++)
{
for(j=0;j<col;j++)
{
*(ar+j)=*(prdata+i*col+j);
*(ai+j)=*(pidata+i*col+j);
}
fft(col,ar,ai);
for(j=0;j<col;j++)
{
*(prdata+i*col+j)=*(ar+j);
*(pidata+i*col+j)=*(ai+j);
}
}
for(i=0;i<col;i++)
{
for(j=0;j<row;j++)
{
*(ar+j)=*(prdata+j*col+i);
*(ai+j)=*(pidata+j*col+i);
}
fft(row,ar,ai);
for(j=0;j<row;j++)
{
*(prdata+j*col+i)=*(ar+j);
*(pidata+j*col+i)=*(ai+j);
}
}
//printf("the end\n");
gmax=0;
for(i=0;i<row;i++)
for(j=1;j<col;j++)
{
*(ap+col*i+j)=sqrt(sqr(*(prdata+i*col+j))+sqr(*(pidata+i*col+j)));
if(gmax<*(ap+col*i+j))
gmax=*(ap+col*i+j);
}
for(i=0;i<row;i++)
for(j=1;j<col;j++)
{
*(ap+i*col+j)=(unsigned char)(*(ap+col*i+j)*LEVEL/gmax);
}
free(ar);
free(ai);
return 1;
}
// 函数名: 快速傅立叶变换（来源《C常用算法集》）
// 本函数测试OK,可以在TC2.0,VC++6.0,Keil C51测试通过。
// 如果你的MCS51系统有足够的RAM时,可以验证一下用单片机处理FFT有多么的慢。
//
// 入口参数：
// l: l = 0, 傅立叶变换; l = 1, 逆傅立叶变换
// il: il = 0,不计算傅立叶变换或逆变换模和幅角；il = 1,计算模和幅角
// n: 输入的点数，为偶数，一般为32，64，128，...,1024等
// k: 满足n=2^k(k>0),实质上k是n个采样数据可以分解为偶次幂和奇次幂的次数
// pr[]: l=0时，存放N点采样数据的实部
// l=1时, 存放傅立叶变换的N个实部
// pi[]: l=0时，存放N点采样数据的虚部
// l=1时, 存放傅立叶变换的N个虚部
//
// 出口参数：
// fr[]: l=0, 返回傅立叶变换的实部
// l=1, 返回逆傅立叶变换的实部
// fi[]: l=0, 返回傅立叶变换的虚部
// l=1, 返回逆傅立叶变换的虚部
// pr[]: il = 1,i = 0 时，返回傅立叶变换的模
// il = 1,i = 1 时，返回逆傅立叶变换的模
// pi[]: il = 1,i = 0 时，返回傅立叶变换的辐角
// il = 1,i = 1 时，返回逆傅立叶变换的辐角
// data: 2005.8.15,Mend Xin Dong
void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
for (it=0; it<=n-1; it++)
{
m = it;
is = 0;
for(i=0; i<=k-1; i++)
{
j = m/2;
is = 2*is+(m-2*j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
//----------------------------
pr[0] = 1.0;
pi[0] = 0.0;
p = 6.283185306/(1.0*n);
pr[1] = cos(p);
pi[1] = -sin(p);
if (l!=0)
{
pi[1]=-pi[1];
}
for (i=2; i<=n-1; i++)
{
p = pr[i-1]*pr[1];
q = pi[i-1]*pi[1];
s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
pr[i] = p-q;
pi[i] = s-p-q;
}
for (it=0; it<=n-2; it=it+2)
{
vr = fr[it];
vi = fi[it];
fr[it] = vr+fr[it+1];
fi[it] = vi+fi[it+1];
fr[it+1] = vr-fr[it+1];
fi[it+1] = vi-fi[it+1];
}
m = n/2;
nv = 2;
for (l0=k-2; l0>=0; l0--)
{
m = m/2;
nv = 2*nv;
for(it=0; it<=(m-1)*nv; it=it+nv)
for (j=0; j<=(nv/2)-1; j++)
{
p = pr[m*j]*fr[it+j+nv/2];
q = pi[m*j]*fi[it+j+nv/2];
s = pr[m*j]+pi[m*j];
s = s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr = p-q;
poddi = s-p-q;
fr[it+j+nv/2] = fr[it+j]-poddr;
fi[it+j+nv/2] = fi[it+j]-poddi;
fr[it+j] = fr[it+j]+poddr;
fi[it+j] = fi[it+j]+poddi;
}
}
if(l!=0)
for(i=0; i<=n-1; i++)
{
fr[i] = fr[i]/(1.0*n);
fi[i] = fi[i]/(1.0*n);
}
if(il!=0)
{
for(i=0; i<=n-1; i++)
{
pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
if(fabs(fr[i])<0.000001*fabs(fi[i]))
{
if ((fi[i]*fr[i])>0)
pi[i] = 90.0;
else
pi[i] = -90.0;
}
else
pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;
}
}
return;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -