徐士良老师编写的c语言算法程序下载链接：https://pan.baidu.com/s/1zDV6iLeYeXmZaoZlP4yRAA
提取码：8opo
一：徐士良老师编写的FFT程序比较
1.FFT源程序【kfft.c】
抽样点为：64
CLOCKS_PER_SEC 单位1000ms
#include "math.h"
void kfft(pr,pi,n,k,fr,fi)
int n,k;
double pr[],pi[],fr[],fi[];
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
for (it=0; it<=n-1; it++) //将pr[0]和pi[0]循环赋值给fr[]和fi[]
{
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); //将w=e^-j2pi/n用欧拉公式表示
pi[1]=-sin(p);
for (i=2; i<=n-1; i++) //计算pr[]
{
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;
}
}
for (i=0; i<=n-1; i++)
{
pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]); //幅值计算
}
return;
}
2.源程序【FFT2.c】
输入信号：0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <D:\temp\kfft.c>
#include <time.h>
#define PI 3.1415926
#define N 64
typedef struct
{
double real;
double imag;
}complex;
void Wn(int ,int ,int ,complex *);//定义旋转因子
void c_jiafa(complex ,complex ,complex *); //复数加法运算
void c_chengfa(complex ,complex ,complex *); //复数乘法运算
//void kfft(preal,pimag,n,k,freal,fimag);
int main()
{
complex f[N], A[N], a[N];
double y[N];
int i,j;
double pr[64],pi[64],fr[64],fi[64],t[64];
clock_t begin, end;
double cost1 ,cost2, persent;
for (i=0; i<=63; i++) //生成输入信号
{
t[i] = i*0.001;
pr[i]=0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i); pi[i]=0.0; //0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i)
}
begin = clock(); //开始记录时间
kfft(pr,pi,64,6,fr,fi); //调用FFT函数
end = clock(); //结束记录时间
cost1 = (double)(end - begin) / CLOCKS_PER_SEC;
printf("*****************************FFT变换输出*****************************\n");
for (i=0; i<64; i++)
{
printf("%d\t%lf\n",i,pr[i]); //输出结果
}
printf("*****************************DFT变换输出*****************************\n");
for(int i=0;i<N;i++) //提供初始信号
{
a[i].real=0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i);//0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)
a[i].imag=0;
}
begin = clock(); //开始记录时间
for(int k=0;k<N;k++)
{
//int n;
complex sum={0,0};
for(int n=0; n<N;n++)
{
complex wn,t;
Wn(N,n,k,&wn);
c_chengfa(a[n],wn,&t);
c_jiafa(sum,t,&sum);
}
A[k].real=sum.real;
A[k].imag=sum.imag;
}
end = clock(); //结束记录时间
cost2 = (double)(end - begin) / CLOCKS_PER_SEC;
for(int i=0; i<N;i++)//幅值计算
{
y[i]=sqrt(A[i].real*A[i].real+A[i].imag*A[i].imag);
printf("%d\t%lf\n",i,y[i]);
}
printf("在抽样点为%d的情况下，FFT运算一共花费时间为：%lf秒\n",N, cost1); //CLOCKS_PER_SEC 单位1000ms
printf("在抽样点为%d的情况下，DFT运算一共花费时间为：%lf秒\n",N, cost2); //CLOCKS_PER_SEC 单位1000ms
//persent = (double)(cost1 / cost2);
//printf("在抽样点为%d的情况下，FFT算法运算时间为DFT算法运算时间的%lf倍\n", N, persent);
system("pause");
return 0;
}
void Wn(int n1,int i,int k,complex *Wn) //定义旋转因子
{
//用欧拉公式分成实部和虚部
Wn->real=cos(2*PI*i*k/n1);
Wn->imag=-sin(2*PI*i*k/n1);
}
void c_jiafa(complex a,complex b,complex *c) //复数加法运算
{
c->real=a.real+b.real;
c->imag=a.imag+b.imag;
}
void c_chengfa(complex a,complex b,complex *c) //复数乘法运算
{
c->real=a.real*b.real-a.imag*b.imag;
c->imag=a.real*b.imag+a.imag*b.real;
}
3.TCC下编译
4.结果显示
二、自己理解得FFT程序比较(有不足之处还望为我指出，谢谢大家了)
1.程序源码
输入信号：0.6*sin(2*pi*500*i)+0.6*sin(2*pi*50*i)
抽样点为：64
#include <stdio.h>
#include <math.h>
#include <time.h>
#define N 64 //设置抽样点数
#define PI 3.1514926 //定义圆周率
typedef struct //定义复数结构体变量
{
double real;
double imag;
}complex;
void c_jiafa(complex, complex, complex *); //复数加运算
void c_jianfa(complex, complex, complex *); //复数减运算
void c_chengfa(complex, complex, complex *); //复数乘运算
void Wn_i(int, int, complex *); //FFT旋转因子
void Wn_ik(int, int, int, complex *); //DFT旋转因子
int main()
{
complex f[N], A[N], a[N]; //f[N]为输出的快速傅里叶变换序列 A[N]为DFT变换输出序列 a[N]为初始序列
int LH, K, J, B, L, k, N1, P, M, K1; //L表示第L级蝶形 p旋转因子指数 B两序列间隔点数k 第k个序列
double T;
double y[N];
double pr[N];
clock_t begin, end;
double cost1, cost2, persent; //定义时间
double t_r,t_i;
M =6;
LH = N / 2;
J = LH;
N1 = N - 1;
for (int i = 0; i < N; i++) //为DFT运算提供初始序列
{
a[i].real =0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i) ;
a[i].imag = 0;
}
printf("******************************* 级数为：%d *******************************\n", M);
//.......................................................................................................................
for (int I = 1; I < N1; I++) //定义倒序序列函数
{
if (I < J)
{
t_r=a[I].real;
t_i=a[I].imag;
a[I].real=a[J].real;
a[I].imag=a[J].imag;
a[J].real=t_r;
a[J].imag=t_i;
}
K = LH;
if (J >= K)
{
do
{
J = J - K;
K = K / 2;
} while (J >= K);
}
J = J + K;
}
for (int i = 0; i < N; i++) //将序列赋给结构体函数进行运算
{
f[i].real = 0.6*sin(2*PI*500*i)+0.6*sin(2*PI*50*i);
f[i].imag = 0;
}
begin = clock(); //开始记录时间
for (L = 1; L <= M; L++) //FFT运算
{
B = (int)(pow(2, L - 1));
for (J = 0; J < B; J++)
{
P = (int)(J*pow(2, M - L));
for (k = J; k < N; k = (int)(k + pow(2, L)))
{
K1 = k + B;
complex wn, t;
Wn_i(N, P, &wn);
c_chengfa(f[K1], wn, &t); //。。。。。。。。。。。。
c_jianfa(f[k], t, &(f[K1])); //蝶形运算
c_jiafa(f[k], t, &(f[k])); //。。。。。。。。。。。。
}
}
}
end = clock(); //结束记录时间
cost1 = (double)(end - begin) / CLOCKS_PER_SEC;
//.......................................................................................................................
printf("*****************************快速傅里叶变换输出*****************************\n");
for (int i=0; i<=N-1; i++)
{
pr[i]=sqrt(f[i].real*f[i].real+f[i].imag*f[i].imag); //幅值计算
printf("%d\t%lf\n", i, pr[i] );
}
//.......................................................................................................................
begin = clock(); //开始记录时间
for (int k = 0; k < N; k++) //DFT运算
{
complex sum = { 0, 0 };
for (int n = 0; n < N; n++)
{
complex wn, t;
Wn_ik(N, n, k, &wn);
c_chengfa(a[n], wn, &t);
c_jiafa(sum, t, &sum);
}
A[k].real = sum.real;
A[k].imag = sum.imag;
}
end = clock(); //结束记录时间
cost2 = (double)(end - begin) / CLOCKS_PER_SEC;
//.......................................................................................................................
printf("*****************************DFT变换输出*****************************\n");
for (int i = 0; i < N; i++) //DFT变换输出
{
y[i]=sqrt(A[i].real*A[i].real+A[i].imag*A[i].imag);
printf("%d\t%lf\n",i,y[i]);
//printf("%d\t%lf\n", i, A[i].real, A[i].imag);
}
//.......................................................................................................................
printf("FFT算法一共花费时间为：%lf秒\n", cost1); //CLOCKS_PER_SEC 单位1000ms
printf("DFT算法一共花费时间为：%lf秒\n", cost2); //CLOCKS_PER_SEC 单位1000ms
persent = (float)(cost1 / cost2);
printf("在抽样点为%d的情况下，FFT算法运算时间为DFT算法运算时间的%f倍\n", N, persent);
return 0;
}
void c_jiafa(complex a, complex b, complex *c) //复数加法
{
c->real = a.real + b.real;
c->imag = a.imag + b.imag;
}
void c_jianfa(complex a, complex b, complex *c) //复数减法
{
c->real = a.real - b.real;
c->imag = a.imag - b.imag;
}
void c_chengfa(complex a, complex b, complex *c) //复数乘法
{
c->real = a.real*b.real - a.imag*b.imag;
c->imag = a.real*b.imag + a.imag*b.real;
}
void Wn_i(int n1, int i, complex *Wn) //定义FFT旋转因子
{
Wn->real = cos(2 * PI*i / n1);
Wn->imag = -sin(2 * PI*i / n1);
}
void Wn_ik(int n1, int i, int k, complex *Wn) //定义DFT旋转因子
{
Wn->real = cos(2 * PI*i*k / n1);
Wn->imag = -sin(2 * PI*i*k / n1);
}
2.TCC编译
3.结果显示
三、总结
FFT是DFT的快速算法，可以节省大量的计算时间，快速傅里叶变换（FFT）是一种能在O(nlogn)的时间内将一个多项式转换成它的点值表示的算法
913
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
