這是我寫的1024點的快速傅里葉變換程序，下面有驗證，你把數(shù)組
讓客戶滿意是我們工作的目標(biāo)，不斷超越客戶的期望值來自于我們對這個行業(yè)的熱愛。我們立志把好的技術(shù)通過有效、簡單的方式提供給客戶，將通過不懈努力成為客戶在信息化領(lǐng)域值得信任、有價值的長期合作伙伴，公司提供的服務(wù)項目有：國際域名空間、虛擬主機、營銷軟件、網(wǎng)站建設(shè)、威海網(wǎng)站維護、網(wǎng)站推廣。
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)
的程序可以針對任何點數(shù)，只要是2的n次方
具體程序如下：
#include
iostream.h
#include
"math.h"
#includestdio.h
#includestring.h
#include
stdlib.h
#include
fstream.h
#include
afx.h
void
FFT(double
data[],
int
nn,
int
isign)
{
//復(fù)數(shù)的快速傅里葉變換
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)
//這個循環(huán)進行的是碼位倒置。
{
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;
//這里表示一次的數(shù)字的變化。也體現(xiàn)了級數(shù)，若第一級時，也就是書是的第0級，其為兩個虛數(shù)，所以對應(yīng)數(shù)組應(yīng)該增加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];//這兩句表示蝶形因子的下一個數(shù)乘以W因子所得的實部和虛部。
tempi=double(wr)*data[j+1]+double(wi)*data[j];
data[j]
=
data[i]
-
tempr;
//蝶形單元計算后下面單元的實部，下面為虛部，注意其變換之后的數(shù)組序號與書上蝶形單元是一致的
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()
{
//本程序已經(jīng)和MATLAB運算結(jié)果對比，準(zhǔn)確無誤，需要注意的的是，計算中數(shù)組都是從1開始取得，丟棄了A[0]等數(shù)據(jù)
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中的數(shù)據(jù),
其中的數(shù)據(jù)格式見該文件
FILE
*fp
=
fopen("data1024.txt","r");
if(!fp)
{
cout"Open
file
is
failing!"endl;
return;
}
while(!feof(fp))
//feof(fp)有兩個返回值:如果遇到文件結(jié)束，函數(shù)feof（fp）的值為1，否則為0。
{
memset(line,0,50);
//清空為0
memset(dataA,0,20);
memset(dataB,0,20);
fgets(line,50,fp);
//函數(shù)的功能是從fp所指文件中讀入n-1個字符放入line為起始地址的空間內(nèi)
sscanf(line,
"%s%s",
dataA,
dataB);
//我同時讀入了兩列值，但你要求1024個，那么我就只用了第一列的1024個值
//dataA讀入第一列，dataB讀入第二列
B[ij]=atof(dataA);
//將字符型的dataA值轉(zhuǎn)化為float型
ij++;
}
for
(int
mm=1;mm1025;mm++)//A[2*mm-1]是實部，A[2*mm]是虛部，當(dāng)只要輸入實數(shù)時，那么保證虛部A[mm*2]為零即可
{
A[2*mm-1]=B[mm];
A[2*mm]=0;
}
//*******************************************正式計算FFT
FFT(A,1024,1);
//********************************************寫入數(shù)據(jù)到workout.txt文件中
for
(int
k=1;k2049;k=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+只能在文件最后補充，光標(biāo)在結(jié)尾。沒有則創(chuàng)建
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計算的頻譜，是復(fù)數(shù)頻譜
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來說，是要輸出數(shù)據(jù)的地址。
//
size：要寫入的字節(jié)數(shù)；
//
count:要進行寫入size字節(jié)的數(shù)據(jù)項的個數(shù)；
//
fp:目標(biāo)文件指針。
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查看結(jié)果"endl;
}
你好，這是我的回答，希望可以幫到你。
1）結(jié)果討論
一，如果對信號進行同樣點數(shù)N的FFT變換，采樣頻率fs越高，則可以分析越高頻的信號；與此同時，采樣頻率越低，對于低頻信號的頻譜分辨率則越好。
二，假設(shè)采樣點不在正弦信號的波峰、波谷、以及0電壓處，頻譜則會產(chǎn)生泄露（leakage）。
三，對于同樣的采樣率fs，提高FFT的點數(shù)N，則可提高頻譜的分辨率。
四，如果采樣頻率fs小于2倍信號頻率2*fs（奈圭斯特定理），則頻譜分析結(jié)果會出錯。
五，對于（二）中泄露現(xiàn)象，可以通過在信號后面補零點解決。
2）程序及注解如下
%清除命令窗口及變量
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/2;4
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;
祝你好運！
我可以幫助你，你先設(shè)置我最佳答案后，我百度Hii教你。
float ar[1024],ai[1024];/* 原始數(shù)據(jù)實部，虛部 */
float a[2050];
void fft(int nn) /* nn數(shù)據(jù)長度 */
{
int n1,n2,i,j,k,l,m,s,l1;
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,};
switch(nn)
{
case 1024: s=10; break;
case 512: s=9; break;
case 256: s=8; break;
}
n1=nn/2; n2=nn-1;
j=1;
for(i=1;i=nn;i++)
{
a[2*i]=ar[i-1];
a[2*i+1]=ai[i-1];
}
for(l=1;ln2;l++)
{
if(lj)
{
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 (kj)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(i=1;i=s;i++)
{
u1=1;
u2=0;
m=(1i);
k=m1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j=k;j++)
{
for(l=j;lnn;l=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=1;i=nn/2;i++)
{
ar[i]=4*a[2*i+2]/nn; /* 實部 */
ai[i]=-4*a[2*i+3]/nn; /* 虛部 */
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); /* 幅值 */
}
}
int DFT(int dir,int m,double *x1,double *y1)
{
long i,k;
double arg;
double cosarg,sinarg;
double *x2=NULL,*y2=NULL;
x2=malloc(m*sizeof(double));
y2=malloc(m*sizeof(double));
if(x2==NULL||y2==NULL)return(FALSE);
for(i=0;im;i++)
{
x2[i]=0;
y2[i]=0;
arg=-dir*2.0*3.141592654*(double)i/(double)m;
for(k=0;km;k++)
{
cosarg=cos(k*arg);
sinarg=sin(k*arg);
x2[i]+=(x1[k]*cosarg-y1[k]*sinarg);
y2[i]+=(x1[k]*sinarg+y1[k]*cosarg);
}
}
/*Copythedataback*/
if(dir==1)
{
for(i=0;im;i++)
{
x1[i]=x2[i]/(double)m;
y1[i]=y2[i]/(double)m;
}
}
else
{
for(i=0;im;i++)
{
x1[i]=x2[i];
y1[i]=y2[i];
}
}
free(x2);
free(y2);
return(TRUE);
}
#includestdio.h
#include math.h
class complex //定義一個類,實現(xiàn)復(fù)數(shù)的所有操作
{
double Real,Image; //實部與虛部
public:
complex(double r="0",double i="0"){Real=r;Image=i;}
double GetR(){return Real;} //取出實部
double GetI(){return Image;} //取出虛部
complex operator + (complex ); //復(fù)數(shù)加法
complex operator - (complex ); //復(fù)數(shù)減法
complex operator * (complex ); //復(fù)數(shù)乘法
void operator =(complex ); //復(fù)數(shù) 賦值
};
complex complex::operator + (complex c) //復(fù)數(shù)加法
{
complex t;
t.Real=Real+c.Real;
t.Image=Image+c.Image;
return t;
}
complex complex::operator - (complex c) //復(fù)數(shù)減法
{
complex t;
t.Real=Real-c.Real;
t.Image=Image-c.Image;
return t;
}
complex complex::operator * (complex c) //復(fù)數(shù)乘法
{
complex t;
t.Real=Real*c.Real-Image*c.Image;
t.Image=Real*c.Image+Image*c.Real;
return t;
}
void complex::operator = (complex c) //復(fù)數(shù) 賦值
{
Real=c.Real;
Image=c.Image;
}
void fft(complex a[],int length,int jishu) //實現(xiàn)fft的函數(shù)
{
const double PI="3".141592653589793;
complex u,Wn,t;
int i,j,k,m,kind,distance,other;
double tmp;
for(i=0;ilength;i++) //實現(xiàn)倒敘排列
{
k="i";
j=0;
for(m=0;mjishu;m++)
{
j="j"*2+k%2;
k/=2;
}
if(ij)
{
t="a";
a=a[j];
a[j]=t;
}
}
for(m=1;m=jishu;m++) //第m級蝶形運算,總級數(shù)為jishu
{
kind = (int)pow(2,m-1); //第m級有2^(m-1)種蝶形運算
distance = 2*kind; //同種蝶形結(jié)相鄰距離為2^m
u=complex(1,0); //旋轉(zhuǎn)因子初始值為 1
tmp=PI/kind;
Wn=complex(cos(tmp),-sin(tmp));//旋轉(zhuǎn)因子Wn
for(j=0;jkind;j++) //每種蝶形運算的起始點為j,共有kind種
{
for(i=j;ilength;i+=distance) //同種蝶形運算
{
other=i+kind;//蝶形運算的兩個因子對應(yīng)單元下標(biāo)的距離為2^(m-1)
t=a[other]*u; // 蝶形運算的乘積項
a[other]=a-t; //蝶形運算
a=a+t; //蝶形運算
}
u="u"*Wn; //修改旋轉(zhuǎn)因子,多乘一個基本DFT因子WN
}
}
}
void main(void)
{
double a,b;
complex x[8]; //此程序以8點序列測試
printf("8點序列:\n");
for(int i="0";i8;i++) //初始化并輸出原始序列
{
x=complex(i,i+1);
printf("x(%d) = %lf + %lf i\n",i+1,x.GetR(),x.GetI());
}
fft(x,8,3); //調(diào)用fft函數(shù)
printf("fft變換的結(jié)果為:\n");
for(i=0;i8;i++) //輸出結(jié)果
printf("X(%d)= %lf + %lf i\n",i+1,x.GetR(),x.GetI());
}
快速傅氏變換，是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改進獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計算機系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說是進了一大步。
設(shè)x(n)為N項的復(fù)數(shù)序列，由DFT變換，任一X（m）的計算都需要N次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法，而一次復(fù)數(shù)乘法等于四次實數(shù)乘法和兩次實數(shù)加法，一次復(fù)數(shù)加法等于兩次實數(shù)加法，即使把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次“運算”（四次實數(shù)乘法和四次實數(shù)加法），那么求出N項復(fù)數(shù)序列的X（m）, 即N點DFT變換大約就需要N2次運算。當(dāng)N=1024點甚至更多的時候，需要N2=1048576次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列（設(shè)N=2k,k為正整數(shù)），分為兩個N/2項的子序列，每個N/2點DFT變換需要（N/2）2次運算，再用N次運算把兩個N/2點的DFT 變換組合成一個N點的DFT變換。這樣變換以后，總的運算次數(shù)就變成N+2（N/2）2=N+N2/2。繼續(xù)上面的例子，N=1024時，總的運算次數(shù)就變成了525312次，節(jié)省了大約50%的運算量。而如果我們將這種“一分為二”的思想不斷進行下去，直到分成兩兩一組的DFT運算單元，那么N點的 DFT變換就只需要Nlog2N次的運算，N在1024點時，運算量僅有10240次，是先前的直接算法的1%，點數(shù)越多，運算量的節(jié)約就越大，這就是 FFT的優(yōu)越性.
傅里葉變換（Transformée de Fourier）是一種積分變換。因其基本思想首先由法國學(xué)者傅里葉系統(tǒng)地提出，所以以其名字來命名以示紀(jì)念。
應(yīng)用
傅里葉變換在物理學(xué)、數(shù)論、組合數(shù)學(xué)、信號處理、概率論、統(tǒng)計學(xué)、密碼學(xué)、聲學(xué)、光學(xué)、海洋學(xué)、結(jié)構(gòu)動力學(xué)等領(lǐng)域都有著廣泛的應(yīng)用（例如在信號處理中，傅里葉變換的典型用途是將信號分解成幅值分量和頻率分量）。
概要介紹
傅里葉變換能將滿足一定條件的某個函數(shù)表示成三角函數(shù)（正弦和/或余弦函數(shù)）或者它們的積分的線性組合。在不同的研究領(lǐng)域，傅里葉變換具有多種不同的變體形式，如連續(xù)傅里葉變換和離散傅里葉變換。最初傅里葉分析是作為熱過程的解析分析的工具被提出的（參見：林家翹、西格爾著《自然科學(xué)中確定性問題的應(yīng)用數(shù)學(xué)》，科學(xué)出版社，北京。原版書名為 C. C. Lin L. A. Segel, Mathematics Applied to Deterministic Problems in the Natural Sciences, Macmillan Inc., New York, 1974）。
傅里葉變換屬于諧波分析。
傅里葉變換的逆變換容易求出,而且形式與正變換非常類似;
正弦基函數(shù)是微分運算的本征函數(shù),從而使得線性微分方程的求解可以轉(zhuǎn)化為常系數(shù)的代數(shù)方程的求解.在線性時不變的物理系統(tǒng)內(nèi),頻率是個不變的性質(zhì),從而系統(tǒng)對于復(fù)雜激勵的響應(yīng)可以通過組合其對不同頻率正弦信號的響應(yīng)來獲取;
卷積定理指出:傅里葉變換可以化復(fù)雜的卷積運算為簡單的乘積運算,從而提供了計算卷積的一種簡單手段;
離散形式的傅里葉變換可以利用數(shù)字計算機快速的算出(其算法稱為快速傅里葉變換算法(FFT)).
重慶分公司
重慶分公司
