求FFT的c语言程序

快速傅里叶变换 要用C++ 才行吧 你可以用MATLAB来实现更方便点啊

此FFT 是用VC6.0编写，由FFT.CPP；STDAFX.H和STDAFX.CPP三个文件组成，编译成功。程序可以用文件输入和输出为文件。文件格式为TXT文件。测试结果如下：

输入文件：8.TXT 或手动输入

8 //N

1

2

3

4

5

6

7

8

输出结果为：或保存为TXT文件。（8OUT.TXT）

8

(36,0)

(-4,9.65685)

(-4,4)

(-4,1.65685)

(-4,0)

(-4,-1.65685)

(-4,-4)

(-4,-9.65685)

下面为FFT.CPP文件：

// FFT.cpp : 定义控制台应用程序的入口点。

#include "stdafx.h"

#include <iostream>

#include <complex>

#include <bitset>

#include <vector>

#include <conio.h>

#include <string>

#include <fstream>

using namespace std

bool inputData(unsigned long &, vector<complex<double>>&) //手工输入数据

void FFT(unsigned long &, vector<complex<double>>&) //FFT变换

void display(unsigned long &, vector<complex<double>>&) //显示结果

bool readDataFromFile(unsigned long &, vector<complex<double>>&) //从文件中读取数据

bool saveResultToFile(unsigned long &, vector<complex<double>>&) //保存结果至文件中

const double PI = 3.1415926

int _tmain(int argc, _TCHAR* argv[])

{

vector<complex<double>>vecList //有限长序列

unsigned long ulN = 0 //N

char chChoose = ' ' //功能选择

//功能循环

while(chChoose != 'Q' &&chChoose != 'q')

{

//显示选择项

cout <<"\nPlease chose a function" <<endl

cout <<"\t1.Input data manually, press 'M':" <<endl

cout <<"\t2.Read data from file, press 'F':" <<endl

cout <<"\t3.Quit, press 'Q'" <<endl

cout <<"Please chose:"

//输入选择

chChoose = getch()

//判断

switch(chChoose)

{

case 'm': //手工输入数据

case 'M':

if(inputData(ulN, vecList))

{

FFT(ulN, vecList)

display(ulN, vecList)

saveResultToFile(ulN, vecList)

}

break

case 'f'://从文档读取数据

case 'F':

if(readDataFromFile(ulN, vecList))

{

FFT(ulN, vecList)

display(ulN, vecList)

saveResultToFile(ulN, vecList)

}

break

}

}

return 0

}

bool Is2Power(unsigned long ul) //判断是否是2的整数次幂

{

if(ul <2)

return false

while( ul >1 )

{

if( ul % 2 )

return false

ul /= 2

}

return true

}

bool inputData(unsigned long &ulN, vector<complex<double>>&vecList)

{

//题目

cout<<"\n\n\n==============================Input Data===============================" <<endl

//输入N

cout<<"\nInput N:"

cin>>ulN

if(!Is2Power(ulN)) //验证N的有效性

{

cout<<"N is invalid (N must like 2, 4, 8, .....), please retry." <<endl

return false

}

//输入各元素

vecList.clear() //清空原有序列

complex<double>c

for(unsigned long i = 0i <ulNi++)

{

cout <<"Input x(" <<i <<"):"

cin >> c

vecList.push_back(c)

}

return true

}

bool readDataFromFile(unsigned long &ulN, vector<complex<double>>&vecList) //从文件中读取数据

{

//题目

cout<<"\n\n\n===============Read Data From File==============" <<endl

//输入文件名

string strfilename

cout <<"Input filename:"

cin >>strfilename

//打开文件

cout <<"open file " <<strfilename <<"......." <<endl

ifstream loadfile

loadfile.open(strfilename.c_str())

if(!loadfile)

{

cout <<"\tfailed" <<endl

return false

}

else

{

cout <<"\tsucceed" <<endl

}

vecList.clear()

//读取N

loadfile >>ulN

if(!loadfile)

{

cout <<"can't get N" <<endl

return false

}

else

{

cout <<"N = " <<ulN <<endl

}

//读取元素

complex<double>c

for(unsigned long i = 0i <ulNi++)

{

loadfile >>c

if(!loadfile)

{

cout <<"can't get enough infomation" <<endl

return false

}

else

cout <<"x(" <<i <<") = " <<c <<endl

vecList.push_back(c)

}

//关闭文件

loadfile.close()

return true

}

bool saveResultToFile(unsigned long &ulN, vector<complex<double>>&vecList) //保存结果至文件中

{

//询问是否需要将结果保存至文件

char chChoose = ' '

cout <<"Do you want to save the result to file? (y/n):"

chChoose = _getch()

if(chChoose != 'y' &&chChoose != 'Y')

{

return true

}

//输入文件名

string strfilename

cout <<"\nInput file name:"

cin >>strfilename

cout <<"Save result to file " <<strfilename <<"......" <<endl

//打开文件

ofstream savefile(strfilename.c_str())

if(!savefile)

{

cout <<"can't open file" <<endl

return false

}

//写入N

savefile <<ulN <<endl

//写入元素

for(vector<complex<double>>::iterator i = vecList.begin()i <vecList.end()i++)

{

savefile <<*i <<endl

}

//写入完毕

cout <<"save succeed." <<endl

//关闭文件

savefile.close()

return true

}

void FFT(unsigned long &ulN, vector<complex<double>>&vecList)

{

//得到幂数

unsigned long ulPower = 0 //幂数

unsigned long ulN1 = ulN - 1

while(ulN1 >0)

{

ulPower++

ulN1 /= 2

}

//反序

bitset<sizeof(unsigned long) * 8>bsIndex //二进制容器

unsigned long ulIndex//反转后的序号

unsigned long ulK

for(unsigned long p = 0p <ulNp++)

{

ulIndex = 0

ulK = 1

bsIndex = bitset<sizeof(unsigned long) * 8>(p)

for(unsigned long j = 0j <ulPowerj++)

{

ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0

ulK *= 2

}

if(ulIndex >p)

{

complex<double>c = vecList[p]

vecList[p] = vecList[ulIndex]

vecList[ulIndex] = c

}

}

//计算旋转因子

vector<complex<double>>vecW

for(unsigned long i = 0i <ulN / 2i++)

{

vecW.push_back(complex<double>(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN)))

}

for(unsigned long m = 0m <ulN / 2m++)

{

cout<<"\nvW[" <<m <<"]=" <<vecW[m]

}

//计算FFT

unsigned long ulGroupLength = 1 //段的长度

unsigned long ulHalfLength = 0 //段长度的一半

unsigned long ulGroupCount = 0//段的数量

complex<double>cw //WH(x)

complex<double>c1 //G(x) + WH(x)

complex<double>c2 //G(x) - WH(x)

for(unsigned long b = 0b <ulPowerb++)

{

ulHalfLength = ulGroupLength

ulGroupLength *= 2

for(unsigned long j = 0j <ulNj += ulGroupLength)

{

for(unsigned long k = 0k <ulHalfLengthk++)

{

cw = vecW[k * ulN / ulGroupLength] * vecList[j + k + ulHalfLength]

c1 = vecList[j + k] + cw

c2 = vecList[j + k] - cw

vecList[j + k] = c1

vecList[j + k + ulHalfLength] = c2

}

}

}

}

void display(unsigned long &ulN, vector<complex<double>>&vecList)

{

cout <<"\n\n===========================Display The Result=========================" <<endl

for(unsigned long d = 0d <ulNd++)

{

cout <<"X(" <<d <<")\t\t\t = " <<vecList[d] <<endl

}

}

下面为STDAFX.H文件：

// stdafx.h : 标准系统包含文件的包含文件，

// 或是常用但不常更改的项目特定的包含文件

#pragma once

#include <iostream>

#include <tchar.h>

// TODO: 在此处引用程序要求的附加头文件

下面为STDAFX.CPP文件：

// stdafx.cpp : 只包括标准包含文件的源文件

// FFT.pch 将成为预编译头

// stdafx.obj 将包含预编译类型信息

#include "stdafx.h"

// TODO: 在 STDAFX.H 中

//引用任何所需的附加头文件，而不是在此文件中引用

这个是我上个星期才交的课程设计里面的源代码，绝对可以运行的。本人也觉得应该比较容易看懂吧。

#include<iostream>

#include<math.h>

using namespace std

#define PI 3.1415927

#define MAX 20000

////////////////////////////////

//以指数形式表示的复数

/////////////////////////////////

struct CompExp

{

float AbValue

float Angle

}

//////////////////////////////////////

//以一般形式表示的复数

///////////////////////////////////////

struct Complex

{

float Re

float Im

}A[MAX]

int n

float a[MAX]

//////////////////////////////////////////

//将一个数的二进制数反转之后的新数返回

//例如6->110->011->3

//////////////////////////////////////////

int Rev(int i)

{

int index,s//index为要返回的数

s=(int)log2(n)

index=0

while(s>0)

{

s--

index+=(i%2)*(int)pow(2,s)

i=i/2

}

return index

}

////////////////////////////////////////////

//将输入的离散数据进行位反置转换，并复制在复数数组A中

////////////////////////////////////////////

void bit(float* a, Complex* A)

{

for(int i=0i<ni++)

{

A[Rev(i)].Re=a[i]

A[Rev(i)].Im=0

}

}

/////////////////////////////////////////////

//复数的加法a+b=c

/////////////////////////////////////////////

void Add(Complex* a, Complex* b, Complex* c)

{

c->Re = a->Re + b->Re

c->Im = a->Im + b->Im

}

/////////////////////////////////////////////

//复数的减法a-b=c

/////////////////////////////////////////////

void Sub(Complex* a, Complex* b, Complex* c)

{

c->Re = a->Re - b->Re

c->Im = a->Im - b->Im

}

////////////////////////////////////////////

//复数的乘法a*b=c

////////////////////////////////////////////

void Mul(Complex* a, Complex* b, Complex* c)

{

float tempRe,tempIm

tempRe = a->Re * b->Re - a->Im * b->Im

tempIm = a->Re * b->Im + a->Im * b->Re

c->Re = tempRe

c->Im = tempIm

}

/////////////////////////////////////////

//将复数从指数形式转化为一般形式

/////////////////////////////////////////

void CompTrans(CompExp* a, Complex* b)

{

b->Re = a->AbValue * cos(a->Angle)

b->Im = a->AbValue * sin(a->Angle)

}

///////////////////////////////////////

//基于迭代的快速傅里叶算法

///////////////////////////////////////

void FFT()

{

Complex Wm,W,u,t

CompExp wm,w

bit(a,A)

for(int i=1i<=(int)log2(n)i++)

{

int m=(int)pow(2,i)

wm.AbValue = 1

wm.Angle = -2*PI/m

CompTrans(&wm,&Wm)

for(int k=0k<nk+=m)

{

w.AbValue = 1

w.Angle = 0

CompTrans(&w,&W)

for(int j=0j<m/2j++)

{

Mul(&W,&A[k+j+m/2],&t)

u=A[k+j]//蝴蝶 *** 作

Add(&u,&t,&A[k+j])

Sub(&u,&t,&A[k+j+m/2])

Mul(&W,&Wm,&W)

}

}

}

}

///////////////////////////////////

//测试输入的离散点个数n是否是二的幂

///////////////////////////////////

int test(int i)

{

while(1)

{

if(i==2)break

if(i/2*2!=i)return 0

i=i/2

}

return 1

}

//////////////////////////////////

//管理输入的函数

//////////////////////////////////

void Input()

{

cout<<"请输入离散数据的个数，合法的输入为2的正整数幂："<<endl

cin>>n

while(!test(n))

{

cout<<"序列长度不合法，请重新输入："<<endl

cin>>n

}

for(int i=0i<ni++)

{

cout<<"x["<<i<<"]= "

cin>>a[i]

}

}

//////////////////////////////////

//管理输出的函数

//////////////////////////////////

void Output()

{

for(int i=0i<ni++)

cout<<"X["<<i<<"]= ("<<A[i].Re<<") + j("<<A[i].Im<<")"<<endl

cout<<endl

}

///////////////////////////////////

//程序入口

///////////////////////////////////

int main()

{

while(1)

{

Input()

FFT()

Output()

}

}

clc

clear all

%输入f、N、T、是否补零（补几个零）

f=input('Input frequency of the signal: f\n')

N=input('Input number of pointsl: N\n')

T=input('Input sampling time: T\n')

flag=input('Add zero too sampling signal or not? yes=1 no=0\n')

if(flag)

ZeroNum=input('Input nmber of zeros\n')

else

ZeroNum=0

end

%生成信号，signal是原信号。signal为采样信号。

fs=1/T

t=0:0.00001:T*(N+ZeroNum-1)

signal=sin(2*pi*f*t)

t2=0:T:T*(N+ZeroNum-1)

signal2=sin(2*pi*f*t2)

if (flag)

signal2=[signal2 zeros(1, ZeroNum)]

end

%画出原信号及采样信号。

figure

subplot(2,1,1)

plot(t,signal)

xlabel('Time(s)')

ylabel('Amplitude(volt)')

title('Singnal')

hold on

subplot(2,1,1)

stem(t2,signal2,'r')

axis([0 T*(N+ZeroNum) -1 1])

%作FFT变换，计算其幅值，归一化处理，并画出频谱。

Y = fft(signal2,N)

Pyy = Y.* conj(Y)

Pyy=(Pyy/sum(Pyy))*2

f=0:fs/(N-1):fs/24

subplot(2,1,2)

bar(f,Pyy(1:N/2))

xlabel('Frequency(Hz)')

ylabel('Amplitude')

title('Frequency compnents of signal')

axis([0 fs/2 0 ceil(max(Pyy))])

grid on

---