发信人: soo (狙击手), 信区: HITEA
标 题: 快速付立叶变换c程序
发信站: 哈工大紫丁香 (2001年05月09日19:19:37 星期三), 转信
【 以下文字转载自 soo 的信箱 】
【 原文由 ttq 所发表 】
快速傅立叶变换
程序说明
用C语言自编一个基二的实FFT变换程序，输入需变换的实数，输出变换后的实部和虚部
．FFT程序为
void myfft(double * mytime,int length,complex * spec)
mytime为需变换的时间信号．length为信号的长度，一般取为二的幂次方，本程序中其
值应不大于1024．spec为变换后的结果，为复数结构．
complex为自定义的结构，其中的real代表其实部，im代表其虚部．
c_mul , c_add , c_sub分别为复数乘，复数加，复数减子程序．
complex w_exp(int up,int down) 为计算旋转因子Wkn的子程序,k即程序中的up , n即
down．
程序编号之后，在LabWindows/CVI调用myfft程序，分别对长度为1024,采样频率为1024
Hz的正弦信号和三角波信号进行FFT变换，然后把所得结果转化为极坐标的形式，并乘二
后在示波器上显示出来（显示单边幅值，故应乘二）．
程序
#include <analysis.h>
#include <fftuir.h>
#include <cvirte.h> /* Needed if linking in external compiler; harmless othe
rwise */
#include <userint.h>
#include <ansi_c.h>
static int panelHandle;
int k=1,i;
double pi=3.1415926;
double temp[1024],mytime[1024];
typedef struct _complex
{
double real;
double im;
}complex;
complex ctemp[1024];
complex c_mul(complex c1,complex c2) ;
complex c_add(complex c1,complex c2) ;
complex c_sub(complex c1,complex c2) ;
complex w_exp(int up,int down) ;
//FFT快速算法
void myfft(double * mytime,int length,complex * spec)
{
int i,j,l,k;
int con,con1;
int con,con1;
//排序
for(k=1;pow(2,k)<=length;k++)
{
con=(int)length/pow(2,(k-1));
for (j=0;j<pow(2,(k-1));j++)
{
for(l=0;l<con/2;l++)
{
temp[l+j*con]=mytime[j*con+2*l];
temp[l+j*con+con/2]=mytime[j*con+2*l+1];
}
}
for(i=0;i<length;i++)
{
mytime[i]=temp[i];
}
}
//实数转化为复数
for(i=0;i<length;i++)
{
spec[i].real=mytime[i];
spec[i].im=0;
spec[i].im=0;
}
//进行运算, 乘旋转因子
for(k=1;pow(2,k)<=length;k++)
{
con1=pow(2,k);
for(j=0;j<length/con1;j++)
{
for(l=0;l<con1/2;l++)
{
ctemp[l+j*con1]=c_add (spec[j*con1+l],
c_mul(spec[j*con1+l+con1/2],w_exp(l,con1)));
ctemp[l+con1/2+j*con1]=c_sub (spec[j*con1+l],
c_mul(spec[j*con1+l+con1/2],w_exp(l,con1)));
}
}
for(i=0;i<length;i++)
{
spec[i]=ctemp[i];
}
}
return;
}
}
//复数乘
complex c_mul(complex c1,complex c2)
{
complex result;
result.real=c1.real*c2.real-c1.im*c2.im;
result.im=c1.real*c2.im+c1.im*c2.real;
return result;
}
//复数加
complex c_add(complex c1,complex c2)
{
complex result;
result.real=c1.real+c2.real;
result.im=c1.im+c2.im;
return result;
}
//复数减
complex c_sub(complex c1,complex c2)
{
complex result;
result.real=c1.real-c2.real;
result.im=c1.im-c2.im;
return result;
}
//旋转因子
complex w_exp(int up,int down)
{
complex result;
result.real=cos(-2*pi*up/down);
result.im=sin(-2*pi*up/down) ;
return result;
}
//自编程序的应用
int main (int argc, char *argv[])
{
if (InitCVIRTE (0, argv, 0) == 0)
return -1; /* out of memory */
if ((panelHandle = LoadPanel (0, "fftuir.uir", PANEL)) < 0)
return -1;
DisplayPanel (panelHandle);
RunUserInterface ();
return 0;
}
int CVICALLBACK SIGNAL (int panel, int control, int event,
void *callbackData, int eventData1, int eventData2)
{
int signal_type;
double time[1024],amp[512];
double frequence;
double phase=0;
complex specture[1024];
switch (event)
{
case EVENT_COMMIT:
GetCtrlVal (panelHandle, PANEL_SIGNAL, &signal_type);
GetCtrlVal (panelHandle, PANEL_FREQUENCE, &frequence);
switch (signal_type)
{
case 1:
SinePattern (1024, 1.0, 0.0, frequence, time); //产生正弦波
break;
case 2:
TriangleWave (1024, 1.0, frequence/1024, &phase, time);//产生三角波
break;
default:
default:
return 0;
}
DeleteGraphPlot (panelHandle, PANEL_TIME, -1, VAL_IMMEDIATE_DRAW);
SetCtrlAttribute (panelHandle, PANEL_TIME, ATTR_XAXIS_GAIN, 1.0/1024);
PlotY (panelHandle, PANEL_TIME, time, 128, VAL_DOUBLE,
VAL_THIN_LINE,VAL_EMPTY_SQUARE,VAL_SOLID,1,VAL_RED); //画时域波形
myfft(&time[0],1024,specture); //进行FFT变换
amp[0]=sqrt((specture[0].real*specture[0].real+
specture[0].im*specture[0].im)/(1024*1024)); //转化成极坐标
for(i=1;i<512;i++)
{
amp[i]=2*sqrt(((specture[i].real*specture[i].real)+
(specture[i].im*specture[i].im))/(1024*1024));
}
DeleteGraphPlot ( panelHandle,PANEL_SPECTURE,-1,
VAL_IMMEDIATE_DRAW);
PlotY (panelHandle, PANEL_SPECTURE, amp, 512, VAL_DOUBLE,
VAL_THIN_LINE,
VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_RED); //画频域波形
break;
}
return 0;
return 0;
}
//程序终止
int CVICALLBACK EXIT (int panel, int control, int event,
void *callbackData, int eventData1, int eventData2)
{
switch (event)
{
case EVENT_COMMIT:
QuitUserInterface (0);
break;
}
return 0;
}
--
※ 来源:·哈工大紫丁香 bbs.hit.edu.cn·[FROM: 202.118.235.173]
--
※ 转载:．哈工大紫丁香 bbs.hit.edu.cn．[FROM: xmz621.hit.edu.cn]