快速傅里叶变化:主要有几个难点
1.如果使用DIT首先要倒位序,网上的方法值得参考
思路是和二进制进位的思路相同只是是反着进位
2.要自己封装复数类型
3.难点是公式里的k不好求,因为k不是连续取得的,我用了一个数j来存储,这样可以把k计算出来
4.wn^r中r的取值,直接移位即可,不需要转化为二进制
5.输入不是2的倍数时自动补零
其余我都写在注释里了,感觉还是有很多地方要改进
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "复数类型.c"
//可以用M_E表示自然对数,M_PI表示圆周率
/*首先要导位序,也可以只要一个存储空间A[i],只需要改变A[i]中i
重点1.只改变当i小于j的时候的值,而i大于等于j时因为已经换过,不能再换回来了
重点2.思路是与二进制算法思路,小于N/2时即,最高位为0,此时的结果需要加上N/2
当大于N/2时,最高位为1,此时要向低位进位,意味着要先-N/2代表进位
再判断是否大于N/2,大于再向次低位进位,否则直接加N/4*/
void daoxu(double *x,int n)
{
int i=0,j=0,k=0;//i表示原始值,而j表示倒位序后的值
double temp;
for(;i < n;i++)
{
if(i<j)
{
temp = x[j];
x[j] = x[i];
x[i] = temp;
}
k = n/2;
while(j >= k && k != 0)
{
j = j - k;
k = k/2;
}
j = j + k;
}
/*for(i=0;i<n;i++)
{
printf("x[%d]=%f\n",i,x[i]);
//printf("result: new_k[%d]=%d\n",i,new_k[i]);
}*/
}
int process(int* ori_len)//对输入的数据进行处理,还要对数据进行补零,问题在于如何分配内存,使用链表?
{
int temp = 1,i;
for(i=0;;i++)
{
if(temp >= *ori_len)
break;
temp =temp * 2;
}
*ori_len = temp;
return i;//返回幂次方
}
complex w_n(int k,int m,int L,int length)
//如果这个值是复数意味着要包装复数类型,并且和这个数做四则运算的其他数也必须是复数
//需要新建复数这个数据结构,但是其实最后计算的是幅值
{
int r=k<<(L-m);//不应该取这个k,而应该取倒位序的k
r=r&(length-1);//去除高位
//printf("得到的r为=%d\n",r);
complex result;
result.real=cos(2*M_PI/length*r);//用欧拉公式展开
//printf("result.real=cos(2*pi/%d*%d)=%f\n",length,r,result.real);
result.imag=-sin(2*M_PI/length*r);
//printf("result.imag=-sin(2*pi/%d*%d)=%f\n",length,r,result.imag);
return result;
}
void butterfly_operation(double *x,int powe,int length)
{
int k_num = length/2,k_cal,k,l;
complex r;
complex *back_up=(complex*)calloc(length,sizeof(complex));//需要备份值来保留处理前的数据,定义为虚数类型
complex *y=(complex*)calloc(length,sizeof(complex));
for(int m=1;m<=powe;m++)//pow表示蝶形的次数
{
//printf("第%d次蝶形运算\n",m);
//printf("要做%d次蝶形运算\n",k_num);
k = 0;//运算之前将k清零,这里的k是实际是倒位序的
l = 0;
for(int i=0;i<length;i++)//运算完一次蝶形更新备份值
{
if(m==1)
back_up[i] = cpx_make(x[i],0);
else
back_up[i] = y[i];
//printf("back_up[i]=%f+%fi\n",back_up[i].real,back_up[i].imag);
}
k_cal=pow(2,m-1);//k_cal实际代表的是2^(m-1)的值
//printf("计算2^(m-1)的值=%d\n",k_cal);
for(int j=0;j<k_num;j++)//在做第m次蝶形的时候,要做k_num次蝶形运算
{
//printf("k=%d\n",k);
r = w_n(k,m,powe,length);
//printf("r=%f+%fi\n",r.real,r.imag);
y[k] =cpx_add(back_up[k],cpx_mul(back_up[k+k_cal],r));//x[k]取的是模值
//printf("y[%d]=%f+%fi\n",k,y[k].real,y[k].imag);
y[k+k_cal] = cpx_sub(back_up[k],cpx_mul(back_up[k+k_cal],r));
//printf("cpx_mul(back_up[k+k_cal],r)=%f+%fi\n",cpx_mul(back_up[k+k_cal],r).real,cpx_mul(back_up[k+k_cal],r).imag);
//printf("y[%d]=%f+%fi\n",k+k_cal,y[k+k_cal].real,y[k_cal+k].imag);
if(k_cal>=2&&k<l+k_cal-1)
{
k+=1;
}
else
{
k=k+k_cal+1;
l=k;
}
}
//printf("-------------------------------------------------------------------\n");
}
for(int i=0;i<length;i++)
{
x[i]=sqrt(y[i].real*y[i].real+y[i].imag*y[i].imag);
//printf("y[i].real=%f,y[i].imag=%f\n",y[i].real,y[i].imag);
}
free(back_up);
free(y);
}
int main()
{
int length,power,old_length;//以后要更正
double t;
printf("请输入采样的长度length=");
scanf("%d",&length);//输入采样的长度*/
double *x,*y;
y=(double*)calloc(length,sizeof(double));//可以动态分配内存
for(int i=0;i<length;i++)
{
//y[i] = i+1;//在输入之前就会分配好内存空间因此不能临时输入
t=i*0.001;
y[i] = 1.2+2.7*cos(2*M_PI*33*t)+5*cos(2*M_PI*200*t+M_PI/2);
//y[i]=0.6*sin(2*M_PI*500*i)+0.6*sin(2*M_PI*50*i);
//y[i] = sin(M_PI*i/8);
//printf("%f\n",y[i]);
}
old_length = length;
power = process(&length);//改变了length
x=(double*)calloc(length,sizeof(double));
for(int i=0;i<old_length;i++)
{
x[i] = y[i];
printf("x[i]=%f\n",x[i]);
}
for(int i=old_length;i<length;i++)
{
x[i]=0;
printf("x[i]=%f\n",x[i]);
}
free(y);
daoxu(x,length);
butterfly_operation(x,power,length);
for(int i=0;i<length;i++)
printf("%d\t%f\n",i,x[i]);//在输入之前就会分配好内存空间因此不能临时输入
free(x);
return 0;
}
我自己封装的复数头文件
#include <stdio.h>
#define Real(c) (c).real
#define Imag(c) (c).imag
typedef struct
{
double real;
double imag;
}complex;
complex cpx_make(double real, double imag)
{
complex ret;
ret.real = real;
ret.imag = imag;
return ret;
}
complex cpx_add(complex a,complex b)
{
return cpx_make(Real(a) + Real(b), Imag(a) + Imag(b));
}
complex cpx_sub(complex a,complex b)
{
return cpx_make(Real(a) - Real(b), Imag(a) - Imag(b));
}
complex cpx_mul(complex a,complex b)
{
return cpx_make(Real(a) * Real(b)-Imag(a) * Imag(b),Real(a) * Imag(b) + Imag(a) * Real(b));
}
complex cpx_div(complex a,complex b)
{
Real(b) = 1/Real(b);
Imag(b) = 1/Imag(b);
return cpx_mul(a,b);
}