📄 fft快速傅立叶变换程序.txt
字号:
fft快速傅立叶变换程序关键词： 程序
实验用的头文件 MYFFT.H
作用：为帮助小虎子做实验，这个头文件提供了完整的一维与二维FFT算法，我想应改是够你折腾了吧！
#include <complex> // complex<float>
using namespace std;
typedef complex<float> Comp; // 复数类型定义
const float _2PI_ = 2.0f * 3.14159265f; // 常数2PI定义
const int MAX_N = 256; // 最大DFT点数
/*----*----*----*----*----*----*----*----*----*----*----*----*
FFT算法模块接口定义
*----*----*----*----*----*----*----*----*----*----*----*----*/
///////////////////////////////////////////
// Function name : BitReverse
// Description : 二进制倒序操作
// Return type : int
// Argument : int src 待倒读的数
// Argument : int size 二进制位数
int BitReverse(int src, int size)
{
int tmp = src;
int des = 0;
for (int i=size-1; i>=0; i--)
{
des = ((tmp & 0x1) << i) | des;
tmp = tmp >> 1;
}
return des;
}
//////////////////////////////////////////////////
// Function name : Reorder
// Description : 数据二进制整序
// Return type : void
// Argument : Comp x[MAX_N] 待整序数组
// Argument : int N FFT点数
// Argument : int M 点数的2的幂次
void Reorder(Comp x[MAX_N], int N, int M)
{
Comp new_x[MAX_N];
for (int i=0; i<N; i++)
new_x = x[BitReverse(i, M)];
// 重新存入原数据中(已经是二进制整序过了的数据)
for (i=0; i<N; i++)
x = new_x;
}
//////////////////////////////////////////////////
// Function name : CalcW
// Description : 计算旋转因子
// Return type : void
// Argument : Comp W[MAX_N] 存放因子的数组
// Argument : int N FFT的点数
// Argument : int flag 正逆变换标记
void CalcW(Comp W[MAX_N], int N, int flag)
{
for (int i=0; i<N/2; i++)
{
if (flag == 1)
W = exp(Comp(0, -_2PI_ * i / N)); // 正FFT变换
else
W = exp(Comp(0, _2PI_ * i / N)); // 逆FFT变换
}
}
/////////////////////////////////////////////////
// Function name : FFT_1D_Kernel
// Description : FFT算法核心
// Return type : void
// Argument : Comp* x 数据
// Argument : int M 幂次
// Argument : int flag 正逆变换标记
以下本应由自己完成。
void FFT_1D(Comp* x, int M, int flag)
{
int N = (1 << M);
// 二进制整序
Reorder(x, N, M);
// 旋转因子计算
Comp W[MAX_N];
CalcW(W, N, flag);
// 级内群数
int GroupNum = N/2; // 第一级的群数为N/2
// 群内蝶形单元数
int CellNum = 1; // 第一级的群内蝶形单元数为1
// 处理各级
for (int i=0; i<M; i++)
{
// 处理各群
for (int j=0; j<GroupNum; j++)
{
// 处理各蝶形单元
for (int k=0; k<CellNum; k++)
{
// (1) 计算出当前蝶形单元所含元素在数据数组中的位置
// 第一元素位置
int Pos1 = CellNum * j * 2 + k ; // 已经处理了前 j 群，每群有 CellNum 个单元，
每单元有 2 个元素
// 第二元素位置
int Pos2 = Pos1 + CellNum;
// (2) 计算旋转因子与单元的第二元素的复数乘积
Comp TMP = x[Pos2] * W[k * GroupNum] ;
// (3) 计算最终结果, 并存入到数组的原先位置
x[Pos2] = x[Pos1] - TMP ;
x[Pos1] = x[Pos1] + TMP ;
}
}
GroupNum >>= 1; // 级别增加, 则相应的群数减少一半
CellNum <<= 1; // 级别增加, 则相应的群内单元数增加一倍
}
// 如果是IFFT，各元素还要再除以N
if (flag != 1)
{
for (i=0; i<N; i++)
x /= N;
}
}
//////////////////////////////////////////////////////
// Function name : FFT_2D_Kernel
// Description : 2D FFT核心算法
// Return type : void
// Argument : Comp x[MAX_N][MAX_N] 二维数据
// Argument : int M 幂次
// Argument : int flag 正逆变换标记
以下本应由自己完成。
void FFT_2D(Comp x[MAX_N][MAX_N], int M, int flag)
{
int N = (1 << M);
// 先逐行进行 1D-FFT
for (int i=0; i<N; i++)
FFT_1D(x, M, flag); // <--- 计算结果再存入矩阵x中
// 再逐列进行 1D-FFT
Comp col[MAX_N];
for (int j=0; j<N; j++)
{
// 取得第j列的数据
for (int i=0; i<N; i++)
col = x[j];
// 对第j列数据进行 1D-FFT
FFT_1D(col, M, flag); // <--- 计算结果在数组col中
// 将结果放回矩阵第j列中
for (i=0; i<N; i++)
x[j] = col;
}
}
// <--- End of [FFT.H]
快速傅立叶变换（FFT）C语言源程序
Tags: FFT 快速傅立叶变换 C语言 程序 1年2月前 ，6919 次点击
作者: linker110
语言: 英语
快速傅立叶变化的C语言程序，比起MATLAB里简单的一句话可是差得远了，但是这个东西我可是花了很多的心血啊，虽然做的不怎么。
里面附带了试验的图片，相信学了数字信号的同学都会知道是什么东西了。
对C语言源文件：
1.complex.h定义了复数的加减乘除运算的子程序，用于傅立叶变换
2.dsp.h定义了傅立叶变换及快速傅立叶变换的正反变换，和相关的逆序程序。但本实验并未使用DFT
3.savebmp.h定义了图形模式下保存屏幕的函数
4.main.c为相应的主程序。其中的#include 均使用<>格式
5.进行快速傅立叶变换是输入数据必须合法，只能是2的n次幂，且未做错误处理
6.对虚部部为零的序列进行离散傅立叶变换是相频特性似乎有时会出现错误，尚未解决
7.画图的系列函数没有做成独立的文件所以，所以main函数过于杂乱
8.对复数信号在作图时以顺序自然数做横轴以复数的实数部分做纵轴，并不合理
9.画图部分尚未实现autoscale,是将原点定在画布中心，对最大绝对值进行归一化。坐标原点可变，但并不提倡
10.保存BPM文件的源代码来自网络
11.对大部分判断正确执行的返回值未作处理。
特别：SAVEBMP.H主题来在网络，我只是把它改动了一下成为一个相对独立可调用的模块。
存在一些问题，欢迎交流^_^
还有，<free>message</free>写在哪里才有效呢？
*****************fft programe*********************/
#include "typedef.h"
#include "math.h"
struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}
void FFT(struct compx *xin,int N)
{
int f,m,nv2,nm1,i,k,j=1,l;
/*int f,m,nv2,nm1,i,k,j=N/2,l;*/
struct compx v,w,t;
nv2=N/2;
f=N;
for(m=1;(f=f/2)!=1;m++){;}
nm1=N-1;
/*变址运算*/
for(i=1;i<=nm1;i++)
{
if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=nv2;
while(k<j){j=j-k;k=k/2;}
j=j+k;
}
{
int le,lei,ip;
float pi;
for(l=1;l<=m;l++)
{ le=pow(2,l);// 这里用的是L而不是1 !!!!
lei =le/2;
pi=3.14159;
v.real=1.0;
v.imag=0.0;
w.real=cos(pi/lei);
w.imag=-sin(pi/lei);
for(j=1;j<=lei;j++)
{
/*double p=pow(2,m-l)*j;
double ps=2*pi/N*p;
w.real=cos(ps);
w.imag=-sin(ps);*/
for(i=j;i<=N;i=i+le)
{ /* w.real=cos(ps);
w.imag=-sin(ps);*/
ip=i+lei;
t=EE(xin[ip],v);
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
v=EE(v,w);
}
}
}
return;
}
/*****************main programe********************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "typedef.h"
float result[257];
struct compx s[257];
int Num=256;
const float pp=3.14159;
main()
{
int i=1;
for(;i<0x101;i++)
{
s[i].real=sin(pp*i/32);
s[i].imag=0;
}
FFT(s,Num);
for(i=1;i<0x101;i++)
{
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
}
}
(文章推荐人：丁丁)
fft程序（c语言）收藏
void fftt(struct COMPLEX *a,int l)
{
const double pai=3.141592653589793;
struct COMPLEX u,w,t,tmp1;
unsigned n=1,nv2,nm1,k,le,lei,ip;
unsigned i,j,m;
double tmp;
n<<=l;
nv2=n>>1;
nm1=n-1;
j=0;
for(i=0;i<nm1;i++)
{
if(i<j)
{
t.re=a[j].re;
t.im=a[j].im;
a[j].re=a[i].re;
a[j].im=a[i].im;
a[i].re=t.re;
a[i].im=t.im;
}
k=nv2;
while(k<=j)
{
j-=k;
k>>=1;
}
j+=k;
}
le=1;
for(m=1;m<=l;m++)
{
lei=le;
le<<=1;
u.re=1.0;
u.im=0.0;
tmp=pai/lei;
w.re=cos(tmp);
w.im=-sin(tmp);
for(j=0;j<lei;j++)
{
for(i=j;i<n;i+=le)
{
ip=i+lei;
mul(a[ip],u,&tmp1);
t.re=tmp1.re;
t.im=tmp1.im;
sub(a[i],t,&tmp1);
a[ip].re=tmp1.re;
a[ip].im=tmp1.im;
add(a[i],t,&tmp1);
a[i].re=tmp1.re;
a[i].im=tmp1.im;
}
mul(u,w,&tmp1);
u.re=tmp1.re;
u.im=tmp1.im;
}
}
}
发表于 @ 2004年10月09日 15:28:00|评论(loading...)|编辑
|
快速傅立叶变换程序（转）
magicchip 发表于 2008-11-27 5:55:00 0
推荐＃i nclude <complex> // complex<float>
using namespace std;
typedef complex<float> Comp; // 复数类型定义
const float _2PI_ = 2.0f * 3.14159265f; // 常数2PI定义
const int MAX_N = 256; // 最大DFT点数
/*----*----*----*----*----*----*----*----*----*----*----*----*
FFT算法模块接口定义
*----*----*----*----*----*----*----*----*----*----*----*----*/
///////////////////////////////////////////
// name : BitReverse
// Deion : 二进制倒序操作
// Return type : int
// Argument : int src 待倒读的数
// Argument : int size 二进制位数
int BitReverse(int src, int size)
{
int tmp = src;
int des = 0;
for (int i=size-1; i>=0; i--)
{
des = ((tmp & 0x1) << i) | des;
tmp = tmp >> 1;
}
return des;
}
//////////////////////////////////////////////////
// name : Reorder
// Deion : 数据二进制整序
// Return type : void
// Argument : Comp x[MAX_N] 待整序数组
// Argument : int N FFT点数
// Argument : int M 点数的2的幂次
void Reorder(Comp x[MAX_N], int N, int M)
{
Comp new_x[MAX_N];
for (int i=0; i<N; i++)
new_x = x[BitReverse(i, M)];
// 重新存入原数据中(已经是二进制整序过了的数据)
for (i=0; i<N; i++)
x = new_x;
}
//////////////////////////////////////////////////
// name : CalcW
// Deion : 计算旋转因子
// Return type : void
// Argument : Comp W[MAX_N] 存放因子的数组
// Argument : int N FFT的点数
// Argument : int flag 正逆变换标记
void CalcW(Comp W[MAX_N], int N, int flag)
{
for (int i=0; i<N/2; i++)
{
if (flag == 1)
W = exp(Comp(0, -_2PI_ * i / N)); // 正FFT变换
else
W = exp(Comp(0, _2PI_ * i / N)); // 逆FFT变换
}
}
/////////////////////////////////////////////////
// name : FFT_1D_Kernel
// Deion : FFT算法核心
// Return type : void
// Argument : Comp* x 数据
// Argument : int M 幂次
// Argument : int flag 正逆变换标记
void FFT_1D(Comp* x, int M, int flag)
{
int N = (1 << M);
// 二进制整序
Reorder(x, N, M);
// 旋转因子计算
Comp W[MAX_N];
CalcW(W, N, flag);
// 级内群数
int GroupNum = N/2; // 第一级的群数为N/2
// 群内蝶形单元数
int CellNum = 1; // 第一级的群内蝶形单元数为1
// 处理各级
for (int i=0; i<M; i++)
{
// 处理各群
for (int j=0; j<GroupNum; j++)
{
// 处理各蝶形单元
for (int k=0; k<CellNum; k++)
{
// (1) 计算出当前蝶形单元所含元素在数据数组中的位置
// 第一元素位置
int Pos1 = CellNum * j * 2 + k ; // 已经处理了前 j 群，每群有 CellNum 个单元， 每单元有 2 个元素
// 第二元素位置
int Pos2 = Pos1 + CellNum;
// (2) 计算旋转因子与单元的第二元素的复数乘积
Comp TMP = x[Pos2] * W[k * GroupNum] ;
// (3) 计算最终结果, 并存入到数组的原先位置
x[Pos2] = x[Pos1] - TMP ;
x[Pos1] = x[Pos1] + TMP ;
}
}
GroupNum >>= 1; // 级别增加, 则相应的群数减少一半
CellNum <<= 1; // 级别增加, 则相应的群内单元数增加一倍
}
// 如果是IFFT，各元素还要再除以N
if (flag != 1)
{
for (i=0; i<N; i++)
x /= N;
}
}
//////////////////////////////////////////////////////
// name : FFT_2D_Kernel
// Deion : 2D FFT核心算法
// Return type : void
// Argument : Comp x[MAX_N][MAX_N] 二维数据
// Argument : int M 幂次
// Argument : int flag 正逆变换标记
void FFT_2D(Comp x[MAX_N][MAX_N], int M, int flag)
{
int N = (1 << M);
// 先逐行进行 1D-FFT
for (int i=0; i<N; i++)
FFT_1D(x, M, flag); // <--- 计算结果再存入矩阵x中
// 再逐列进行 1D-FFT
Comp col[MAX_N];
for (int j=0; j<N; j++)
{
// 取得第j列的数据
for (int i=0; i<N; i++)
col = x[j];
// 对第j列数据进行 1D-FFT
FFT_1D(col, M, flag); // <--- 计算结果在数组col中
// 将结果放回矩阵第j列中
for (i=0; i<N; i++)
x[j] = col;
}
}
// <--- End of [FFT.H]
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -