NOW ! It's time to implement "FFT" !

Fast Fourier Transform

              半年前的样子写了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

of 4 points. This pattern continues until there are N signals composed of a single point. An interlaced decomposition is used each time a signal is broken in two, that is, the signal is separated into its even and odd numbered samples.

上面 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 file       : bit_reverse.m
%     Code date       : 2014.09.04
%     e-mail          : jasonleaster@gmail.com
%
%     Function        : ret = get_r_in_Wn(k,m,N,bits)
%
%     Input Parameter : @k, the location of input signal
%                       @m, the current layyer
%                       @N, the total number of inputed signal
%                       @bits, how many bits should be used to represent 'N'
%
%     Output Parameter: @ret , the value of 'r'
%             
%****************************************************************************

function ret = get_r_in_Wn(k,m,N,bits)

         tmp = bitshift(k,bits - m);
         
         ret = bitand(tmp,(2^m) - 1);    
end

%****************************************************************************
%     Code writer     : EOF
%     Code file       : bit_reverse.m
%     Code date       : 2014.09.04
%     e-mail          : jasonleaster@gmail.com
%
%     Function        : ret = bit_reverse(Number,Bits)
%
%     Input Parameter : @Number, that is which number you want to reverse
%                                in bits.
%                       @Bits,   You should tell this function that how
%                                many bits you used to represent 'Number'
%
%     Output Parameter : @ret, that is a bit-reverse number of input 'Number'
%
%     Usage:
%            bit_reverse(1,3) == 4        001 --> 100 
%
%****************************************************************************
function ret = bit_reverse(Number,Bits)

    copy_num = Number;
    tmp_bits = Bits;
    tmp	 = 0;
	
    while	tmp_bits > 0
      
        tmp = tmp + (mod(copy_num,2)) * (2^(tmp_bits-1));
        
        copy_num = floor(copy_num/2);
        
        tmp_bits = tmp_bits - 1;
    end

    ret = tmp;
    
end

%% ***************************************************************************************
% 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