问题描述:
追20分
解析:
快速傅里叶变换 要用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 <plex>
#include <bitset>
#include <vector>
#include <conio.h>
#include <string>
#include <fstream>
using namespace std
bool inputData(unsigned long &, vector<plex<double>>&)手工输入数据
void FFT(unsigned long &, vector<plex<double>>&)FFT变换
void display(unsigned long &, vector<plex<double>>&)显示结果
bool readDataFromFile(unsigned long &, vector<plex<double>>&)从文件中读取数据
bool saveResultToFile(unsigned long &, vector<plex<double>>&)保存结果至文件中
const double PI = 3.1415926
int _tmain(int argc, _TCHAR* argv[])
{
vector<plex<double>>vecList有限长序列
unsigned long ulN = 0N
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<plex<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()清空原有序列
plex<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<plex<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
}
读取元素
plex<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<plex<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<plex<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<plex<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)
{
plex<double>c = vecList[p]
vecList[p] = vecList[ulIndex]
vecList[ulIndex] = c
}
}
计算旋转因子
vector<plex<double>>vecW
for(unsigned long i = 0i <ulN / 2i++)
{
vecW.push_back(plex<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段的数量
plex<double>cwWH(x)
plex<double>c1G(x) + WH(x)
plex<double>c2G(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<plex<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 中
引用任何所需的附加头文件，而不是在此文件中引用
设x（N）为N点有限长离散序列，代入式（8－3）、式（8－4），并令 其傅里叶变换（DFT）为
地球物理数据处理基础
反变换（IDFT）为
地球物理数据处理基础
两者的差异只在于W的指数符号不同，以及差一个常数1/N，因此下面我们只讨论DFT正变换式（8－5）的运算量，其反变换式（8－6）的运算是完全相同的。
一般来说，W是复数，因此，X（j）也是复数，对于式（8－5）的傅里叶变换（DFT），计算一个X（j）值需要N次复数乘法和N－1次复数加法。而X（j）一共有N个值（j＝0，1，…，N－1），所以完成整个DFT运算总共需要N2次复数乘法和N（N－1）次复数加法。
直接计算DFT，乘法次数和加法次数都是与N2成正比的，当N很大时，运算量会很大，例如，当N＝8时，DFT需64次复数乘法；而当N＝1024时，DFT所需乘法为1048576次，即一百多万次的复数乘法运算，对运算速度要求高。所以需要改进DFT的计算方法，以减少运算次数。
分析Wjk，表面上有N2个数值，由于其周期性，实际上仅有N个不同的值W0，W1，…，WN－1。对于N＝2m时，由于其对称性，只有N/2个不同的值W0，W1，…，
地球物理数据处理基础
因此可以把长序列的DFT分解为短序列DFT，而前面已经分析DFT与N2成正比，所以N越小越有利。同时，利用ab＋ac＝a（b＋c）结合律法则，可以将同一个Wr对应的系数x（k）相加后再乘以Wr，就能大大减少运算次数。这就是快速傅里叶变换（FFT）的算法思路。下面，我们来分析N＝2m情况下的FFT算法。
1.N＝4的FFT算法
对于m＝2，N＝4，式（8－5）傅里叶变换为
地球物理数据处理基础
将式（8－7）写成矩阵形式
地球物理数据处理基础
为了便于分析，将上式中的j，k写成二进制形式，即
地球物理数据处理基础
代入式（8－7），得
地球物理数据处理基础
分析Wjk的周期性来减少乘法次数
地球物理数据处理基础
则 代回式（8－9），整理得
地球物理数据处理基础
上式可分层计算，先计算内层，再计算外层时就利用内层计算的结果，可避免重复计算。写成分层形式
地球物理数据处理基础
则X（j1 j0）＝X2（j1 j0）。
上式表明对于N＝4的FFT，利用Wr的周期关系可分为m＝2步计算。实际上，利用Wr的对称性，仍可以对式（8－11）进行简化计算。考虑到
地球物理数据处理基础
式（8－11）可以简化为
地球物理数据处理基础
令j＝j0；k＝k0，并把上式表示为十进制，得
地球物理数据处理基础
可以看到，完成上式N＝4的FFT计算（表8－1）需要N·（m－1）/2＝2次复数乘法和N·m＝8次复数加法，比N＝4的DFT算法的N2＝16次复数乘法和N·（N－1）＝12次复数加法要少得多。
表8－1 N＝4的FFT算法计算过程
注：W0＝1；W1＝－i。
［例1］求N＝4样本序列1，3，3，1的频谱（表8－2）。
表8－2 N＝4样本序列
2.N＝8的FFT算法
类似N＝4的情况，用二进制形式表示，有
地球物理数据处理基础
写成分层计算的形式：
地球物理数据处理基础
则X（j2 j1 j0）＝X3（j2 j1 j0）。
对式（8－14）的X1（k1 k0 j0）进行展开，有
地球物理数据处理基础
还原成十进制，并令k＝2k1＋k0，即k＝0，1，2，3，有
地球物理数据处理基础
用类似的方法对式（8－14）的X2（k0 j1 j0），X3（j2 j1 j0）进行展开，整理得
地球物理数据处理基础
用式（8－16）、式（8－17）逐次计算到X3（j）＝X（j）（j＝0，1，…，7），即完成N＝23＝8的FFT计算，其详细过程见表8－3。
表8－3 N＝8的FFT算法计算过程
注：对于正变换 对于反变换 所
［例2］求N＝8样本序列（表8－4）x（k）＝1，2，1，1，3，2，1，2的频谱。
表8－4 N＝8样本序列
3.任意N＝2m的FFT算法
列出N＝4，N＝8的FFT计算公式，进行对比
地球物理数据处理基础
观察式（8－18）、式（8－19），不难看出，遵循如下规律：
（1）等式左边的下标由1递增到m，可用q＝1，2，…，m代替，则等式右边为q－1；
（2）k的上限为奇数且随q的增大而减小，至q＝m时为0，所以其取值范围为k＝0，1，2，…，（2m－q－1）；
（3）j的上限为奇数且随q的增大而增大，且q＝1时为0，其取值范围为j＝0，1，2，…，（2q－1－1）；
（4）k的系数，在等式左边为2q，等式右边为2q－1（包括W的幂指数）；
（5）等式左边序号中的常数是2的乘方形式，且幂指数比下标q小1，即2q－1；等式右边m对式子序号中的常数都是定值2m－1。
归纳上述规则，写出对于任意正整数m，N＝2m的FFT算法如下：
由X0（p）＝x（p）（p＝0，1，…，N－1）开始：
（1）对q＝1，2，…，m，执行（2）～（3）步；
（2）对k＝0，1，2，…，（2m－q－1）及j＝0，1，2，…，（2q－1－1），执行
地球物理数据处理基础
（3）j，k循环结束；
（4）q循环结束；由Xm（p）（p＝0，1，…，N－1）输出原始序列x（p）的频谱X（p）。
在计算机上很容易实现上述FFT算法程序，仅需要三个复数数组，编程步骤如下：
（1）设置复数数组X1（N－1），X2（N－1）和 （数组下界都从0开始）；
（2）把样本序列x赋给X1，即X1（k）＝x（k）（k＝0，1，…，N－1）；
（3）计算W，即正变换 反变换
（4）q＝1，2，…，m，若q为偶数，执行（6），否则执行第（5）步；
（5）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作
X2（2qk＋j）＝X1（2q－1k＋j）＋X1（2q－1k＋j＋2m－1）
X2（2qk＋j＋2q－1）＝［X1（2q－1k＋j）－X1（2q－1k＋j＋2m－1）］W（2q－1k）
至k，j循环结束；
（6）k＝0，1，2，…，（2m－q－1）和j＝0，1，2，…，（2q－1－1）循环，作
X1（2qk＋j）＝X2（2q－1k＋j）＋X2（2q－1k＋j＋2m－1）
X1（2qk＋j＋2q－1）＝［X2（2q－1k＋j）－X2（2q－1k＋j＋2m－1）］W（2q－1k）
至k，j循环结束；
（7）q循环结束，若m为偶数，输出X1（j），否则输出X2（j）（j＝0，1，…，N－1），即为所求。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）