半年前的样子写了DFT，就想自己实现一个FFT，当时写不出，蝶形算法就是不明白,现在时机成熟了.
Now, it‘s time to implement "FFT".
How the FFT works
The FFT is a complicated algorithm, and its details are usually left to those that specialize in such things. This section describes the general operation of the FFT. Don‘t
worry if the details elude you; few scientists and engineers that use the FFT could write the program from scratch.
In complex notation, the time and frequency domains each contain one signal made up of N complex points. Each of these complex points is composed of two numbers, the real
part and the imaginary part.
这里输入序列是0~15的顺序输入
In this example, a 16 point signal is decomposed through four separate stages. The first stage breaks the 16 point signal into two signals each consisting of 8 points. The second stage decomposes the data into four signals
上面 0 ~ 15 的顺序输入被化做 0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15的输入，有什么依据可循吗？
看下面：
Bit reversal
Bit reversal的实现是很简单的
这里我给出C 的demo
/*******************************************************
Code writer : EOF
Code file : reverse_bits.c
e-mail : jasonleaster@gmail.com
code purpose:
This code is a demo for how to bit-reversal input bits
in FFT.
If there is something wrong with my code, please
touche me by e-mail. Thank you.
*********************************************************/
#include<stdio.h>
#define ARRAY_SIZE 16
int main(void)
{
int array[ARRAY_SIZE] = {0,};
int tmp = 0;
int bits = 0;
/*
** Initialization of array
*/
for(tmp = 0;tmp < ARRAY_SIZE;tmp++)
{
array[tmp] = tmp;
}
for(tmp = ARRAY_SIZE-1;tmp > 0; tmp >>= 1)
{
bits++;
}
for(tmp = 0;tmp < ARRAY_SIZE; tmp++)
{
printf("%4d %4d\n",array[tmp],reverse_bits(array[tmp],bits));
}
return 0 ;
}
/*
** Reverse the bits of input number
*/
int reverse_bits(int num,int bits)
{
int ret = 0;
int copy_num = 0;
for(ret = 0,copy_num = num; bits > 0; bits--)
{
ret += (copy_num % 2) * (1<<(bits-1));
copy_num >>= 1;
}
return ret;
}
在Octave下实现的FFT(没有利用特殊的API，所以应该可以在Matlab中运行，我没有Matlab环境，如果有心人可以去测试一下)
%% ***************************************************************************************
% code writer : EOF
% code date : 2014.09.03
% e-mail : jasonleaster@gmail.com
% code file : fft_EOF.m
% Version : 1.0
%
% code purpose :
%
% It's time to finish my demo for DFT. I would like to share my code with
% someone who is interesting in DSP. If there is something wrong with my code,
% please touche me by e-mail. Thank you!
%
%% ***************************************************************************************
clear all;
%*************************************************
% The number of all the signal that our sensor got
%*************************************************
TotalSample = 8;
% We assume that the preiod of the signal we generated is 'circle';
circle = TotalSample/2;
%**************************************************************
% This varible is used for recording the signal which
% were processed by inverse-DFT in time domain
%**************************************************************
SignalInT = zeros(TotalSample,1);
SignalInT_reversed = zeros(TotalSample,1);
%This varible is used for recording the signal which were processed by inverse-DFT in time domain
OutPutSignal = zeros(TotalSample,1);
OriginalSignal = zeros(TotalSample,1);
%This varible is used for recording the original signal that we got.
%% initialize a square wave
for SampleNumber = -(TotalSample/2):(TotalSample/2)-1
if (mod(abs(SampleNumber),circle) < (circle/2))&&(SampleNumber>0)
OriginalSignal((TotalSample/2)+1+SampleNumber) = 5;
elseif (mod(abs(SampleNumber),circle) >= (circle/2))&&(SampleNumber>0)
OriginalSignal((TotalSample/2)+1+SampleNumber) = 0;
elseif (mod(abs(SampleNumber),circle) < (circle/2))&&(SampleNumber<0)
OriginalSignal((TotalSample/2)+1+SampleNumber) = 0;
elseif (mod(abs(SampleNumber),circle) >= (circle/2))&&(SampleNumber<0)
OriginalSignal((TotalSample/2)+1+SampleNumber) = 5;
end
end
% We show the original signal in time domain.
figure(1);
plot( -(TotalSample/2):(TotalSample/2)-1,OriginalSignal,'.-');
title('The original signal');
tmp = TotalSample - 1;
%%***********************************************************************
% @Bits : describe how many bits should be used to make up the TotalSample
%%***********************************************************************
Bits = 0;
while tmp > 0
%% floor (X) Return the largest integer not greater than X.
tmp = floor(tmp/2);
Bits = Bits + 1;
end
SignalInT = OriginalSignal;
%******************************************
% Reverse the bits of input number
%******************************************
for SampleNumber = 1 : TotalSample
ret = bit_reverse(SampleNumber - 1,Bits);
SignalInT_reversed(SampleNumber) = SignalInT(ret+1);
end
InPutSignal = SignalInT_reversed;
for layyer = 1 : Bits
for SampleNumber = 1 : 2^(layyer) : TotalSample
for pre_half = SampleNumber:(SampleNumber+2^(layyer-1) -1)
r = get_r_in_Wn(pre_half-1,layyer,TotalSample,Bits);
W_rN = exp(-2*pi*j*(r)/TotalSample) ;
OutPutSignal(pre_half) = ...
InPutSignal(pre_half) + ...
W_rN * InPutSignal(pre_half + 2^(layyer-1));
OutPutSignal(pre_half + 2^(layyer-1)) = ...
InPutSignal(pre_half) - ...
W_rN * InPutSignal(pre_half + 2^(layyer-1));
end
end
InPutSignal = OutPutSignal;
end程序的运行结果输出数据OutPutSignal,和Octave自带的FFT函数对比，发现计算结果是正确的 \O/
这个blog仅仅记录FFT的实现结果，之后将补充理论基础以及实现细节分析.
—— EOF
于XTU. 2014.09.04
NOW ! It's time to implement "FFT" !
原文地址：http://blog.csdn.net/cinmyheart/article/details/39042623
