快速傅里葉變換c語言
『壹』 如何用c++快速實現傅里葉變換
/Data為輸入數據指針，Log2N=log2(length),flag=-1時為正變換，flag=1時為反變換，變換結果同樣由指針Data指向的空間返回
void dft(complex<double>*Data,int Log2N,int flag)
...{
int i,j,length;
complex<double> wn;
length=1<<Log2N;
complex<double>*temp=new complex<double>(length);
for(i=0;i<length;i++)
...{
temp[i]=0;
for(j=0;j<length;j++)
...{
wn=complex<double>(cos(2.0*pi/length*i*j),sin(flag*2.0*pi/length*i*j));
temp[i]+=Data[j]*wn;
}
}
if(flag==1)
for(i=0;i<length;i++)
Data[i]=temp[i]/length;
delete[] temp;
}
『貳』 求基2、基4、基8FFT（快速傅里葉變換）的c語言程序，要能運行得出來的
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct{
double r;
double i;
}my_complex
;
//檢查a是否為2的整數次方數
#define NOT2POW(a) (((a)-1)&(a)||(a)<=0)
//pi
#define MYPI 3.14159265358979323846
my_complex* fft(const my_complex* x, unsigned int len){
unsigned int ex=0,t=len;
unsigned int i,j,k;
my_complex *y;
double tr,ti,rr,ri,yr,yi;
if(NOT2POW(len)) return NULL; //如果失敗，返回空指針
for(;!(t&1);t>>=1) ex++; //len應該等於2的ex次方
y=(my_complex*)malloc(len*sizeof(my_complex));
if(!y) return NULL;
//變址計算，庫里-圖基演算法
for(i=0;i<len;i++){
k=i;
j=0;
t=ex;
while((t--)>0){
j<<=1;
j|=k&1;
k>>=1;
}
if(j>=i){
y[i]=x[j];
y[j]=x[i];
}
}
//用變址後的y向量進行計算
for(i=0;i<ex;i++){
t=1<<i;
for(j=0;j<len;j+=t<<1){
for(k=0;k<t;k++){
ti=-MYPI*k/t;
rr=cos(ti);
ri=sin(ti);
tr=y[j+k+t].r;
ti=y[j+k+t].i;
yr=rr*tr-ri*ti;
yi=rr*ti+ri*tr;
tr=y[j+k].r;
ti=y[j+k].i;
y[j+k].r=tr+yr;
y[j+k].i=ti+yi;
y[j+k+t].r=tr-yr;
y[j+k+t].i=ti-yi;
}
}
}
return y;
}
//以下為測試
int main()
{
int i,DATA_LEN;
my_complex *x,*y;
printf("基二FFT測試\n輸入生成序列長度:");
scanf("%d",&DATA_LEN);
x=(my_complex*)malloc(DATA_LEN*sizeof(my_complex));
for(i=0;i<DATA_LEN;i++){
x[i].r=i;
x[i].i=i-1;
}
printf("處理前...\n實部\t\t虛部\n");
for(i=0;i<DATA_LEN;i++)
printf("%lf\t%lf\n",x[i].r,x[i].i);
y=fft(x,DATA_LEN);
if(!y){
printf("序列長度不為2的整數次方!\n");
return 0;
}
printf("處理後...\n實部\t\t虛部\n");
for(i=0;i<DATA_LEN;i++)
printf("%lf\t%lf\n",y[i].r,y[i].i);
free(y);
free(x);
return 0;
}
『叄』 求個快速傅里葉變換的C語言程序
voidfft()
{
intnn,n1,n2,i,j,k,l,m,s,l1;
floatar[1024],ai[1024];//實部虛部
floata[2050];
floatt1,t2,x,y;
floatw1,w2,u1,u2,z;
floatfsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,};//優化
floatfcos[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;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]=a[2*i+2]/nn;
ai[i]=-a[2*i+3]/nn;
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]);//幅值
}
}
『肆』 一個關於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;
}
『伍』 求FFT的c語言程序
快速傅里葉變換 要用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 <complex>
#include <bitset>
#include <vector>
#include <conio.h>
#include <string>
#include <fstream>
using namespace std;
bool inputData(unsigned long &, vector<complex<double> >&); //手工輸入數據
void FFT(unsigned long &, vector<complex<double> >&); //FFT變換
void display(unsigned long &, vector<complex<double> >&); //顯示結果
bool readDataFromFile(unsigned long &, vector<complex<double> >&); //從文件中讀取數據
bool saveResultToFile(unsigned long &, vector<complex<double> >&); //保存結果至文件中
const double PI = 3.1415926;
int _tmain(int argc, _TCHAR* argv[])
{
vector<complex<double> > vecList; //有限長序列
unsigned long ulN = 0; //N
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<complex<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(); //清空原有序列
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
cout << "Input x(" << i << "):";
cin >> c;
vecList.push_back(c);
}
return true;
}
bool readDataFromFile(unsigned long & ulN, vector<complex<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;
}
//讀取元素
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
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<complex<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<complex<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<complex<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 = 0; p < ulN; p++)
{
ulIndex = 0;
ulK = 1;
bsIndex = bitset<sizeof(unsigned long) * 8>(p);
for(unsigned long j = 0; j < ulPower; j++)
{
ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0;
ulK *= 2;
}
if(ulIndex > p)
{
complex<double> c = vecList[p];
vecList[p] = vecList[ulIndex];
vecList[ulIndex] = c;
}
}
//計算旋轉因子
vector<complex<double> > vecW;
for(unsigned long i = 0; i < ulN / 2; i++)
{
vecW.push_back(complex<double>(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN)));
}
for(unsigned long m = 0; m < ulN / 2; m++)
{
cout<< "\nvW[" << m << "]=" << vecW[m];
}
//計算FFT
unsigned long ulGroupLength = 1; //段的長度
unsigned long ulHalfLength = 0; //段長度的一半
unsigned long ulGroupCount = 0; //段的數量
complex<double> cw; //WH(x)
complex<double> c1; //G(x) + WH(x)
complex<double> c2; //G(x) - WH(x)
for(unsigned long b = 0; b < ulPower; b++)
{
ulHalfLength = ulGroupLength;
ulGroupLength *= 2;
for(unsigned long j = 0; j < ulN; j += ulGroupLength)
{
for(unsigned long k = 0; k < ulHalfLength; k++)
{
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<complex<double> >& vecList)
{
cout << "\n\n===========================Display The Result=========================" << endl;
for(unsigned long d = 0; d < ulN;d++)
{
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 中
//引用任何所需的附加頭文件，而不是在此文件中引用
『陸』 c語言實現音樂信號的快速傅里葉變換，為什麼要有周期中斷來ad採集音樂，這個周期採集的周期怎麼確定
頻率和周期互為倒數。 f = 1/T; T=1/f;
f = 40khz = 40000 hz = 40000 ( 1 秒 多少次 叫 多少 赫茲）；
T = 1/f = 1.0 / 40000.0; 采樣的時間間隔。
傅里葉變換 -- 時域到頻域變換，用於研究時序信號的頻域特性
快速傅里葉變換 -- 數據點數 必須是 2 的整數次方，例如 1024，2048，4096 。。。。不足時要補點，補點有多種方法，最常用是補0 或 假定信號從頭再來。
-------
你想用 40000 hz 采樣頻率，1秒就要有 40000 點。
總的信號長度若是幾分鍾，點數就嚇人地多，FFT 耗時也要很長。
你也許可以把信號分段來FFT, 或1次就分析4096點，或 設一個 t0 時步推進，分段信號 重疊 幾分之1，得動態頻譜。。。。看你的需要。
『柒』 傅里葉變換用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;
}