Convert Fast Fourier Transform (FFT) to Fixed Point

This example shows how to convert a textbook version of the Fast Fourier Transform (FFT) algorithm into fixed-point MATLAB® code.

Run the following code to copy functions from the Fixed-Point Designer™ examples directory into a temporary directory so this example doesn't interfere with your own work.

tempdirObj = fidemo.fiTempdir('fi_radix2fft_demo');
copyfile(fullfile(matlabroot,'toolbox','fixedpoint','fidemos','+fidemo',...
                  'fi_m_radix2fft_algorithm1_6_2.m'),'.','f');
copyfile(fullfile(matlabroot,'toolbox','fixedpoint','fidemos',...
                  'fi_m_radix2fft_withscaling.m'),'.','f');

Run the following code to capture current states, and reset the global states.

FIPREF_STATE = get(fipref);
reset(fipref)

Textbook FFT Algorithm

FFT is a complex-valued linear transformation from the time domain to the frequency domain. For example, if you construct a vector as the sum of two sinusoids and transform it with the FFT, you 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
x0 = complex(x0);                           % The textbook algorithm requires
                                            % the input to be complex
y  = fft(x0);                               % Frequency-domain transformation

figure(gcf); clf
subplot(211); plot(t,real(x0),'b.-'); xlabel('Time (s)'); ylabel('Amplitude');legend('x0')
subplot(212); plot(f,abs(y),'m.-'); xlabel('Frequency (Hz)'); ylabel('Magnitude');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, then the output y is conjugate-symmetric:

$$y(k) = \mbox{conj}(y(N-k)).$$

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

This example uses the decimation-in-time unit-stride FFT shown in Algorithm 1.6.2 on page 45 of the book Computational Frameworks for the Fast Fourier Transform by Charles Van Loan.

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

Algorithm 1.6.2. If $x$ is a complex vector of length $n$ and $n = 2^t$, then the following algorithm overwrites $x$ with $F_nx$.

$$\begin{array}{llll}
   \multicolumn{4}{l}{x = P_nx}\\
   \multicolumn{4}{l}{w = w_n^{(long)}\mbox{\hspace*{3em}(See Van Loan \S 1.4.11.)}}\\
   \mbox{for}\ q\ & \multicolumn{3}{l}{ = 1:t}\\
       & \multicolumn{3}{l}{L=2^q;\ r=n/L;\ L_\ast=L/2;}\\
       & \mbox{for}\ k\ & \multicolumn{2}{l}{=0:r-1}\\
       & & \mbox{for}\ j\ & =0:L_\ast-1\\
       & &                & \tau  = w(L_\ast-1+j) \cdot x(kL+j+L_\ast)\\
       & &                & x(kL+j+L_\ast) = x(kL+j)  - \tau\\
       & &                & x(kL+j)    = x(kL+j)  + \tau\\
       & & \mbox{end}\\
       & \mbox{end}\\
  \mbox{end}\\
\end{array}$$

The textbook algorithm uses zero-based indexing. $F_n$ is an n-by-n Fourier-transform matrix, $P_n$ 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:

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;
% Equation 1.4.11, p. 34
while L<=n
  theta = 2*pi/L;
  % Algorithm 1.4.1, p. 23
  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 = fidemo.fi_radix2twiddles(n);
polar(angle(w0),abs(w0),'o')
title('Twiddle Factors: Complex roots of unity')

Verify Floating-Point Code

To implement the algorithm in MATLAB, you can use the fidemo.fi_bitreverse function to bit-reverse the input sequence, and you must add one to the indices to convert them from zero-based to one-based.

function x = fi_m_radix2fft_algorithm1_6_2(x, w)
n = length(x);  t = log2(n);
x = fidemo.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

To verify that you correctly implemented the algorithm in MATLAB, run a known signal through it and compare the results to the results produced by the MATLAB FFT function.

y = fi_m_radix2fft_algorithm1_6_2(x0, w0);

y0 = fft(x0); % MATLAB's built-in FFT for comparison

fidemo.fi_fft_demo_plot(real(x0),y,y0,Fs,'Double data', {'FFT Algorithm 1.6.2','Built-in FFT'});

Because the error is within tolerance of the MATLAB built-in FFT function, you know you have correctly implemented the algorithm.

Identify Fixed-Point Issues

Now, try converting the data to fixed-point and see if the algorithm still looks good. In this first pass, you use all the defaults for signed fixed-point data by using the sfi constructor.

x = sfi(x0);  % Convert to signed fixed-point
w = sfi(w0);  % Convert to signed fixed-point

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

Note that the magnitude plot (center) of the fixed-point FFT does not resemble the plot of the built-in FFT. The error (bottom plot) is much larger than what you would expect to see for round off error, so it is likely that overflow has occurred.

Use Min/Max Instrumentation to Identify Overflows

To instrument the MATLAB® code, you create a MEX function from the MATLAB® function using the buildInstrumentedMexbuildInstrumentedMex command. The inputs to buildInstrumentedMex are the same as the inputs to fiaccelfiaccel, but buildInstrumentedMex has no fi-object restrictions. The output of buildInstrumentedMex is a MEX function with instrumentation inserted, so when the MEX function is run, the simulated minimum and maximum values are recorded for all named variables and intermediate values.

The '-o' option is used to name the MEX function that is generated. If the '-o' option is not used, then the MEX function is the name of the MATLAB® function with '_mex' appended. You can also name the MEX function the same as the MATLAB® function, but you need to remember that MEX functions take precedence over MATLAB® functions and so changes to the MATLAB® function will not run until either the MEX function is re-generated, or the MEX function is deleted and cleared.

Create the input with a scaled double datatype so its values will attain full range and you can identify potential overflows.

x_scaled_double = fi(x0,'DataType','ScaledDouble');
buildInstrumentedMex fi_m_radix2fft_algorithm1_6_2 ...
    -o fft_instrumented -args {x_scaled_double w}

Run the instrumented MEX function to record min/max values.

y_scaled_double = fft_instrumented(x_scaled_double,w);

Show the instrumentation results.

showInstrumentationResults fft_instrumented

You can see from the instrumentation results that there were overflows when assigning into the variable x.

Modify the Algorithm to Address Fixed-Point Issues

The magnitude of an individual bin in the FFT grows, at most, by a factor of n, where n is the length of the FFT. Hence, by scaling your data by 1/n, you can prevent overflow from occurring for any input.

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

However, if you scale the input to each of the stages of the FFT by 1/2, you 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].

An efficient way to scale by 1/2 in fixed-point is to right-shift the data. To do this, you use the bit shift right arithmetic function bitsra. After scaling each stage of the FFT, and optimizing the index variable computation, your algorithm becomes:

function x = fi_m_radix2fft_withscaling(x, w)
n = length(x);  t = log2(n);
x = fidemo.fi_bitreverse(x,n);
% Generate index variables as integer constants so they are not computed in
% the loop.
LL = int32(2.^(1:t)); rr = int32(n./LL); LL2 = int32(LL./2);
for q=1:t
    L = LL(q); r = rr(q); L2 = LL2(q);
    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) = bitsra(x(k*L+j+1) - temp, 1);
            x(k*L+j+1)    = bitsra(x(k*L+j+1) + temp, 1);
        end
    end
end

Run the scaled algorithm with fixed-point data.

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

y = fi_m_radix2fft_withscaling(x,w);
fidemo.fi_fft_demo_plot(real(x), y, y0/n, Fs, 'Fixed-point data', ...
                        {'Fixed-point FFT with scaling','Scaled built-in'});

You can see that the scaled fixed-point FFT algorithm now matches the built-in FFT to a tolerance that is expected for 16-bit fixed-point data.

References

Charles Van Loan, Computational Frameworks for the Fast Fourier Transform, SIAM, 1992.

Cleve Moler, Numerical Computing with MATLAB, SIAM, 2004, Chapter 8 Fourier Analysis.

Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time Signal Processing, Prentice Hall, 1989.

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.

Run the following code to restore the global states.

fipref(FIPREF_STATE);
clearInstrumentationResults fft_instrumented
clear fft_instrumented

Run the following code to delete the temporary directory.

tempdirObj.cleanUp;
Was this topic helpful?