{
int nn,n1,n2,i,j,k,l,m,s,l1
float ar[1024],ai[1024] // 实部 虚部
float a[2050]
float t1,t2,x,y
float w1,w2,u1,u2,z
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,}// 优化
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,}
nn=1024
s=10
n1=nn/2 n2=nn-1
j=1
for(i=1i<=nni++)
{
a[2*i]=ar[i-1]
a[2*i+1]=ai[i-1]
}
for(l=1l<n2l++)
{
if(l<j)
{
t1=a[2*j]
t2=a[2*j+1]
a[2*j]=a[2*l]
a[2*j+1]=a[2*l+1]
a[2*l]=t1
a[2*l+1]=t2
}
斗友 k=n1
则物 while (k<j)
{
j=j-k
k=k/2
}
j=j+k
}
for(i=1i<=si++)
{
u1=1
u2=0
m=(1<<i)
k=m>>1
w1=fcos[i-1]
w2=-fsin[i-1]
for(j=1j<=kj++)
{
for(l=jl<nnl=l+m)
{
孙销液 l1=l+k
t1=a[2*l1]*u1-a[2*l1+1]*u2
t2=a[2*l1]*u2+a[2*l1+1]*u1
a[2*l1]=a[2*l]-t1
a[2*l1+1]=a[2*l+1]-t2
a[2*l]=a[2*l]+t1
a[2*l+1]=a[2*l+1]+t2
}
z=u1*w1-u2*w2
u2=u1*w2+u2*w1
u1=z
}
}
for(i=1i<=nn/2i++)
{
ar[i]=a[2*i+2]/nn
ai[i]=-a[2*i+3]/nn
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]) // 幅值
}
}
这是我写的1024点的快速傅里叶变换程序，下面有验证，你把数组double
A[2049]={0}
double
B[1100]={0}
double
powerA[1025]={0}
改成
A[256]={0}
B[130]={0}
power[129]={0}就行了，
void
FFT(double
data[],
int
nn,
int
isign)
的程序可以针对任何点数，只要是2的n次方
具体程序如下：
#include
<iostream.h>
#include
"math.h"
#include<stdio.h>
#include<string.h>
#include
<stdlib.h>
#include
<fstream.h>
#include
<afx.h>
void
FFT(double
data[],
int
nn,
int
isign)
{
//复数的快速傅里叶变换
int
n,j,i,m,mmax,istep
double
tempr,tempi,theta,wpr,wpi,wr,wi,wtemp
n
=
2
*
nn
j
=
1
for
(i
=
1
i<=n
i=i+2)
//这个循环进行的是码位倒置。
{
if(
j
>
i)
{
tempr
=
data[j]
tempi
=
data[j
+
1]
data[j]
=
data[i]
data[j
+
1]
=
data[i
+
1]
data[i]
=
tempr
data[i
+
1]
=
tempi
}
m
=
n
/
2
while
(m
>=
2
&&
j
>
m)
{
j
=
j
-
m
m
=
m
/
2
}
j
=
j
+
m
}
mmax
=
2
while(
n
>
mmax
)
{
istep
=
2
*
mmax
//这里表示一次的数字的变化。也体现了级数，若第一级时，也就是书是的第0级，其为两个虚数，所以对应数组应该增加4，这样就可以进入下一组运算
theta
=
-6.28318530717959
/
(isign
*
mmax)
wpr
=
-2.0
*
sin(0.5
*
theta)*sin(0.5
*
theta)
wpi
=
sin(theta)
wr
=
1.0
wi
=
0.0
for(
m
=
1
m<=mmax
m=m+2)
{
for
(i
=
m
i<=n
i=i+istep)
{
j
=
i
+
mmax
tempr=double(wr)*data[j]-double(wi)*data[j+1]//这两句表示蝶形因子的下一个数乘以W因子所得的实部和虚部。
tempi=double(wr)*data[j+1]+double(wi)*data[j]
data[j]
=
data[i]
-
tempr
//蝶形单元计算后下面单元的实部，下面为虚部，注意其变换之后的数组序号与书上蝶形单元是一致的
data[j
+
1]
=
data[i
+
1]
-
tempi
data[i]
=
data[i]
+
tempr
data[i
+
1]
=
data[i
+
1]
+
tempi
}
wtemp
=
wr
wr
=
wr
*
wpr
-
wi
*
wpi
+
wr
wi
=
wi
*
wpr
+
wtemp
*
wpi
+
wi
}
mmax
=
istep
}
}
void
main()
{
//本程序已经和MATLAB运算结果对比，准确无误，需要注意的的是，计算告扒中数组都是从1开始取得，丢弃了A[0]等数据
double
A[2049]={0}
double
B[1100]={0}
double
powerA[1025]={0}
char
line[50]
char
dataA[20],
dataB[20]
int
ij
char
ch1[3]="巧带\t"
char
ch2[3]="\n"
int
strl1,strl2
CString
str1,str2
ij=1
//********************************读入文件data1024.txt中的数据,
其中的数据格式见该文件
FILE
*fp
=
fopen("data1024.txt","r")
if(!fp)
{
cout<<"Open
file
is
failing!"<<endl
return
}
while(!feof(fp))
//feof(fp)有两个返回值:如果遇到文件结束，函数feof（fp）的值为1，否则为0。
{
memset(line,0,50)
//清空为0
memset(dataA,0,20)
memset(dataB,0,20)
fgets(line,50,fp)
//函数的功能是从fp所指文件中读入n-1个字符放入line为起始地址的空间内
sscanf(line,
"%s%s",
dataA,
dataB)
//我同时读入孝友芦了两列值，但你要求1024个，那么我就只用了第一列的1024个值
//dataA读入第一列，dataB读入第二列
B[ij]=atof(dataA)
//将字符型的dataA值转化为float型
ij++
}
for
(int
mm=1mm<1025mm++)//A[2*mm-1]是实部，A[2*mm]是虚部，当只要输入实数时，那么保证虚部A[mm*2]为零即可
{
A[2*mm-1]=B[mm]
A[2*mm]=0
}
//*******************************************正式计算FFT
FFT(A,1024,1)
//********************************************写入数据到workout.txt文件中
for
(int
k=1k<2049k=k+2)
{
powerA[(k+1)/2]=sqrt(pow(A[k],2.0)+pow(A[k+1],2.0))//求功率谱
FILE
*pFile=fopen("workout.txt","a+")
//?a+只能在文件最后补充，光标在结尾。没有则创建
memset(ch1,0,15)
str1.Format("%.4f",powerA[(k+1)/2])
if
(A[k+1]>=0)
str2.Format("%d\t%6.4f%s%6.4f
%s",(k+1)/2,A[k],"+",A[k+1],"i")//保存fft计算的频谱，是复数频谱
else
str2.Format("%d\t%6.4f%6.4f
%s",(k+1)/2,A[k],A[k+1],"i")
strl1=strlen(str1)
strl2=strlen(str2)
//
用
法:fwrite(buffer,size,count,fp)
//
buffer：是一个指针，对fwrite来说，是要输出数据的地址。
//
size：要写入的字节数；
//
count:要进行写入size字节的数据项的个数；
//
fp:目标文件指针。
fwrite(str2,1,strl2,pFile)
fwrite(ch1,1,3,pFile)
fwrite(ch1,1,3,pFile)
fwrite(str1,1,strl1,pFile)
fwrite(ch2,1,3,pFile)
fclose(pFile)
}
cout<<"计算完毕，到fft_test\workout.txt查看结果"<<endl
}
分类: 教育/科学 >>学习帮助问题描述:
追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 中
引用任何所需的附加头文件，而不是在此文件中引用
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）