新聞中心
一個關于128點的快速傅立葉的C語言程序
這是我寫的1024點的快速傅里葉變換程序,下面有驗證,你把數(shù)組
公司主營業(yè)務:成都網(wǎng)站制作、做網(wǎng)站、移動網(wǎng)站開發(fā)等業(yè)務。幫助企業(yè)客戶真正實現(xiàn)互聯(lián)網(wǎng)宣傳,提高企業(yè)的競爭能力。成都創(chuàng)新互聯(lián)公司是一支青春激揚、勤奮敬業(yè)、活力青春激揚、勤奮敬業(yè)、活力澎湃、和諧高效的團隊。公司秉承以“開放、自由、嚴謹、自律”為核心的企業(yè)文化,感謝他們對我們的高要求,感謝他們從不同領域給我們帶來的挑戰(zhàn),讓我們激情的團隊有機會用頭腦與智慧不斷的給客戶帶來驚喜。成都創(chuàng)新互聯(lián)公司推出南潯免費做網(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)
{
//復數(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ù),所以對應數(shù)組應該增加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運算結果對比,準確無誤,需要注意的的是,計算中數(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)有兩個返回值:如果遇到文件結束,函數(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值轉化為float型
ij++;
}
for
(int
mm=1;mm1025;mm++)//A[2*mm-1]是實部,A[2*mm]是虛部,當只要輸入實數(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+只能在文件最后補充,光標在結尾。沒有則創(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計算的頻譜,是復數(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:目標文件指針。
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;
}
傅里葉變換常用公式有哪些?
1、門函數(shù)F(w)=2w w sin=Sa() w。
2、指數(shù)函數(shù)(單邊)f(t)=e-atu(t) F(w)=1,實際上是一個低通濾波器a+jw。
3、單位沖激函數(shù)F(w)=1,頻帶無限寬,是一個均勻譜。
4、常數(shù)1 常數(shù)1是一個直流信號,所以它的頻譜當然只有在w=0的時候才有值,體現(xiàn)為(w)。F(w)=2(w) 可以由傅里葉變換的對稱性得到。
5、正弦函數(shù)F(ejw0t)=2(w-w0),相當于是直流信號的移位。F(sinw0t)=F((ejw0t-e-jw0t)/2)=((w-w0)-(w+w0))F(sinw0t)=F((e。
6、單位沖擊序列jw0t-e-jw0t)/2j)=j((w-w0)-(w+w0)) T(t)=(t-Tn) -這是一個周期函數(shù),每隔T出現(xiàn)一個沖擊,周期函數(shù)的傅里葉變換是離散的F(T(t))=w0(w-nw0)=w0,w0(w) n=-單位沖擊序列的傅里葉變換仍然是周期序列,周期是w0=2T。
傅立葉變換:
傅立葉變換是指將滿足一定條件的某個函數(shù)表示成三角函數(shù)的積分。傅立葉變換是在對傅立葉級數(shù)的研究中產(chǎn)生的。在不同的研究領域,傅立葉變換具有不同的作用。
在分析信號的時候 主要考慮的頻率、幅值、相位。
傅里葉變換的作用主要是將函數(shù)轉化成多個正弦組合(或e指數(shù))的形式,本質(zhì)上變換之后信號還是原來的信號只是換了一種表達方式 這樣可以更直觀的分析一個函數(shù)里的頻率、幅值、相位成分。
所以分析一個復雜的信號只需經(jīng)過傅里葉變換后可以輕易的看出其頻率和相位、幅度分量。
傅里葉變換的性質(zhì)及常用函數(shù)
基本性質(zhì)
線性性質(zhì)線性linear,指量與量之間按比例、成直線的關系,在數(shù)學上可以理解為一階導數(shù)為常數(shù)的函數(shù);非線性non-linear則指不按比例、不成直線的關系,一階導數(shù)不為常數(shù)。 如問:兩個眼睛的視敏度是一個眼睛的幾倍?很容易想到的是兩倍,可實際是 6-10倍!這就是非線性。激光也是非線性的!天體運動存在混沌;電、光與聲波的振蕩,會突陷混沌;地磁場在400萬年間,方向突變16次,也是由于混沌。甚至人類自己,原來都是非線性的:與傳統(tǒng)的想法相反,健康人的腦電圖和心臟跳動并不是規(guī)則的,而是混沌的,混沌正是生命力的表現(xiàn),混沌系統(tǒng)對外界的刺激反應,比非混沌系統(tǒng)快?!珊瘮?shù)之和的傅里葉變換等于各自變換之和。數(shù)學描述是:若函數(shù)f \left( x\right )和g \left(x \right)的傅里葉變換\mathcal[f]和\mathcal[g]都存在,α 和 β 為任意常系數(shù),則\mathcal[\alpha f+\beta g]=\alpha\mathcal[f]+\beta\mathcal[g];傅里葉變換算符\mathcal可經(jīng)歸一化成為么正算符;
頻移性質(zhì)
若函數(shù)f \left( x\right )存在傅里葉變換,則對任意實數(shù) ω0,函數(shù)f(x) e^{i \omega_ x}也存在傅里葉變換,且有\(zhòng)mathcal[f(x)e^{i \omega_ x}]=F(\omega + \omega _0 ) 。式中花體\mathcal是傅里葉變換的作用算子,平體F表示變換的結果(復函數(shù)),e 為自然對數(shù)的底,i 為虛數(shù)單位\sqrt;
微分關系
若函數(shù)f \left( x\right )當|x|\rightarrow\infty時的極限為0,而其導函數(shù)f'(x)的傅里葉變換存在,則有\(zhòng)mathcal[f'(x)]=-i \omega \mathcal[f(x)] ,即導函數(shù)的傅里葉變換等于原函數(shù)的傅里葉變換乘以因子 ? iω 。更一般地,若f(\pm\infty)=f'(\pm\infty)=\ldots=f^{(k-1)}(\pm\infty)=0,且\mathcal[f^{(k)}(x)]存在,則\mathcal[f^{(k)}(x)]=(-i \omega)^ \mathcal[f] ,即 k 階導數(shù)的傅里葉變換等于原函數(shù)的傅里葉變換乘以因子( ? iω)k。
卷積特性
若函數(shù)f \left( x\right )及g \left( x\right )都在(-\infty,+\infty)上絕對可積,則卷積函數(shù)f*g=\int_{-\infty}^{+\infty} f(x-\xi)g(\xi)d\xi的傅里葉變換存在,且\mathcal[f*g]=\mathcal[f]\cdot\mathcal[g] 。卷積性質(zhì)的逆形式為\mathcal^[F(\omega)G(\omega)]=\mathcal^[F(\omega)]*\mathcal^[G(\omega)] ,即兩個函數(shù)乘積的傅里葉逆變換等于它們各自的傅里葉逆變換的卷積,同時還有兩個函數(shù)卷積的傅里葉逆變換等于它們各自的傅里葉逆變換的乘積。
Parseval定理
若函數(shù)f \left( x\right )可積且平方可積,則\int_{-\infty}^{+\infty} f^2 (x)dx = \frac{2\pi}\int_{-\infty}^{+\infty} |F(\omega)|^d\omega 。其中 F(ω) 是 f(x) 的傅里葉變換。
常見函數(shù)傅里葉變換基本公式是什么?
傅里葉變換是:F(ω)=∫(∞,-∞) f(t)e^(-iωt)dt f(t) = (1/2π) ∫(∞,-∞) F(ω)e^(iωt)dω 令:f(t)=δ(t),那么:∫(∞,-∞) δ(t)e^(-iωt)dt = 1 而上式的反變換。
傅立葉變換的主要作用就是讓函數(shù)在時域和頻域可以相互轉化。最顯而易見的應用就是:當輸入函數(shù)和單位沖激響應函數(shù)都被轉化為頻域函數(shù)后,兩個頻域函數(shù)直接做乘法,就可以得到輸出的頻域函數(shù)。最后再反變換回時域,就可以得到輸出的時域函數(shù)。
簡介
因FFT是為時序電路而設計的,因此,控制信號要包括時序的控制信號及存儲器的讀寫地址,并產(chǎn)生各種輔助的指示信號。
同時在計算模塊的內(nèi)部,為保證高速,所有的乘法器都須始終保持較高的利用率。這意味著在每一個時鐘來臨時都要向這些單元輸入新的操作數(shù),而這一切都需要控制信號的緊密配合。
C語言 1024點快速傅里葉變換(FFT)程序,最好經(jīng)過優(yōu)化,執(zhí)行速度快
void fft()
{
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,};// 優(yōu)化
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=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]=a[2*i+2]/nn;
ai[i]=-a[2*i+3]/nn;
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); // 幅值
}
}
傅里葉變換用C語言程序怎么實現(xiàn)?
#include math.h
#include stdio.h
#define N 8
void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il);
void main()
{
double xr[N],xi[N],Yr[N],Yi[N],l=0,il=0;
int i,j,n=N,k=3;
for(i=0;iN;i++)
{
xr[i]=i;
xi[i]=0;
}
printf("------FFT------\n");
l=0;
kkfft(xr,xi,n,k,Yr,Yi,l,il);
for(i=0;iN;i++)
{
printf("%-11lf + j* %-11lf\n",Yr[i],Yi[i]);
}
printf("-----DFFT-------\n");
l=1;
kkfft(Yr,Yi,n,k,xr,xi,l,il);
for(i=0;iN;i++)
{
printf("%-11lf + j* %-11lf\n",xr[i],xi[i]);
}
getch();
}
void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
for (it=0; it=n-1; it++)
{
m = it;
is = 0;
for(i=0; i=k-1; i++)
{
j = m/2;
is = 2*is+(m-2*j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
pr[0] = 1.0;
pi[0] = 0.0;
p = 6.283185306/(1.0*n);
pr[1] = cos(p);
pi[1] = -sin(p);
if (l!=0)
pi[1]=-pi[1];
for (i=2; i=n-1; i++)
{
p = pr[i-1]*pr[1];
q = pi[i-1]*pi[1];
s = (pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
pr[i] = p-q;
pi[i] = s-p-q;
}
for (it=0; it=n-2; it=it+2)
{
vr = fr[it];
vi = fi[it];
fr[it] = vr+fr[it+1];
fi[it] = vi+fi[it+1];
fr[it+1] = vr-fr[it+1];
fi[it+1] = vi-fi[it+1];
}
m = n/2;
nv = 2;
for (l0=k-2; l0=0; l0--)
{
m = m/2;
nv = 2*nv;
for(it=0; it=(m-1)*nv; it=it+nv)
for (j=0; j=(nv/2)-1; j++)
{
p = pr[m*j]*fr[it+j+nv/2];
q = pi[m*j]*fi[it+j+nv/2];
s = pr[m*j]+pi[m*j];
s = s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr = p-q;
poddi = s-p-q;
fr[it+j+nv/2] = fr[it+j]-poddr;
fi[it+j+nv/2] = fi[it+j]-poddi;
fr[it+j] = fr[it+j]+poddr;
fi[it+j] = fi[it+j]+poddi;
}
}
/*逆傅立葉變換*/
if(l!=0)
{
for(i=0; i=n-1; i++)
{
fr[i] = fr[i]/(1.0*n);
fi[i] = fi[i]/(1.0*n);
}
}
/*是否計算模和相角*/
if(il!=0)
{
for(i=0; i=n-1; i++)
{
pr[i] = sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
if(fabs(fr[i])0.000001*fabs(fi[i]))
{
if ((fi[i]*fr[i])0)
pi[i] = 90.0;
else
pi[i] = -90.0;
}
else
pi[i] = atan(fi[i]/fr[i])*360.0/6.283185306;
}
}
return;
}
分享名稱:c語言傅立葉變換函數(shù) c語言實現(xiàn)傅立葉變換
當前網(wǎng)址:http://www.dlmjj.cn/article/doocchh.html