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
}
傅立叶分析在工程中经常用到，通常可以研究变量的频率成分。
使用Excel生成待研究数据，公式为：y==2*SIN(x*2)+COS(x)+RAND()*0.1
其中包括两种不同的频率成分，且强度为两倍，还有一个强度为0.1的噪声信号分量。
使用傅立叶分析工具
得到分析的模数
使用IMABS函数计算功率密度
分析结果：
比较原信号和频谱图，可以看到在0.3和0.6处分别有两个峰值，且强度为60和120，满足输入函数。
傅立叶变换分为四种类别：
1、非周期性连续信号傅立叶变换(Fourier Transform, FT)。
2、周期性连续信号傅立叶级数(Fourier Series, FS)。
3、非周期性离散信号离散时域傅立叶变换(Discrete Time Fourier Transform, DTFT)。
4、周期性离散信号离散傅立叶变换(Discrete Fourier Transform, DFT)。
Excel的傅立叶分析是快速傅里叶变换(Fast Fourier Transform)及其逆变换。快速傅里叶变换是利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称。
快速傅里叶变换有广泛的应用：数字信号处理、计算大整数乘法、求解偏微分方程、用于判断时间序列周期性。
扩展资料：傅里叶变换的应用：
傅里叶变换是线性算子，若赋予适当的范数，它还是酉算子。傅里叶变换的逆变换容易求出，而且形式与正变换非常类似。
正弦基函数是微分运算的本征函数，从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变的物理系统内，频率是个不变的性质，从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取。
著名的卷积定理指出：傅里叶变换可以化复杂的卷积运算为简单的乘积运算，从而提供了计算卷积的一种简单手段。
离散形式的傅里叶变换可以利用数字计算机快速的算出（其算法称为快速傅里叶变换算法（FFT))。
f(t)是t的周期函数，如果t满足狄里赫莱条件：在一个以2T为周期内f(X)连续或只有有限个第一类间断点，附f（x）单调或可划分成有限个单调区间。
则F（x）以2T为周期的傅里叶级数收敛，和函数S（x）也是以2T为周期的周期函数，且在这些间断点上，函数是有限值；在一个周期内具有有限个极值点；绝对可积。
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）