此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 <stdio.h>#include <math.h>
#include <stdlib.h>
#define N 1000
/*定义复数类型*/
typedef struct{
double real
double img
}complex
complex x[N], *W/*输入序列,变换核*/
int size_x=0 /*输入序列的大小，在本程序中仅限2的次幂*/
double PI/*圆周率*/
void fft()/*快速傅里叶变换*/
void initW() /*初始化变换核*/
void change()/*变址*/
void add(complex ,complex ,complex *)/*复数加法*/
void mul(complex ,complex ,complex *)/*复数乘法*/
void sub(complex ,complex ,complex *)/*复数减法*/
void output()
int main(){
int i/*输出结果*/
system("cls")
PI=atan(1)*4
printf("Please input the size of x:\n")
scanf("%d",&size_x)
printf("Please input the data in x[N]:\n")
for(i=0i<size_xi++)
scanf("%lf%lf",&x[i].real,&x[i].img)
initW()
fft()
output()
return 0
}
/*快速傅里叶变换*/
void fft(){
int i=0,j=0,k=0,l=0
complex up,down,product
change()
for(i=0i<log(size_x)/log(2) i++){ /*一级蝶形运算*/
l=1<<i
for(j=0j<size_xj+= 2*l ){ /*一组蝶形运算*/
for(k=0k<lk++){/*一个蝶形运算*/
mul(x[j+k+l],W[size_x*k/2/l],&product)
add(x[j+k],product,&up)
sub(x[j+k],product,&down)
x[j+k]=up
x[j+k+l]=down
}
}
}
}
/*初始化变换核*/
void initW(){
int i
W=(complex *)malloc(sizeof(complex) * size_x)
for(i=0i<size_xi++){
W[i].real=cos(2*PI/size_x*i)
W[i].img=-1*sin(2*PI/size_x*i)
}
}
/*变址计算，将x(n)码位倒置*/
void change(){
complex temp
unsigned short i=0,j=0,k=0
double t
for(i=0i<size_xi++){
k=ij=0
t=(log(size_x)/log(2))
while( (t--)>0 ){
j=j<<1
j|=(k &1)
k=k>>1
}
if(j>i){
temp=x[i]
x[i]=x[j]
x[j]=temp
}
}
}
/*输出傅里叶变换的结果*/
void output(){
int i
printf("The result are as follows\n")
for(i=0i<size_xi++){
printf("%.4f",x[i].real)
if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img)
else if(fabs(x[i].img)<0.0001)printf("\n")
else printf("%.4fj\n",x[i].img)
}
}
void add(complex a,complex b,complex *c){
c->real=a.real+b.real
c->img=a.img+b.img
}
void mul(complex a,complex b,complex *c){
c->real=a.real*b.real - a.img*b.img
c->img=a.real*b.img + a.img*b.real
}
void sub(complex a,complex b,complex *c){
c->real=a.real-b.real
c->img=a.img-b.img
}
function y=myditfft(x)%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次
%------------------------------------
%
myditfft.c
%------------------------------------
m=nextpow2(x)
%求的x长度对应的2的最低幂次m
N=2^m
if length(x)<N
x=[x,zeros(1,N-length(x))]
%若的长度不是2的幂，补0到2的整数幂
end
nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1
%求1：2^m数列的倒序
y=x(nxd)
%将倒序排列作为的初始值
for mm=1:m
%将DFT做m次基2分解，从左到右，对每次分解作DFT运算
Nmr=2^mm
u=1
%旋转因子u初始化
WN=exp(-i*2*pi/Nmr)
%本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
for j=1:Nmr/2
%本次跨越间隔内的各次碟形运算
for k=j:Nmr:N
%本次碟形运算的跨越间隔为Nmr=2^mm
kp=k+Nmr/2
%确定碟形运算的对应单元下标
t=y(kp)*u
%碟形运算的乘积项
y(kp)=y(k)-t
%碟形运算的加法项
y(k)=y(k)+t
end
u=u*WN
%修改旋转因子，多乘一个基本DFT因子WN
end
end
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）