fft演算法代碼

① matlab中fft2的計算代碼！！！

function f = fft2(x, mrows, ncols)
%FFT2 Two-dimensional discrete Fourier Transform.
% FFT2(X) returns the two-dimensional Fourier transform of matrix X.
% If X is a vector, the result will have the same orientation.
%
% FFT2(X,MROWS,NCOLS) pads matrix X with zeros to size MROWS-by-NCOLS
% before transforming.
%
% Class support for input X:
% float: double, single
%
% See also FFT, FFTN, FFTSHIFT, FFTW, IFFT, IFFT2, IFFTN.

% Copyright 1984-2010 The MathWorks, Inc.
% $Revision: 5.13.4.3 $ $Date: 2010/02/25 08:08:32 $

if ismatrix(x)
if nargin==1
f = fftn(x);
else
f = fftn(x,[mrows ncols]);
end
else
if nargin==1
f = fft(fft(x,[],2),[],1);
else
f = fft(fft(x,ncols,2),mrows,1);
end
end

② fft演算法matlab的實現代碼！完整版的！

function result = MyFFT(vector)
result = fft(vector);

③ 基2—fft演算法的軟體實現(MATLAB代碼)

參考網路： clc; clear all; close all; x=ones(1,128); %輸入的信號，自己可以改變 %整體運用原位計算 m=nextpow2(x);N=2^m; % 求x的長度對應的2的最低冪次m if length(x)<N x=[x,zeros(1,N-length(x))]; % 若x的長度不是2的冪，補零到2的整數冪 end nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; % 求1:2^m數列序號的倒序 y=x(nxd); % 將x倒序排列作為y的初始值 for mm=1:m % 將DFT作m次基2分解,從左到右，對每次分解作DFT運算,共做m級蝶形運算，每一級都有2^(mm-1)個蝶形結 Nz=2^mm;u=1; % 旋轉因子u初始化為WN^0=1 WN=exp(-i*2*pi/Nz); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nz) for j=1:Nz/2 % 本次跨越間隔內的各次蝶形運算,在進行第mm級運算時需要2^(mm-1)個 蝶形 for k=j:Nz:N % 本次蝶形運算的跨越間隔為Nz=2^mm kp=k+Nz/2; % 蝶形運算的兩個因子對應單元下標的關系 t=y(kp)*u; % 蝶形運算的乘積項 y(kp)=y(k)-t; % 蝶形運算 y(k)=y(k)+t; % 蝶形運算 end u=u*WN; % 修改旋轉因子,多乘一個基本DFT因子WN end end y y1=fft(x) %自己編的FFT跟直接調用的函數運算以後的結果進行對比

④ 一維復數序列的快速傅里葉變換（FFT）

設x（N）為N點有限長離散序列，代入式（8－3）、式（8－4），並令 其傅里葉變換（DFT）為

地球物理數據處理基礎

反變換（IDFT）為

地球物理數據處理基礎

兩者的差異只在於W的指數符號不同，以及差一個常數1/N，因此下面我們只討論DFT正變換式（8－5）的運算量，其反變換式（8－6）的運算是完全相同的。

一般來說，W是復數，因此，X（j）也是復數，對於式（8－5）的傅里葉變換（DFT），計算一個X（j）值需要N次復數乘法和N－1次復數加法。而X（j）一共有N個值（j＝0，1，…，N－1），所以完成整個DFT運算總共需要N²次復數乘法和N（N－1）次復數加法。

直接計算DFT，乘法次數和加法次數都是與N²成正比的，當N很大時，運算量會很大，例如，當N＝8時，DFT需64次復數乘法；而當N＝1024時，DFT所需乘法為1048576次，即一百多萬次的復數乘法運算，對運算速度要求高。所以需要改進DFT的計算方法，以減少運算次數。

分析W^jk，表面上有N²個數值，由於其周期性，實際上僅有N個不同的值W⁰，W¹，…，W^N－1。對於N＝2^m時，由於其對稱性，只有N/2個不同的值W⁰，W¹，…，

地球物理數據處理基礎

因此可以把長序列的DFT分解為短序列DFT，而前面已經分析DFT與N²成正比，所以N越小越有利。同時，利用ab＋ac＝a（b＋c）結合律法則，可以將同一個W^r對應的系數x（k）相加後再乘以W^r，就能大大減少運算次數。這就是快速傅里葉變換（FFT）的演算法思路。

下面，我們來分析N＝2^m情況下的FFT演算法。

1.N＝4的FFT演算法

對於m＝2，N＝4，式（8－5）傅里葉變換為

地球物理數據處理基礎

將式（8－7）寫成矩陣形式

地球物理數據處理基礎

為了便於分析，將上式中的j，k寫成二進制形式，即

地球物理數據處理基礎

代入式（8－7），得

地球物理數據處理基礎

分析W^jk的周期性來減少乘法次數

地球物理數據處理基礎

則 代回式（8－9），整理得

地球物理數據處理基礎

上式可分層計算，先計算內層，再計算外層時就利用內層計算的結果，可避免重復計算。寫成分層形式

地球物理數據處理基礎

則X（j₁ j₀）＝X₂（j₁ j₀）。

上式表明對於N＝4的FFT，利用W^r的周期關系可分為m＝2步計算。實際上，利用W^r的對稱性，仍可以對式（8－11）進行簡化計算。考慮到

地球物理數據處理基礎

式（8－11）可以簡化為

地球物理數據處理基礎

令j＝j₀；k＝k₀，並把上式表示為十進制，得

地球物理數據處理基礎

可以看到，完成上式N＝4的FFT計算（表8－1）需要N·（m－1）/2＝2次復數乘法和N·m＝8次復數加法，比N＝4的DFT演算法的N²＝16次復數乘法和N·（N－1）＝12次復數加法要少得多。

表8－1 N＝4的FFT演算法計算過程

註：W⁰＝1；W¹＝－i。

［例1］求N＝4樣本序列1，3，3，1的頻譜（表8－2）。

表8－2 N＝4樣本序列

2.N＝8的FFT演算法

類似N＝4的情況，用二進制形式表示，有

地球物理數據處理基礎

寫成分層計算的形式：

地球物理數據處理基礎

則X（j₂ j₁ j₀）＝X₃（j₂ j₁ j₀）。

對式（8－14）的X₁（k₁ k₀ j₀）進行展開，有

地球物理數據處理基礎

還原成十進制，並令k＝2k₁＋k₀，即k＝0，1，2，3，有

地球物理數據處理基礎

用類似的方法對式（8－14）的X₂（k₀ j₁ j₀），X₃（j₂ j₁ j₀）進行展開，整理得

地球物理數據處理基礎

用式（8－16）、式（8－17）逐次計算到X₃（j）＝X（j）（j＝0，1，…，7），即完成N＝2³＝8的FFT計算，其詳細過程見表8－3。

表8－3 N＝8的FFT演算法計算過程

註：對於正變換 對於反變換 所

［例2］求N＝8樣本序列（表8－4）x（k）＝1，2，1，1，3，2，1，2的頻譜。

表8－4 N＝8樣本序列

3.任意N＝2^m的FFT演算法

列出N＝4，N＝8的FFT計算公式，進行對比

地球物理數據處理基礎

觀察式（8－18）、式（8－19），不難看出，遵循如下規律：

（1）等式左邊的下標由1遞增到m，可用q＝1，2，…，m代替，則等式右邊為q－1；

（2）k的上限為奇數且隨q的增大而減小，至q＝m時為0，所以其取值范圍為k＝0，1，2，…，（2^m－q－1）；

（3）j的上限為奇數且隨q的增大而增大，且q＝1時為0，其取值范圍為j＝0，1，2，…，（2^q－1－1）；

（4）k的系數，在等式左邊為2^q，等式右邊為2^q－1（包括W的冪指數）；

（5）等式左邊序號中的常數是2的乘方形式，且冪指數比下標q小1，即2^q－1；等式右邊m對式子序號中的常數都是定值2^m－1。

歸納上述規則，寫出對於任意正整數m，N＝2^m的FFT演算法如下：

由X₀（p）＝x（p）（p＝0，1，…，N－1）開始：

（1）對q＝1，2，…，m，執行（2）～（3）步；

（2）對k＝0，1，2，…，（2^m－q－1）及j＝0，1，2，…，（2^q－1－1），執行

地球物理數據處理基礎

（3）j，k循環結束；

（4）q循環結束；由X_m（p）（p＝0，1，…，N－1）輸出原始序列x（p）的頻譜X（p）。

在計算機上很容易實現上述FFT演算法程序，僅需要三個復數數組，編程步驟如下：

（1）設置復數數組X₁（N－1），X₂（N－1）和 （數組下界都從0開始）；

（2）把樣本序列x賦給X₁，即X₁（k）＝x（k）（k＝0，1，…，N－1）；

（3）計算W，即正變換 反變換

（4）q＝1，2，…，m，若q為偶數，執行（6），否則執行第（5）步；

（5）k＝0，1，2，…，（2^m－q－1）和j＝0，1，2，…，（2^q－1－1）循環，作

X₂（2^qk＋j）＝X₁（2^q－1k＋j）＋X₁（2^q－1k＋j＋2^m－1）

X₂（2^qk＋j＋2^q－1）＝［X₁（2^q－1k＋j）－X₁（2^q－1k＋j＋2^m－1）］W（2^q－1k）

至k，j循環結束；

（6）k＝0，1，2，…，（2^m－q－1）和j＝0，1，2，…，（2^q－1－1）循環，作

X₁（2^qk＋j）＝X₂（2^q－1k＋j）＋X₂（2^q－1k＋j＋2^m－1）

X₁（2^qk＋j＋2^q－1）＝［X₂（2^q－1k＋j）－X₂（2^q－1k＋j＋2^m－1）］W（2^q－1k）

至k，j循環結束；

（7）q循環結束，若m為偶數，輸出X₁（j），否則輸出X₂（j）（j＝0，1，…，N－1），即為所求。

⑤ 求用C++實現的FFT演算法

很早以前的，如果管用別忘了給我加分呀
/*
This computes an in-place complex-to-complex FFT
x and y are the real and imaginary arrays of 2^m points.
dir = 1 gives forward transform
dir = -1 gives reverse transform
*/
short FFT(short int dir,long m,double *x,double *y)
{
long n,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;

/* Calculate the number of points */
n = 1;
for (i=0;i<m;i++)
n *= 2;

/* Do the bit reversal */
i2 = n >> 1;
j = 0;
for (i=0;i<n-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}

/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<n;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}

/* Scaling for forward transform */
if (dir == 1) {
for (i=0;i<n;i++) {
x[i] /= n;
y[i] /= n;
}
}

return(TRUE);
}

---------------------------------------------------------------------------------
/*****************fft programe*********************/
#include "typedef.h"
#include "math.h"

struct compx EE(struct compx b1,struct compx b2)
{
struct compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}

void FFT(struct compx *xin,int N)
{
int f,m,nv2,nm1,i,k,j=1,l;
/*int f,m,nv2,nm1,i,k,j=N/2,l;*/
struct compx v,w,t;
nv2=N/2;
f=N;
for(m=1;(f=f/2)!=1;m++){;}
nm1=N-1;

/*變址運算*/
for(i=1;i <=nm1;i++)
{
if(i <j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
k=nv2;
while(k <j){j=j-k;k=k/2;}
j=j+k;
}

{
int le,lei,ip;
float pi;
for(l=1;l <=m;l++)
{ le=pow(2,l);// 這里用的是L而不是1 !!!!
lei =le/2;
pi=3.14159;
v.real=1.0;
v.imag=0.0;
w.real=cos(pi/lei);
w.imag=-sin(pi/lei);
for(j=1;j <=lei;j++)
{
/*double p=pow(2,m-l)*j;
double ps=2*pi/N*p;
w.real=cos(ps);
w.imag=-sin(ps);*/
for(i=j;i <=N;i=i+le)
{ /* w.real=cos(ps);
w.imag=-sin(ps);*/
ip=i+lei;
t=EE(xin[ip],v);
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
v=EE(v,w);
}
}
}
return;
}

/*****************main programe********************/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "typedef.h"

float result[257];
struct compx s[257];
int Num=256;
const float pp=3.14159;

main()
{

int i=1;
for(;i <0x101;i++)
{
s[i].real=sin(pp*i/32);
s[i].imag=0;
}

FFT(s,Num);

for(i=1;i <0x101;i++)
{
result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
}

}

-----------------------------------------------------------------------------------
FFT變換 C源代碼

FFT C source code (Simple radix-2)

void fft_float (
unsigned NumSamples,
int InverseTransform,
float *RealIn,
float *ImagIn,
float *RealOut,
float *ImagOut )
{
unsigned NumBits; /* Number of bits needed to store indices */
unsigned i, j, k, n;
unsigned BlockSize, BlockEnd;
double angle_numerator = 2.0 * DDC_PI;
double tr, ti; /* temp real, temp imaginary */
if ( !IsPowerOfTwo(NumSamples) )
{
fprintf (
stderr,
"Error in fft(): NumSamples=%u is not power of two\n",
NumSamples );
exit(1);
}
if ( InverseTransform )
angle_numerator = -angle_numerator;
CHECKPOINTER ( RealIn );
CHECKPOINTER ( RealOut );
CHECKPOINTER ( ImagOut );
NumBits = NumberOfBitsNeeded ( NumSamples );
/*
** Do simultaneous data and bit-reversal ordering into outputs...
*/
for ( i=0; i < NumSamples; i++ )
{
j = ReverseBits ( i, NumBits );
RealOut[j] = RealIn;
ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn;
}
/*
** Do the FFT itself...
*/
BlockEnd = 1;
for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
{
double delta_angle = angle_numerator / (double)BlockSize;
double sm2 = sin ( -2 * delta_angle );
double sm1 = sin ( -delta_angle );
double cm2 = cos ( -2 * delta_angle );
double cm1 = cos ( -delta_angle );
double w = 2 * cm1;
double ar[3], ai[3];
double temp;
for ( i=0; i < NumSamples; i += BlockSize )
{
ar[2] = cm2;
ar[1] = cm1;
ai[2] = sm2;
ai[1] = sm1;
for ( j=i, n=0; n < BlockEnd; j++, n++ )
{
ar[0] = w*ar[1] - ar[2];
ar[2] = ar[1];
ar[1] = ar[0];
ai[0] = w*ai[1] - ai[2];
ai[2] = ai[1];
ai[1] = ai[0];
k = j + BlockEnd;
tr = ar[0]*RealOut[k] - ai[0]*ImagOut[k];
ti = ar[0]*ImagOut[k] + ai[0]*RealOut[k];
RealOut[k] = RealOut[j] - tr;
ImagOut[k] = ImagOut[j] - ti;
RealOut[j] += tr;
ImagOut[j] += ti;
}
}
BlockEnd = BlockSize;
}
/*
** Need to normalize if inverse transform...
*/
if ( InverseTransform )
{
double denom = (double)NumSamples;
for ( i=0; i < NumSamples; i++ )
{
RealOut /= denom;
ImagOut /= denom;
}
}
}

int IsPowerOfTwo ( unsigned x )
{
if ( x < 2 )
return FALSE;
if ( x & (x-1) ) // Thanks to 'byang' for this cute trick!
return FALSE;
return TRUE;
}

unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo )
{
unsigned i;
if ( PowerOfTwo < 2 )
{
fprintf (
stderr,
">>> Error in fftmisc.c: argument %d to NumberOfBitsNeeded is too small.\n",
PowerOfTwo );
exit(1);
}
for ( i=0; ; i++ )
{
if ( PowerOfTwo & (1 << i) )
return i;
}
}

unsigned ReverseBits ( unsigned index, unsigned NumBits )
{
unsigned i, rev;
for ( i=rev=0; i < NumBits; i++ )
{
rev = (rev << 1) | (index & 1);
index >>= 1;
}
return rev;
}

double Index_to_frequency ( unsigned NumSamples, unsigned Index )
{
if ( Index >= NumSamples )
return 0.0;
else if ( Index <= NumSamples/2 )
return (double)Index / (double)NumSamples;
return -(double)(NumSamples-Index) / (double)NumSamples;
}

⑥ FFT的公式是什麼和演算法是怎樣實現

二維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演算法的實現。下面討論一維FFT的演算法原理。

【1D-FFT的演算法實現】
設序列h(n)長度為N，將其按下標的奇偶性分成兩組，即he和ho序列，它們的長度都是N/2。這樣，可以將h(n)的FFT計算公式改寫如下 ：

(A)
由於

所以，(A)式可以改寫成下面的形式：

按照FFT的定義，上面的式子實際上是：

其中，k的取值范圍是 0~N-1。
我們注意到He(k)和Ho(k)是N/2點的DFT，其周期是N/2。因此，H(k)DFT的前N/2點和後N/2點都可以用He(k)和Ho(k)來表示

⑦ FFT 演算法 C編程 求救~~

#include<graphics.h>
#include<stdio.h>
#include<math.h>
#include<conio.h>
#include<dos.h>
#define pi 3.14159

float ar[512],ai[512];
float a[1025];

void testdata()
{
int i;
for(i=0;i<512;i++)
{
ar[i]=50*cos(i*pi/32)+60*cos(i*pi/16)+53*cos(i*pi/64)+24*cos(i*pi/8)+10*cos(i*pi/128);
ai[i]=0;
}
}

void fft(int nn,int t)
{
int n1,n2,i,j,k,l,m,s,l1;
float t1,t2,x,y;
float w1,w2,u1,u2,z;

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=cos(pi/k);
w2=t*sin(pi/k);
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]);
}
}

int main ()
{
int i;
int gdriver=DETECT,gmode;
initgraph(&gdriver,&gmode,"");
setfillstyle(SOLID_FILL,WHITE);
bar(0,0,639,479);
//模擬測試數據
testdata();
//波形顯示
setcolor(BLACK);
line(10,119,550,119); // x軸
line(10,10,10,239); // y軸
setcolor(BLUE);
for(i=0;i<511;i++)
line(i+10,119-ar[i],i+10+1,119-ar[i+1]);
//FFT
fft(512,-1);

//頻譜顯示
setcolor(BLACK);
line(10,459,550,459); // x軸
line(10,259,10,459); // y軸
setcolor(RED);
for(i=0;i<256;i++)
line(2*i+10,459-a[i],2*i+10,459);
getch();
closegraph();
return 0;
}

⑧ 基於FFT的演算法優化 要C語言完整程序（利用旋轉因子的性質），有的請留言，答謝！！！(有核心代碼，望指教

實現(C描述)

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

//#include "complex.h"

// --------------------------------------------------------------------------

#define N 8 //64

#define M 3 //6 //2^m=N

#define PI 3.1415926

// --------------------------------------------------------------------------

float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};

float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};

float x_i[N]; //N=8

/*

float twiddle[N/2] = {1, 0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,

0.7075, 0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,

0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,

-0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951}; //N=64

float x_r[N]={1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,};

float x_i[N];

*/

FILE *fp;

// ----------------------------------- func -----------------------------------

/**

* 初始化輸出虛部

*/

static void fft_init( void )

{

int i;

for(i=0; i<N; i++) x_i[i] = 0.0;

}

/**

* 反轉演算法.將時域信號重新排序.

* 這個演算法有改進的空間

*/

static void bitrev( void )

{

int p=1, q, i;

int bit_rev[ N ]; //

float xx_r[ N ]; //

bit_rev[ 0 ] = 0;

while( p < N )

{

for(q=0; q<p; q++)

{

bit_rev[ q ] = bit_rev[ q ] * 2;

bit_rev[ q + p ] = bit_rev[ q ] + 1;

}

p *= 2;

}

for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];

for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];

}

/* ------------ add by sshc625 ------------ */

static void bitrev2( void )

{

return ;

}

/* */

void display( void )

{

printf("\n\n");

int i;

for(i=0; i<N; i++)

printf("%f\t%f\n", x_r[i], x_i[i]);

}

/**

*

*/

void fft1( void )

{ fp = fopen("log1.txt", "a+");

int L, i, b, j, p, k, tx1, tx2;

float TR, TI, temp; // 臨時變數

float tw1, tw2;

/* 深M. 對層進行循環. L為當前層, 總層數為M. */

for(L=1; L<=M; L++)

{

fprintf(fp,"----------Layer=%d----------\n", L);

/* b的意義非常重大,b表示當前層的顆粒具有的輸入樣本點數 */

b = 1;

i = L - 1;

while(i > 0)

{

b *= 2;

i--;

}

// -------------- 是否外層對顆粒循環, 內層對樣本點循環邏輯性更強一些呢! --------------

/*

* outter對參與DFT的樣本點進行循環

* L=1, 循環了1次(4個顆粒, 每個顆粒2個樣本點)

* L=2, 循環了2次(2個顆粒, 每個顆粒4個樣本點)

* L=3, 循環了4次(1個顆粒, 每個顆粒8個樣本點)

*/

for(j=0; j<b; j++)

{

/* 求旋轉因子tw1 */

p = 1;

i = M - L; // M是為總層數, L為當前層.

while(i > 0)

{

p = p*2;

i--;

}

p = p * j;

tx1 = p % N;

tx2 = tx1 + 3*N/4;

tx2 = tx2 % N;

// tw1是cos部分, 實部; tw2是sin部分, 虛數部分.

tw1 = ( tx1>=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ];

tw2 = ( tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];

/*

* inner對顆粒進行循環

* L=1, 循環了4次(4個顆粒, 每個顆粒2個輸入)

* L=2, 循環了2次(2個顆粒, 每個顆粒4個輸入)

* L=3, 循環了1次(1個顆粒, 每個顆粒8個輸入)

*/

for(k=j; k<N; k=k+2*b)

{

TR = x_r[k]; // TR就是A, x_r[k+b]就是B.

TI = x_i[k];

temp = x_r[k+b];

/*

* 如果復習一下 (a+j*b)(c+j*d)兩個復數相乘後的實部虛部分別是什麼

* 就能理解為什麼會如下運算了, 只有在L=1時候輸入才是實數, 之後層的

* 輸入都是復數, 為了讓所有的層的輸入都是復數, 我們只好讓L=1時候的

* 輸入虛部為0

* x_i[k+b]*tw2是兩個虛數相乘

*/

fprintf(fp, "tw1=%f, tw2=%f\n", tw1, tw2);

x_r[k] = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;

x_i[k] = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;

x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;

x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;

fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k, x_r[k], x_i[k]);

fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f\n", k+b, x_r[k+b], x_i[k+b]);

} //

} //

} //

}

/**

* ------------ add by sshc625 ------------

* 該實現的流程為

* for( Layer )

* for( Granule )

* for( Sample )

*

*

*

*

*/

void fft2( void )

{ fp = fopen("log2.txt", "a+");

int cur_layer, gr_num, i, k, p;

float tmp_real, tmp_imag, temp; // 臨時變數, 記錄實部

float tw1, tw2;// 旋轉因子,tw1為旋轉因子的實部cos部分, tw2為旋轉因子的虛部sin部分.

int step; // 步進

int sample_num; // 顆粒的樣本總數(各層不同, 因為各層顆粒的輸入不同)

/* 對層循環 */

for(cur_layer=1; cur_layer<=M; cur_layer++)

{

/* 求當前層擁有多少個顆粒(gr_num) */

gr_num = 1;

i = M - cur_layer;

while(i > 0)

{

i--;

gr_num *= 2;

}

/* 每個顆粒的輸入樣本數N' */

sample_num = (int)pow(2, cur_layer);

/* 步進. 步進是N'/2 */

step = sample_num/2;

/* */

k = 0;

/* 對顆粒進行循環 */

for(i=0; i<gr_num; i++)

{

/*

* 對樣本點進行循環, 注意上限和步進

*/

for(p=0; p<sample_num/2; p++)

{

// 旋轉因子, 需要優化...

tw1 = cos(2*PI*p/pow(2, cur_layer));

tw2 = -sin(2*PI*p/pow(2, cur_layer));

tmp_real = x_r[k+p];

tmp_imag = x_i[k+p];

temp = x_r[k+p+step];

/*(tw1+jtw2)(x_r[k]+jx_i[k])

*

* real : tw1*x_r[k] - tw2*x_i[k]

* imag : tw1*x_i[k] + tw2*x_r[k]

* 我想不抽象出一個

* typedef struct {

* double real; // 實部

* double imag; // 虛部

* } complex; 以及針對complex的操作

* 來簡化復數運算是否是因為效率上的考慮!

*/

/* 蝶形演算法 */

x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );

x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );

/* X[k] = A(k)+WB(k)

* X[k+N/2] = A(k)-WB(k) 的性質可以優化這里*/

// 旋轉因子, 需要優化...

tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));

tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));

x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );

x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );

printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);

printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);

}

/* 開跳!:) */

k += 2*step;

}

}

}

/*

* 後記:

* 究竟是顆粒在外層循環還是樣本輸入在外層, 好象也差不多, 復雜度完全一樣.

* 但以我資質愚鈍花費了不少時間才弄明白這數十行代碼.

* 從中我發現一個於我非常有幫助的教訓, 很久以前我寫過一部分演算法, 其中絕大多數都是遞歸.

* 將數據量減少, 減少再減少, 用歸納的方式來找出數據量加大代碼的規律

* 比如FFT

* 1. 先寫死LayerI的代碼; 然後再把LayerI的輸出作為LayerII的輸入, 又寫死代碼; ......

* 大約3層就可以統計出規律來. 這和遞歸也是一樣, 先寫死一兩層, 自然就出來了!

* 2. 有的功能可以寫偽代碼, 不急於求出結果, 降低復雜性, 把邏輯結果定出來後再添加.

* 比如旋轉因子就可以寫死, 就寫1.0. 流程出來後再寫旋轉因子.

* 寥寥數語, 我可真是流了不少汗! Happy!

*/

void dft( void )

{

int i, n, k, tx1, tx2;

float tw1,tw2;

float xx_r[N],xx_i[N];

/*

* clear any data in Real and Imaginary result arrays prior to DFT

*/

for(k=0; k<=N-1; k++)

xx_r[k] = xx_i[k] = x_i[k] = 0.0;

// caculate the DFT

for(k=0; k<=(N-1); k++)

{

for(n=0; n<=(N-1); n++)

{

tx1 = (n*k);

tx2 = tx1+(3*N)/4;

tx1 = tx1%(N);

tx2 = tx2%(N);

if(tx1 >= (N/2))

tw1 = -twiddle[tx1-(N/2)];

else

tw1 = twiddle[tx1];

if(tx2 >= (N/2))

tw2 = -twiddle[tx2-(N/2)];

else

tw2 = twiddle[tx2];

xx_r[k] = xx_r[k]+x_r[n]*tw1;

xx_i[k] = xx_i[k]+x_r[n]*tw2;

}

xx_i[k] = -xx_i[k];

}

// display

for(i=0; i<N; i++)

printf("%f\t%f\n", xx_r[i], xx_i[i]);

}

// ---------------------------------------------------------------------------

int main( void )

{

fft_init( );

bitrev( );

// bitrev2( );

//fft1( );

fft2( );

display( );

system( "pause" );

// dft();

return 1;

}

本文來自CSDN博客，轉載請標明出處：http://blog.csdn.net/sshcx/archive/2007/06/14/1651616.aspx

⑨ MATLAB中fft的頻率軸怎麼計算

假設你的信號是 8個點，采樣頻率是 100Hz。

那麼，該信號的頻率是50Hz，那麼頻率軸每個間隔是 50/(8-1), 設為df
那麼，頻率軸是 0 df 2*df 3*df 4*df
也就是說，對於8個點的信號，你會得到頻率間隔是 50/(8-1), 可以得到 8/2+1個頻率點。
也就是說，對於N個點的信號，你會得到頻率間隔是 50/(N-1), 可以得到 N/2+1個頻率點。
注意，N是2的某次冪

⑩ 怎麼用C語言實現FFT演算法 呀

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]); /* 幅值 */
}
}

(http://..com/question/284943905.html?an=0&si=2)