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

Fast Fourier Transform

Half a year ago, I wrote DFT and wanted to implement an FFT myself. At that time, I couldn't write it, and the butterfly algorithm just didn't understand it. Now the time is ripe.

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 cocould 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.

An important operation-backward Input

The input sequence is 0 ~ 15. sequential input

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.

Above 0 ~ In the order of 15, the input is converted into 0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15. Is there any evidence to follow?

See the following:

Bit reversal

The implementation of bit reversal is very simple.

Here is a demo of C.

/*******************************************************Code writer: EOFCode file        : reverse_bits.ce-mail : [email protected]code purpose:This code is a demo for how to bit-reversal input bitsin FFT.If there is something wrong with my code, pleasetouche me by e-mail. Thank you.*********************************************************/#include<stdio.h>#define ARRAY_SIZE 16int 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;}

FFT implemented under Ave ave (no special API is used, so it should be able to run in MATLAB. I don't have a MATLAB environment. If you are interested, you can test it)

%% ***************************************************************************************% code writer : EOF% code date     : 2014.09.03% e-mail       : [email protected]% 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 domainOutPutSignal = zeros(TotalSample,1);OriginalSignal = zeros(TotalSample,1);%This varible is used for recording the original signal that we got.%% initialize a square wavefor 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;    endend%  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;whiletmp > 0    %%  floor (X) Return the largest integer not greater than X.    tmp = floor(tmp/2);    Bits = Bits + 1;endSignalInT = OriginalSignal;%******************************************%   Reverse the bits of input number%******************************************for SampleNumber = 1 : TotalSample        ret = bit_reverse(SampleNumber - 1,Bits);        SignalInT_reversed(SampleNumber) = SignalInT(ret+1);endInPutSignal = 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, the output data of the program running result, is compared with the FFT function provided by octave, and it is found that the calculation result is correct \ O/

This blog only records the implementation results of FFT, and will supplement the theoretical basis and implementation details analysis.

-- EOF

On xtu. 2014.09.04

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