下载的几个C语言版的FFT程序仿真移植成功，特分享
/*********************************************************************快速福利叶变换C函数
函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依
赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序的复
数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的
复数
使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的
应该为2的N次方，不满足此条件时应在后面补0
函数调用：FFT(s);
作 者：吉帅虎
时 间：2010-2-20
版 本：Ver1.0
参考文献：
**********************************************************************/
#include<math.h>
#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值
#define FFT_N 128 //定义福利叶变换的点数
struct compx {float real,imag;}; //定义一个复数结构
struct compx s; //FFT输入和输出：从S开始存放，根据大小自己定义
/*******************************************************************
函数原型：struct compx EE(struct compx b1,struct compx b2)
函数功能：对两个复数进行乘法运算
输入参数：两个以联合体定义的复数a,b
输出参数：a和b的乘积，以联合体的形式输出
*******************************************************************/
struct compx EE(struct compx a,struct compx b)
{
struct compx c;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}
/*****************************************************************
函数原型：void FFT(struct compx *xin,int N)
函数功能：对输入的复数组进行快速傅里叶变换（FFT）
输入参数：*xin复数结构体组的首地址指针，struct型
*****************************************************************/
void FFT(struct compx *xin)
{
int f,m,nv2,nm1,i,k,l,j=0;
struct compx u,w,t;
nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nm1=FFT_N-1;
for(i=0;i<nm1;i++)
{
if(i<j) //如果i<j,即进行变址
{
t=xin;
xin=xin;
xin=t;
}
k=nv2; //求j的下一个倒位序
while(k<=j) //如果k<=j,表示j的最高位为1
{
j=j-k; //把最高位变成0
k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k; //把0改为1
}
{
int le,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算
f=FFT_N;
for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数
;
for(m=1;m<=l;m++) // 控制蝶形结级数
{ //m表示第m级蝶形，l为蝶形级总数l=log（2）N
le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点
lei=le/2; //同一蝶形结中参加运算的两点的距离
u.real=1.0; //u为蝶形结运算系数，初始值为1
u.imag=0.0;
w.real=cos(PI/lei); //w为系数商，即当前系数与前一个系数的商
w.imag=-sin(PI/lei);
for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结
{
for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结
{
ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点
t=EE(xin,u); //蝶形运算，详见公式
xin.real=xin.real-t.real;
xin.imag=xin.imag-t.imag;
xin.real=xin.real+t.real;
xin.imag=xin.imag+t.imag;
}
u=EE(u,w); //改变系数，进行下一个蝶形运算
}
}
}
}
/************************************************************
函数原型：void main()
函数功能：测试FFT变换，演示函数使用方法
输入参数：无
输出参数：无
************************************************************/
void main()
{
int i;
for(i=0;i<FFT_N;i++) //给结构体赋值
{
s.real=sin(2*3.141592653589793*i/FFT_N); //实部为正弦波FFT_N点采样，赋值为1
s.imag=0; //虚部为0
}
FFT(s); //进行快速福利叶变换
for(i=0;i<FFT_N;i++) //求变换后结果的模值，存入复数的实部部分
s.real=sqrt(s.real*s.real+s.imag*s.imag);
while(1);
}
下面还有几个版本的。
FFT1.0-1.2ourdev_655640K1EWAP.zip(文件大小:7K) (原文件名:FFT1.0-1.2.zip) 先记一下，以后用的着 mark mark 标记一下 fft有何作用 回复【5楼】SNOOKER 山寨王
fft有何作用
-----------------------------------------------------------------------
时域2频域 mark 顶顶，昨天在AVR那边也看见一个。 mark mark 先记一下，以后用的着 mark mark Mark 看看效果视频就好了 记下 mark 学习 好东西 学习了 其实代码还好，关键是要理解FFT是怎么得来的。干什么用，怎么用。理解了这个，代码就好写了。
我没有具体看，不知道你这个是基几的，多少点的fft。 爪机mark 记下。楼主好人 mark 马马克克 马克 FFT 谢谢了，正需要这个 mark mark mark mark mark 前段时间在找这类的例子 谢谢楼主了 mark 51单片机才刚入门，先mark 先记一下，以后用的着 mark mark mark 不是吧！浮点的，这在51上跑很费力啊。 MAKR c fft mark mark mark I took one of the code and rewrote it slightly to make it more portable.
it has three pieces:
1) fft.h - the header file to be included with any code that utilizes the fft code;
2) fft.c - the source file that provides the actual fft computation
3) main.c - a test file that calls fft.c/fft.h.
all you need to do is to include fft.h in your source file and fft.c in your project.
two customization, all in fft.h:
1) FFT_N: sample size. needs to be power of 2.
2) FFT_TYPE: can be either "double" or "float", depending on your compiler.
the code should be entirely portable. here is fft.h:
#ifndef FFT_H_INCLUDED
#define FFT_H_INCLUDED
/*********************************************************************
快速福利叶变换C函数
函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依
赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序的复
数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的
复数
使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的
应该为2的N次方，不满足此条件时应在后面补0
函数调用：FFT(s);
作 者：吉帅虎
时 间：2010-2-20
版 本：Ver1.0
参考文献：
**********************************************************************/
//fft configuration
#define FFT_N 32 //定义福利叶变换的点数
#define FFT_TYPE double
//end fft configuration
#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值
typedef struct {
FFT_TYPE real,imag;
} compx; //定义一个复数结构
//struct compx s; //FFT输入和输出：从S开始存放，根据大小自己定义
/*****************************************************************
函数原型：void FFT(struct compx *xin,int N)
函数功能：对输入的复数组进行快速傅里叶变换（FFT）
输入参数：*xin复数结构体组的首地址指针，struct型
*****************************************************************/
void fft(compx *xin);
#endif // FFT_H_INCLUDED here is fft.c:
/*********************************************************************
快速福利叶变换C函数
函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依
赖硬件。此函数采用联合体的形式表示一个复数，输入为自然顺序的复
数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的
复数
使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的
应该为2的N次方，不满足此条件时应在后面补0
函数调用：FFT(s);
作 者：吉帅虎
时 间：2010-2-20
版 本：Ver1.0
参考文献：
**********************************************************************/
#include <math.h>
#include "fft.h"
/*******************************************************************
函数原型：struct compx EE(struct compx b1,struct compx b2)
函数功能：对两个复数进行乘法运算
输入参数：两个以联合体定义的复数a,b
输出参数：a和b的乘积，以联合体的形式输出
*******************************************************************/
compx EE(compx a, compx b) {
compx c;
c.real=a.real*b.real-a.imag*b.imag;
c.imag=a.real*b.imag+a.imag*b.real;
return(c);
}
/*****************************************************************
函数原型：void FFT(struct compx *xin,int N)
函数功能：对输入的复数组进行快速傅里叶变换（FFT）
输入参数：*xin复数结构体组的首地址指针，struct型
*****************************************************************/
void fft(compx *xin) {
int f,m,nv2,nm1,i,k,l,j=0;
compx u,w,t;
nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nm1=FFT_N-1;
for(i=0; i<nm1; i++) {
if(i<j) { //如果i<j,即进行变址
t=xin;
xin=xin;
xin=t;
}
k=nv2; //求j的下一个倒位序
while(k<=j) { //如果k<=j,表示j的最高位为1
j=j-k; //把最高位变成0
k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
}
j=j+k; //把0改为1
}
{
int le,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算
f=FFT_N;
for(l=1; (f=f/2)!=1; l++) //计算l的值，即计算蝶形级数
;
for(m=1; m<=l; m++) { // 控制蝶形结级数
//m表示第m级蝶形，l为蝶形级总数l=log（2）N
le=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点
lei=le/2; //同一蝶形结中参加运算的两点的距离
u.real=1.0; //u为蝶形结运算系数，初始值为1
u.imag=0.0;
w.real=cos(PI/lei); //w为系数商，即当前系数与前一个系数的商
w.imag=-sin(PI/lei);
for(j=0; j<=lei-1; j++) { //控制计算不同种蝶形结，即计算系数不同的蝶形结
for(i=j; i<=FFT_N-1; i=i+le) { //控制同一蝶形结运算，即计算系数相同蝶形结
ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点
t=EE(xin,u); //蝶形运算，详见公式
xin.real=xin.real-t.real;
xin.imag=xin.imag-t.imag;
xin.real=xin.real+t.real;
xin.imag=xin.imag+t.imag;
}
u=EE(u,w); //改变系数，进行下一个蝶形运算
}
}
}
}
/************************************************************
函数原型：void main()
函数功能：测试FFT变换，演示函数使用方法
输入参数：无
输出参数：无
************************************************************/
/*
int main(void) {
int i;
//generate the data
for(i=0; i<FFT_N; i++) { //给结构体赋值
s.real = 1.0+ \
2.0*cos(2*PI*1*(3.6/360+1.0*i/FFT_N))+ \
3.0*cos(2*PI*4*(1.8/360+1.0*i/FFT_N))+ \
4.0*cos(2*PI*6*(1.8/360+1.0*i/FFT_N))+ \
1.0e-1*rand()/RAND_MAX; //实部为正弦波FFT_N点采样，赋值为1
s.imag=0; //虚部为0
}
//print out the time-domain data
printf("time-domain data...\n");
for (i=0; i<FFT_N; i++)
printf("i=%4d, real = %10.7f, imag = %10.7f, amp= %10.7f\n", i, s.real, s.imag, sqrt(s.real*s.real + s.imag*s.imag));
//perform fft
fft(s); //进行快速福利叶变换
//printf the fft data
printf("post-FFT data...\n");
for (i=0; i<FFT_N/2; i++)
printf("i=%4d, real = %10.7f, imag = %10.7f, amp= %10.2f\n", i, s.real/FFT_N, s.imag/FFT_N, sqrt(s.real*s.real + s.imag*s.imag)/((i==0)?FFT_N:(FFT_N/2)));
//for(i=0; i<FFT_N; i++) //求变换后结果的模值，存入复数的实部部分
//s.real=sqrt(s.real*s.real+s.imag*s.imag);
//printf the fft data
//printf("mod-FFT data...\n");
//for (i=0; i<FFT_N; i++)
//printf("i=%4d, amp = %10.7f\n", i, s.real, s.imag);
//while(1);
return 0;
}
*/ here is main.c (the test program):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "fft.h"
compx s; //FFT输入和输出：从S开始存放，根据大小自己定义
/************************************************************
函数原型：void main()
函数功能：测试FFT变换，演示函数使用方法
输入参数：无
输出参数：无
************************************************************/
int main(void) {
int i;
//generate the data
for(i=0; i<FFT_N; i++) { //给结构体赋值
s.real = 1.0+ \
2.0*cos(2*PI*1*(3.6/360+1.0*i/FFT_N))+ \
3.0*cos(2*PI*4*(1.8/360+1.0*i/FFT_N))+ \
4.0*cos(2*PI*13*(1.8/360+1.0*i/FFT_N))+ \
1.0e-1*rand()/RAND_MAX; //实部为正弦波FFT_N点采样，赋值为1
s.imag= 1.0e-1*rand()/RAND_MAX; //虚部为0
}
//print out the time-domain data
printf("time-domain data...\n");
for (i=0; i<FFT_N; i++)
printf("i=%4d, real = %10.3f, imag = %10.3f, amp= %10.3f\n", i, s.real, s.imag, sqrt(s.real*s.real + s.imag*s.imag));
//perform fft
fft(s); //进行快速福利叶变换
//printf the fft data
printf("post-FFT data...\n");
for (i=0; i<FFT_N/2; i++)
printf("i=%4d, real = %10.3f, imag = %10.3f, amp= %10.3f\n", i, s.real/FFT_N, s.imag/FFT_N, sqrt(s.real*s.real + s.imag*s.imag)/((i==0)?FFT_N:(FFT_N/2)));
//for(i=0; i<FFT_N; i++) //求变换后结果的模值，存入复数的实部部分
//s.real=sqrt(s.real*s.real+s.imag*s.imag);
//printf the fft data
//printf("mod-FFT data...\n");
//for (i=0; i<FFT_N; i++)
//printf("i=%4d, amp = %10.3f\n", i, s.real, s.imag);
//while(1);
return 0;
} what main.c does is to generate a simulated series of sample data -
s.real = 1.0+ \
2.0*cos(2*PI*1*(3.6/360+1.0*i/FFT_N))+ \
3.0*cos(2*PI*4*(1.8/360+1.0*i/FFT_N))+ \
4.0*cos(2*PI*13*(1.8/360+1.0*i/FFT_N))+ \
1.0e-1*rand()/RAND_MAX; //实部为正弦波FFT_N点采样，赋值为1
s.imag= 1.0e-1*rand()/RAND_MAX; //虚部为0
so s.real has a dc of 1, 1*fs at amptitude of 2.0, 4*fs at amplitude of 3.0, 13*fs at amplitude of 4.0, and a random factor (to simulate noise).
what we expect from the output would be this:
0fs (dc): 1
1fs: 2.0
2fs: 0.0
3fs: 0.0
4fs: 3.0
...
13fs:4.0
...
you can run your code and you will see the output in the "post-fft data..." section.
in my case, this is what I got:
0fs (dc): 1.049 vs. 1.0
1fs: 1.985 vs. 2.0
2fs: 0.007 vs. 0.0
3fs: 0.016 vs. 0.0
4fs: 2.998 vs. 3.0
...
13fs:3.999 vs. 4.0
fairly accurate. because of its use of floating point math, the code is very slow: 32 points fft takes 170k cycles, and 128 point fft takes 800k cycles. the code compiled and tested on turbo c, but should work on any other c compilers as well. the problem with floating point math on a 8-bit mcu is its size / speed. so I adapted a fixed point fft routine.
it completes 32-point fft in 20k cycles and 128-point fft in 107k cycles (on a 24f chip, no optimization and no hardware acceleration) - a ~10x improvement over the floating point math version.
again, the package consists of three files:
1) fft_fp.h: the header file. you need to configure again the buffer size + data type.
2) fft_fp.c: the source file. you need to include it in your project.
3) main.c: the test program.
the calling convention is (almost) identical to that of the floating math version. here is the package:
点击此处下载 ourdev_716318BSJNK0.rar(文件大小:9K) (原文件名:24f fft_fp.rar) I am on a roll!
here is a kalman filter, for 1d variables (aka only one state variable). It is useful for filtering, and or tracking.
the package consists of two files:
1) kalman.h: the header file declaring the functions;
2) kalman.c: the source file laying out the actual code / execution.
in additional, I put together a main.c calling both functions to show how it could be used.
here is the kalman.h:
#ifndef KALMAN_H_INCLUDED
#define KALMAN_H_INCLUDED
//kalman filtering for 1d variables
//source: http://interactive-matter.eu/2009/12/filtering-sensor-data-with-a-kalman-filter/
typedef struct {
double q; //process noise covariance
double r; //measurement noise covariance
double x; //value
double p; //estimation error covariance
double k; //kalman gain
} kalman_state;
//initialize the kalman filter
kalman_state kalman_init(double q, double r, double p, double intial_value);
void kalman_update(kalman_state* state, double measurement);
#endif // KALMAN_H_INCLUDED here is kalman.c:
//kalman filtering for 1d variables
//source: http://interactive-matter.eu/2009/12/filtering-sensor-data-with-a-kalman-filter/
#include "kalman.h" //we use kalman filter
kalman_state kalman_init(double q, double r, double p, double intial_value)
{
kalman_state result;
result.q = q;
result.r = r;
result.p = p;
result.x = intial_value;
return result;
}
void kalman_update(kalman_state* state, double measurement)
{
//prediction update
//omit x = x
state->p = state->p + state->q;
//measurement update
state->k = state->p / (state->p + state->r);
state->x = state->x + state->k * (measurement - state->x);
state->p = (1 - state->k) * state->p;
} here is main.c: it uses a function to generate sampled data.
#include <stdio.h>
#include <stdlib.h>
#include "kalman.h" //we use kalman filtering
//generate simulated data, with a noise term
double x_actual(void) {
static double xt_1=-10, xt_2=20;
double xt = 1.2*xt_1 - 1.0*xt_2 + 00.0 * rand()/RAND_MAX;
xt_2 = xt_1;
xt_1 = xt;
return xt;
}
int main()
{
FILE *fptr; //file ptr
int i;
kalman_state k_state;
double xt;
//open the file
fptr = fopen("c:\kalman.dat", "w");
//fprintf(fptr, "abcdef\n");
k_state = kalman_init(1, 1, 1, 1); //initialize the filter
printf("Hello world!\n");
for (i=0; i<200; i++) {
xt = x_actual(); //generate xt's measurement
kalman_update(&k_state, xt); //estimate xt
fprintf(stdout, "i=%4d x_actual = %10.4f x_est = %10.4f\n", i, xt, k_state.x);
}
fclose(fptr); //close the file
return 0;
} in this particular implementation, you can declare multiple kalman_state variables and run them on the same filter without compromising the data. This allows, for example, tracking of multiple noisy variables using the same code base, particularly helpful for multi-dimensional applications, like flight control based on a 3d sensor. mark 数字信号处理，MARK 很不错学习资料。
快速福利叶变换C函数 mark 果断mark 学习，学习！！！ mark 51的不用考虑做FFT了 biaojiyixia mark{:handshake:} mark！ FFT 学习MARK 大家难道没有看代码么，写得这么粗糙 代码有很多校错误的 主体部分看上去还行 学习了， MARK 下载先存着 mark with FFT function 赞一个，呵呵 好东西啊 谢谢分享 谢谢楼主，用过st例程的 mark! 兼挖 mark,谢谢楼主分享！ 谢谢楼主
标记一下，谢谢楼主 记号一个备用 xin=xin;
这行编译不过！！
FFT.C(56): error C193: '=': incompatible operand 谢谢，楼主。
页:
[1]