❶ 怎樣用C語言實現FFT演算法啊
1、二維FFT相當於對行和列分別進行一維FFT運算。具體的實現辦法如下：
先對各行逐一進行一維FFT，然後再對變換後的新矩陣的各列逐一進行一維FFT。相應的偽代碼如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N);
for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩陣的第i行。注意這只是一個簡單的記法，並不能完全照抄。還需要通過一些語句來生成各行的數據。同理，COL[i]是對矩陣的第i列的一種簡單表示方法。
所以，關鍵是一維FFT演算法的實現。
2、常式：
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#defineN1000
/*定義復數類型*/
typedefstruct{
doublereal;
doubleimg;
}complex;
complexx[N],*W;/*輸入序列,變換核*/
intsize_x=0;/*輸入序列的大小，在本程序中僅限2的次冪*/
doublePI;/*圓周率*/
voidfft();/*快速傅里葉變換*/
voidinitW();/*初始化變換核*/
voidchange();/*變址*/
voidadd(complex,complex,complex*);/*復數加法*/
voidmul(complex,complex,complex*);/*復數乘法*/
voidsub(complex,complex,complex*);/*復數減法*/
voidoutput();
intmain(){
inti;/*輸出結果*/
system("cls");
PI=atan(1)*4;
printf("Pleaseinputthesizeofx: ");
scanf("%d",&size_x);
printf("Pleaseinputthedatainx[N]: ");
for(i=0;i<size_x;i++)
scanf("%lf%lf",&x[i].real,&x[i].img);
initW();
fft();
output();
return0;
}
/*快速傅里葉變換*/
voidfft(){
inti=0,j=0,k=0,l=0;
complexup,down,proct;
change();
for(i=0;i<log(size_x)/log(2);i++){/*一級蝶形運算*/
l=1<<i;
for(j=0;j<size_x;j+=2*l){/*一組蝶形運算*/
for(k=0;k<l;k++){/*一個蝶形運算*/
mul(x[j+k+l],W[size_x*k/2/l],&proct);
add(x[j+k],proct,&up);
sub(x[j+k],proct,&down);
x[j+k]=up;
x[j+k+l]=down;
}
}
}
}
/*初始化變換核*/
voidinitW(){
inti;
W=(complex*)malloc(sizeof(complex)*size_x);
for(i=0;i<size_x;i++){
W[i].real=cos(2*PI/size_x*i);
W[i].img=-1*sin(2*PI/size_x*i);
}
}
/*變址計算，將x(n)碼位倒置*/
voidchange(){
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
for(i=0;i<size_x;i++){
k=i;j=0;
t=(log(size_x)/log(2));
while((t--)>0){
j=j<<1;
j|=(k&1);
k=k>>1;
}
if(j>i){
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
/*輸出傅里葉變換的結果*/
voidoutput(){
inti;
printf("Theresultareasfollows ");
for(i=0;i<size_x;i++){
printf("%.4f",x[i].real);
if(x[i].img>=0.0001)printf("+%.4fj ",x[i].img);
elseif(fabs(x[i].img)<0.0001)printf(" ");
elseprintf("%.4fj ",x[i].img);
}
}
voidadd(complexa,complexb,complex*c){
c->real=a.real+b.real;
c->img=a.img+b.img;
}
voidmul(complexa,complexb,complex*c){
c->real=a.real*b.real-a.img*b.img;
c->img=a.real*b.img+a.img*b.real;
}
voidsub(complexa,complexb,complex*c){
c->real=a.real-b.real;
c->img=a.img-b.img;
}
❷ 傅里葉變換用C語言程序怎麼實現
#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;i<N;i++)
{
xr[i]=i;
xi[i]=0;
}
printf("------FFT------\n");
l=0;
kkfft(xr,xi,n,k,Yr,Yi,l,il);
for(i=0;i<N;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;i<N;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;
}
❸ 一個關於128點的快速傅立葉的C語言程序
這是我寫的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=1;mm<1025;mm++)//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=1;k<2049;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+只能在文件最後補充，游標在結尾。沒有則創建
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;
}
❹ 求傅立葉變換的c語言程序
//希望我的回答對你的學習有幫助
#include<math.h>
#include<stdio.h>
#defineN8
voidkkfft(doublepr[],doublepi[],intn,intk,doublefr[],doublefi[],intl,intil);
voidmain()
{
doublexr[N],xi[N],Yr[N],Yi[N],l=0,il=0;
inti,j,n=N,k=3;
for(i=0;i<N;i++)
{
xr[i]=i;
xi[i]=0;
}
printf("------FFT------ ");
l=0;
kkfft(xr,xi,n,k,Yr,Yi,l,il);
for(i=0;i<N;i++)
{
printf("%-11lf+j*%-11lf ",Yr[i],Yi[i]);
}
printf("-----DFFT------- ");
l=1;
kkfft(Yr,Yi,n,k,xr,xi,l,il);
for(i=0;i<N;i++)
{
printf("%-11lf+j*%-11lf ",xr[i],xi[i]);
}
getch();
}
voidkkfft(doublepr[],doublepi[],intn,intk,doublefr[],doublefi[],intl,intil)
{
intit,m,is,i,j,nv,l0;
doublep,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語言或匯編語言實現FFT(快速傅里葉）變換，並寫出C語言或匯編代碼，萬分感謝。
float ar[1024],ai[1024];/* 原始數據實部伏帶，虛部 */
float a[2050];
void fft(int nn) /* nn數據長蘆橘度 */
{
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;l<n2;l++)
{
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=1;i<=s;i++)
{
u1=1;
u2=0;
m=(1<<i);
k=m>>1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j<缺嘩蘆=k;j++)
{
for(l=j;l<nn;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]); /* 幅值 */
}
}