Fixed-Point Toolbox 3.0

Fixed-Point Fast Fourier Transform (FFT)

This is an example of converting a textbook Fast Fourier Transform (FFT) algorithm into fixed-point MATLAB® code, and then into fixed-point C code.

There is a fixed-point FFT in the Signal Processing Blockset™ product, so we include that implementation at the end for comparison.

Contents

Development Process

We will follow these steps, with the emphasis on doing all of the fixed-point algorithm design and analysis in M before converting to fixed-point C. We will then use the fixed-point M to debug and validate the fixed-point C.

1. Implement textbook algorithm in M.

2. Verify with floating-point in M.

3. Identify fixed-point issues in M.

4. Modify algorithm in M to address fixed-point issues.

5. Convert fixed-point M to C.

6. Debug C code by comparing with fixed-point M using integer display.

FFT Definition

The FFT is a complex-valued linear transformation from the time domain to the frequency domain. For example, if we construct a vector as the sum of two sinusoids and transform it with the FFT, then we can see the peaks of the frequencies in the FFT magnitude plot.

n = 64;                                     % Number of points
Fs = 4;                                     % Sampling frequency in Hz
t  = (0:(n-1))/Fs;                          % Time vector
f  = linspace(0,Fs,n);                      % Frequency vector
f0 = .2; f1 = .5;                           % Frequencies, in Hz
x0 = cos(2*pi*f0*t) + 0.55*cos(2*pi*f1*t);  % Time-domain signal
y  = fft(x0);                               % Frequency-domain transformatio
n

figure(gcf); clf
subplot(211); plot(t,x0,'b.-'); xlabel('Time (s)'); ylabel('Amplitude');lege
nd('x0')
subplot(212); plot(f,abs(y),'m.-'); xlabel('Frequency (Hz)'); ylabel('Magnit
ude');legend('abs(fft(x0))')

The peaks at 0.2 and 0.5 Hz in the frequency plot correspond to the two sinusoids of the time-domain signal at those frequencies.

Note the reflected peaks at 3.5 and 3.8 Hz. When the input to an FFT is real-valued as it is in this case, the output y is conjugate-symmetric:

The Algorithm Pseudo-Code

There are many implementations of the FFT with competing costs and benefits. You may find that a different algorithm is better for your application than the one given here. This algorithm is given as an example of how you might begin your own exploration.

Suppose that we're reading the book Computational Frameworks for the Fast Fourier Transform by Charles Van Loan http://www.mathworks.com/support/books/book1384.html and we're interested in writing a fixed-point version of the decimation-in-time unit-stride FFT in Algorithm 1.6.2 on page 45.

In pseudo-code, the algorithm in the textbook is as follows.

In the textbook algorithm, indexing is zero-based, Fn is an n-by-n Fourier-transform matrix, Pn is an n-by-n bit-reversal permutation matrix, and w is a complex vector of twiddle factors. The twiddle factors w are complex roots of unity computed by the following algorithm:

type_nocomments fi_radix2twiddles

function w = fi_radix2twiddles(n)
t = log2(n);
if floor(t) ~= t
  error('N must be an exact power of two.');
end
w = zeros(n-1,1);
k=1;
L=2;
while L<=n
  theta = 2*pi/L;
  for j=0:(L/2 - 1)
    w(k) = complex( cos(j*theta), -sin(j*theta) );
    k = k + 1;
  end
  L = L*2;
end

figure(gcf);clf
w0 = fi_radix2twiddles(n);
polar(angle(w0),abs(w0),'o')
title('Twiddle Factors: Complex roots of unity')

The Algorithm in M

In MATLAB, we write the algorithm as follows, using the bitreverse function, and adding 1 to the indices to make them 1-based.

type_nocomments fi_m_radix2fft_algorithm1_6_2

function x = fi_m_radix2fft_algorithm1_6_2(x, w)
n = length(x);  t = log2(n);
x = fi_bitreverse(x,n);
for q=1:t
  L = 2^q; r = n/L; L2 = L/2;
  for k=0:(r-1)
    for j=0:(L2-1)
      temp          = w(L2-1+j+1) * x(k*L+j+L2+1);
      x(k*L+j+L2+1) = x(k*L+j+1)  - temp;
      x(k*L+j+1)    = x(k*L+j+1)  + temp;
    end
  end
end

Now we run a known signal through our algorithm to verify that our algorithm is correct.

y = fi_m_radix2fft_algorithm1_6_2(x0, w0);

y0 = fft(x0); % MATLAB's builtin FFT for comparison

fi_fft_demo_plot(x0,y,y0,Fs,'Double data', {'FFT Algorithm 1.6.2','Buil
tin FFT'});

The error is within tolerance of MATLAB's builtin FFT function.

The Algorithm in Fixed-Point M

Now let's convert the data to fixed-point and see if it still looks good. In this first pass, we use all defaults for fixed-point.

x = fi(x0);
w = fi(w0);

% Re-run the same algorithm with fixed-point inputs
y  = fi_m_radix2fft_algorithm1_6_2(x,w);
fi_fft_demo_plot(x,y,y0,Fs,'Fixed-point data', ...
                 {'Fixed-point FFT Algorithm 1.6.2','Builtin'});

Something has gone terribly wrong. The fixed-point FFT is very far from the built-in FFT. The error is much larger than what roundoff error would suggest, so we suspect that overflow has occurred.

Scaling to Prevent Overflow

By using an overall scaling of 1/n, we can avoid conditional scaling and still prevent overflow from occurring for any input.

If we scale only the input to the first stage of a length-n FFT by 1/n, we obtain a noise-to-signal ratio proportional to n^2 [Oppenheim & Schafer 1989, equation 9.101], [Welch 1969].

However, by scaling the input to each of the log2(n) stages by 1/2, we can obtain an overall scaling of 1/n and produce a noise-to-signal ratio proportional to n [Oppenheim & Schafer 1989, equation 9.105], [Welch 1969].

Hence, our algorithm becomes:

type_nocomments_nosubfunctions fi_m_radix2fft_withscaling

function xc = fi_m_radix2fft_withscaling(x, w)
n = length(x);  t = log2(n);
x = fi_bitreverse(x,n);
xc = complex(x,0);
for q=1:t
    L = 2^q; r = n/L; L2 = L/2;
    for k=0:(r-1)
        for j=0:(L2-1)
            temp           = w(L2-1+j+1) * xc(k*L+j+L2+1);
            xc(k*L+j+L2+1) = bitsra(xc(k*L+j+1) - temp, 1);
            xc(k*L+j+1)    = bitsra(xc(k*L+j+1) + temp, 1);
        end
    end
end

Run the scaled algorithm with fixed-point data.

x = fi(x0);
w = fi(w0);

y = fi_m_radix2fft_withscaling(x,w);

fi_fft_demo_plot(x, y, y0/n, Fs, 'Fixed-point data', ...
                 {'Fixed-point FFT with scaling','Scaled builtin'});

Conversion to C

Once we have a working fixed-point algorithm with no overflows, we work on converting it to C with the following constraints:

type_nocomments fi_c_radix2fft_withscaling.h

void fi_c_radix2fft_withscaling(int16_T* xr, int16_T* xi, int16_T* wr, int16
_T* wi,
                                int n, int nw, int t,
                                int Wfraclen)
{
    int32_T tempr, tempi;
    int q, i, j, k;
    int n1, n2, n3;
    int L, kL, r, L2;
    bitreverse(xr,xi,n);
    for (q=1; q<=t; q++) {
        L = 1; L <<= q;
        r = 1; r <<= (t-q);
        L2 = L>>1;
        kL = 0;
        for (k=0; k<r; k++) {
            for (j=0; j<L2; j++) {
                n3     = kL + j;
                n2     = n3 + L2;
                n1     = L2 - 1 + j;
                tempr  = (int32_T)wr[n1]*(int32_T)xr[n2] - (int32_T)wi[n1]*(
int32_T)xi[n2];
                tempi  = (int32_T)wr[n1]*(int32_T)xi[n2] + (int32_T)wi[n1]*(
int32_T)xr[n2];
                xr[n2] = ((((int32_T)xr[n3])<<Wfraclen) - tempr)>&g
t;(Wfraclen+1);
                xi[n2] = ((((int32_T)xi[n3])<<Wfraclen) - tempi)>&g
t;(Wfraclen+1);
                xr[n3] = ((((int32_T)xr[n3])<<Wfraclen) + tempr)>&g
t;(Wfraclen+1);
                xi[n3] = ((((int32_T)xi[n3])<<Wfraclen) + tempi)>&g
t;(Wfraclen+1);
            }
            kL += L;
        }
    }
}

There are no built-in complex data types in C, so we work directly with the real and imaginary parts of all variables.

When we compute the product of the twiddle factors and data, we store the result in a 32-bit integer. When integers are multiplied, they keep their least-significant bits; so to align the product with the data before addition, we must left-shift the data by the fraction length of the twiddle factors, then right-shift by the same amount before overwriting the data.

We have taken care to simplify the computation of indices in the inner loop.

Setting Up Fixed-Point M to Behave Like Integer C

To set up our fixed-point M code to behave like our C code, we use the fimath object with the following settings. We keep the least-significant bits of products and sums. We compute products and sums in 32 bits. Although the overflow characteristics of signed integers in C is not specified by the ANSII standard, it typically wraps on overflow. Since we are right-shifting and letting the bits fall off the end in our C code, we set the rounding mode in M to be floor rounding:

F = fimath;
F.ProductMode       = 'KeepLSB';
F.ProductWordLength = 32;
F.SumMode           = 'KeepLSB';
F.SumWordLength     = 32;
F.OverflowMode      = 'wrap';
F.RoundMode         = 'floor';
F.CastBeforeSum     = true;

Also, we specify the word length for the data. We let the scaling be automatically set for the data x. We specify the twiddle factors w to use fractional scaling since we know that they are complex sinusoids:

wordlength = 16;
x = fi(x0, true, wordlength);
w = fi(w0, true, wordlength, wordlength-1);

Set the fixed-point math attributes for the data and twiddle factors.

x.fimath = F;
w.fimath = F;

y = fi_m_radix2fft_withscaling(x,w);
fi_fft_demo_plot(x, y, y0/n, Fs, 'Fixed-point data', ...
                 {'Fixed-point FFT with C attributes','Scaled builtin'}
);

Validate the C Program

Now, let's compile and test that our algorithm in C is bit-faithful to our fixed-point algorithm in M.

The source files for the C MATLAB Executable (C-MEX) are:

 fi_c_radix2fft_withscaling.c
 fi_c_radix2fft_withscaling.h

To compile the C-MEX file, type

 mex fi_c_radix2fft_withscaling.c

at the MATLAB command line.

The input to the C program expects 16-bit signed integers for the data, so we convert the input to int16. The C algorithm also needs to know the fraction length of the twiddle factors:

yint = fi_c_radix2fft_withscaling(int16(x), int16(w), w.fractionlength);

fi_fft_demo_plot(int16(x),yint,int16(y),Fs,'Integer data', ...
                 {'C algorithm','M algorithm'});

Tips for Debugging the C Program

In the above plots, we see that the output of our C program exactly matches the output of our fixed-point M program. But if the results are different, then here are some suggestions for debugging the code:

      p = fipref;
      p.NumberDisplay      = 'int';
      p.NumericTypeDisplay = 'short';
      p.FimathDisplay      = 'none';

      x = ones(1,n)

     or

      x = [1 zeros(1,n-1)]

Alternate Variations

An optimization that is often made is to skip multiplies by the twiddle factor W^0 = exp(2*pi*0) = 1, which is the initial multiply in the innermost loop. An example of skipping the multiply by W^0 can be found in these files.

 fi_m_radix2fft_skip_w0.m
 fi_c_radix2fft_skip_w0.c

Scaling after every stage is sufficient to guarantee that no overflow will occur for any input; however it is not always necessary. Another method of scaling is ''block floating-point'' in which scaling is only done when an overflow is detected. An example of block floating-point can be found in these files.

 fi_m_radix2fft_blockfloatingpoint.m
 fi_c_radix2fft_blockfloatingpoint.c

close all

Compiling the Fixed-Point M-File into a C-MEX Function Using Embedded MATLAB™ MEX

The fixed-point M-File can also be compiled into a C-MEX function using the Embedded MATLAB™ C-MEX generation function (emlmex). The resulting fixed-point C-MEX function typically runs faster than the M-File, especially for large inputs.

The emlmex function compiles the M-File fi_m_radix2fft_withscaling into a C-MEX function. This step requires the creation of a temporary directory and write permissions in this directory.

emlmexdir = [tempdir filesep 'emlmexdir'];
if ~exist(emlmexdir,'dir')
    mkdir(emlmexdir);
end
emlcurdir = pwd;
cd(emlmexdir);

Compile the M-File into a MEX File

emlmex fi_m_radix2fft_withscaling -eg {x,w} -o fi_m_radix2fft_withscali
ng_mex;

Run the MEX Function with the Fixed-Point Inputs

y1 = fi_m_radix2fft_withscaling_mex(x,w);
fi_fft_demo_plot(x, y, y1, Fs, 'Fixed-point data', ...
                 {'M Function','C-MEX Function'});

You can also generate embeddable C code from the M-File using the Embedded MATLAB™ code generation utility (emlc).

if license('test','Simulink') && license('test','Real-Time_Workshop'
)
    emlc fi_m_radix2fft_withscaling -eg {x,w} -o fi_m_radix2fft_withsca
ling_mex;
end

Clean up the Temporary Directory

cd(emlcurdir);
clear fi_m_radix2fft_withscaling_mex;
status = rmdir(emlmexdir,'s');
close all;

Fixed-Point FFT in Signal Processing Blockset™

A fixed-point FFT is available in the Signal Processing Blockset product, so we demonstrate it here for comparison. Note that fixed-point objects can be passed in and out of Simulink® via the Signal From Workspace and Signal To Workspace blocks. The algorithm in the blockset optimizes away multiplication by 1 and orders the computations to minimize the number of accesses to the twiddle-factor table.

if license('test','Signal_Blocks') && license('test','Fixed-Point_Bl
ocks')

% Run the simulation
simopts = simset('SrcWorkspace','current');
sim('fi_mdl_radix2fft_withscaling',[],simopts)

% Open the model
fi_mdl_radix2fft_withscaling

end

Compare Simulink results with double:

if license('test','Signal_Blocks') && license('test','Fixed-Point_Bl
ocks')

fi_fft_demo_plot(x, y_sim, y0/n ,Fs, 'data',{'Simulink','Double'});
close_system('fi_mdl_radix2fft_withscaling');

end

Real-Time Workshop® Code

To see the Generated Code Report that was produced by Real-Time Workshop® for this model, click on this link: fi_mdl_radix2fft_withscaling_codegen_rpt.html

References

Charles Van Loan, Computational Frameworks for the Fast Fourier Transform, SIAM, 1992, ISBN 0-89871-285-8, http://www.mathworks.com/support/books/book1384.html.

Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time Signal Processing, Prentice Hall, 1989, ISBN 0-13-216292-X.

Peter D. Welch, A Fixed-Point Fast Fourier Transform Error Analysis, IEEE Transactions on Audio and Electroacoustics, Vol. AU-17, No. 2, June 1969, pp. 151-157.