傅里葉變換c語言實現
❶ 求個快速傅里葉變換的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;
}
❸ 怎樣用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;
}
❹ 求基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;
}
❺ 音頻演算法入門-傅里葉變換
上一篇文章中講了一個時域處理的演算法wsola，接下來會學習頻域處理演算法，在這之前必須得對頻域有所了解，這就不得不提傅里葉變換了，本文的目的是讓大家學會用傅里葉變換公式和傅里葉逆變換公式進行計算。數學公式是人們對世界中的現象的描述，我們學習數學公式也不該只停留在使用公式來解決問題的層次，得明白公式到底在描述什麼現象，從這些天才數學家的角度來看世界。懂的地方可跳過。項目地址在文章末尾給出。
我直接說結論，傅里葉級數公式包含了傅里葉變換和傅里葉逆變換（不嚴謹的說就是這么回事）。
先簡單說下具體關系，法國數學家傅里葉發現，任何周期函數都可以用正弦函數和餘弦函數構成的無窮級數來表示，這種表示方式就是傅里葉級數。假如有個波形比較復雜的周期函數，那麼找出能用來構成這個周期函數的正弦函數和餘弦函數的頻率的方法就叫做傅里葉變換，用這些頻率的正弦函數和餘弦函數疊加起來表示這個周期函數的方法就叫做傅里葉逆變換。
再從公式中看下他們的關系，首先介紹傅里葉級數到底是什麼，首先級數是指將數列的項依次用加號連接起來的函數。這么說可能大家還不理解，舉個例子：e^x=1+x/1!+x^2/2!+...x^n/n!....，等號左邊是指數函數，等號右邊就是級數。傅里葉級數公式如下：
我們主要看這個指數形式的傅里葉級數公式，把求和符號去掉，展開一下就是f(t)=Fa*e^jaω0t+Fb*e^jbω0t+Fc*e^jcω0t+Fd*e^jdω0.....。現在看下面的周期函數疊加效果圖，圖中顯示的是3個周期函數分別在坐標軸（橫軸時間，縱軸幅度）的圖像，寫成傅里葉級數形式就是f(t)=fa(t)+fb(t)+0+0....，這就是傅里葉級數公式要描述的現象。其中Fa*e^jaω0t=fa(t),Fb*e^jbω0t=fb(t),Fc*e^jcω0t=0....。
看下圖的傅里葉變換和逆變換公式，你會發現傅里葉逆變換公式和傅里葉級數公式極其相似，而傅里葉級數系數公式Fn又和傅里葉變換公式極其相似。所以對一個周期函數進行傅里葉級數展開的過程可以認為是先做傅里葉變換再做傅里葉逆變換的過程。
上圖就是傅里葉變換公式也叫連續傅里葉變換公式，有個很重要的事情，就是傅里葉變換公式和逆變換公式一定要一起給出，不然就會讓人誤解，你們在網上會看到各種各樣的寫法，但這些寫法都是對的，常見的如下圖所示。
為了方便後面的講解我把角頻率ω換成2πf，如上圖所示，ω是希臘字母讀作Omega，大寫是Ω，小寫是ω，以後這兩個字母會經常看到，都是等於2πf。不要和電學中的電阻單位搞混了，要明白字母只不過是一個符號而已，在不同學科領域都是混著用的，只要不和自己公式中其他字母沖突就行，例如上圖傅里葉變換公式中的j其實就是虛數單位i，一般時候我們會把虛數單位寫成i，但因為傅立葉變換經常用於電學解決一些問題，為了不和電流符號i混淆，所以公式就把i寫成j 。
要想了解傅里葉變換公式，首先要了解歐拉公式e^ix=cosx+isinx在圖像中的含義。以實部的值cosx作為橫坐標值，虛部sinx的值作為縱坐標值，x的取值從負無窮到正無窮，畫出所有的e^ix點後，你會發現這些點會形成一個周期為2π的圓。如下圖1所示（如果不理解，建議看3Blue1Brown的視頻，視頻連接：https://www.bilibili.com/video/BV1pW411J7s8）
所以歐拉公式e^ix其實就是隨著x的增大而在坐標繫上逆時針畫圓的過程，那麼e^-ix就表示順時針畫圓，e^-i2πx就表示畫圓的速度提高2π倍，也就是說x從0到1的過程就是順時針畫出一個完整圓的過程（當然x從1到2或者2到3等等，都能畫出一個完整的圓），把x換成t後，e^-i2πt表示每秒都會順時針畫出一個圓。e^-i2πft表示每秒都會順時針畫出f個圓。f(t)表示t時刻的振幅，f(t)函數畫出來就是時域波形圖。f(t)*e^-i2πft表示每經過1秒會順時針畫出f個圓，並在畫圓的同時，t時刻的圓半徑要乘上t時刻的振幅，其實就是以每秒的音頻振幅數據繞f圈的速度進行旋轉纏繞（為了方便理解，沒有用復雜的音頻數據，用的是一個頻率為3的正弦波音頻做的實驗，請看下圖2，圖的上半部分是時域波形圖，圖的左下角是f等於0.4的時候，用公式f(t)*e^-i2πft在實部和虛部構成的坐標系畫的圖，圖的右下角是頻譜圖，頻譜圖的橫坐標是頻率，縱坐標是振幅，振幅的值就是左下角圖中數據形成的圖案的質心（圖中的紅點）到坐標系原點的距離的2倍）。當改變f的值，你會發現數據大多數時候是和我們想的一樣，以坐標系原點為圓心環繞著，也就是振幅一直都是0，但是當f的值，也就每秒的圈數等於該音頻數據的頻率時，你會發現一個神奇的現象，那就是所有的數據會在實部或虛部坐標軸的一側形成一個圓（如下圖3所示，如此一來就知道這段音頻數據包含了一個頻率為3振幅為0.5的正弦波）。所以將多個正弦波疊加的音頻數據用傅里葉公式，f從負無窮到正無窮遍歷一遍，就可以把這個音頻數據里包含的正弦波都一一找出來。（如果不理解，建議看3Blue1Brown的視頻，視頻連接：https://www.bilibili.com/video/BV1pW411J7s8）
平時我們說的對音頻進行傅里葉變換處理，其實說的是短時離散傅里葉變換。短時離散傅里葉變換的公式（也可以直接叫做離散傅里葉變換公式）如下。
下面將教大家如何理解這個公式。上面說的連續傅里葉變換公式中有兩個原因導致我們無法使用，第一點要求是音頻數據的時間從負無窮到正無窮，第二點要求是任意時間t都要有幅度值x(t)才能代入公式進行計算。所以為了解決這兩個問題，把公式變為短時且離散的傅里葉變換公式，這個公式可以把一段時間（時間假設為Ts秒）的離散音頻數據（有N個采樣數據）進行傅里葉變換。你可以把離散傅里葉變換公式理解成連續傅里葉變換的變形，最重要的一點是連續傅里葉變換公式的f和離散傅里葉變換公式的k不是一個意思，他們的關系是k=f*Ts。所以離散傅里葉變換公式也可以寫成F(f)=1/n*∑f(t)*e^-j2πf*Ts*n/N，其中的Ts*n/N對應的就是連續傅里葉變換公式的t，只不過這個t沒辦法取任意時間了，t的取值也就隨著n的取值成為了離散的時間點，所以前面的系數由1/2π變為1/N。這樣這兩個公式就對應起來了。下面將進一步詳細介紹這個公式。
上一段說了k=f*Ts，這段我來解釋下為什麼，其實離散傅里葉變換公式中k表示的是這段Ts秒的音頻數據環繞坐標系原點的圈數，所以k並不是連續傅里葉變換公式里的頻率f，而頻率f指的是1秒鍾震盪的次數，在這個公式中頻率f也對應著1秒的音頻數據環繞的圈數，所以真正的頻率f=k/Ts。
有人可能會好奇，那為什麼不把離散傅里葉變換公式的自變數k換成f呢，這樣不是更好理解嗎？是會更好理解，但是沒有必要，用f的話還要做一次無用的換算。因為采樣點只有N個的原因，k的取值范圍就被限制住了，k的取值范圍只能是0~N-1的整數，這也是為什麼用k來做自變數而不是用f的原因。
還有人可能會好奇，傅里葉逆變換到底是怎麼把頻域的信息還原回時域的，其實公式計算出來的F(k)是一個復數，這個復數包含了這個頻率的周期函數的振幅和相位的信息，假設F(k)=a+ib，，F(k)的模|F(k)|=(a^2+b^2)^1/2，頻率f=k/Ts時的振幅為|F(k)|*2（因為求出來的值相當於圓心，但實際上振幅是圓離圓心最遠點到坐標原點的距離，所以要乘2），頻率f=k/Ts時的相位為arctan(b/a)。所以如果你知道一個周期函數包含了哪些頻率的周期函數，並且你這到這些周期函數的振幅和相位，你就可以像下圖一樣把fa(t)和fb(t)疊加在一起還原回f(t)。傅里葉逆變換的做法略有不同，但意思就是這么個意思，理解了離散傅里葉變換公式的計算，逆變換其實也是差不多代入數值計算就是了。（如果不理解怎麼用離散傅里葉變換公式計算，建議看視頻，視頻里有離散傅里葉變換完整的計算過程，視頻連接：https://www.hu.com/zvideo/1276595628009377792）
快速傅里葉變換推薦看下面兩個視頻
https://www.bilibili.com/video/BV1za411F76U
https://www.bilibili.com/video/BV1Jh411d7CN
下面是我用java實現的離散傅里葉變換及逆變換和快速傅里葉變換及逆變換，從他們的運行時間就可以看出來快速傅里葉變換快得多。（學完快速傅里葉變換再想想頻譜為何Y軸對稱？為何N/2對稱？）
❻ 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，得動態頻譜。。。。看你的需要。
❼ 求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 中
//引用任何所需的附加頭文件，而不是在此文件中引用