FixedPoint Designer 

This example shows how to convert a textbook version of the Fast Fourier Transform (FFT) algorithm into fixedpoint MATLAB® code.
Run the following code to copy functions from the FixedPoint 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)
On this page… 

Use Min/Max Instrumentation to Identify Overflows 
FFT is a complexvalued 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:(n1))/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); % Timedomain signal x0 = complex(x0); % The textbook algorithm requires % the input to be complex y = fft(x0); % Frequencydomain 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 timedomain signal at those frequencies.
Note the reflected peaks at 3.5 and 3.8 Hz. When the input to an FFT is realvalued, as it is in this case, then the output y is conjugatesymmetric:
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 decimationintime unitstride 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 pseudocode, the algorithm in the textbook is as follows.
Algorithm 1.6.2. If is a complex vector of length and , then the following algorithm overwrites with .
The textbook algorithm uses zerobased indexing. is an nbyn Fouriertransform matrix, is an nbyn bitreversal permutation matrix, and is a complex vector of twiddle factors. The twiddle factors, , 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(n1,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')
To implement the algorithm in MATLAB, you can use the fidemo.fi_bitreverse function to bitreverse the input sequence, and you must add one to the indices to convert them from zerobased to onebased.
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:(r1) for j=0:(L21) temp = w(L21+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 builtin FFT for comparison fidemo.fi_fft_demo_plot(real(x0),y,y0,Fs,'Double data', {'FFT Algorithm 1.6.2','Builtin FFT'});
Because the error is within tolerance of the MATLAB builtin FFT function, you know you have correctly implemented the algorithm.
Now, try converting the data to fixedpoint and see if the algorithm still looks good. In this first pass, you use all the defaults for signed fixedpoint data by using the sfi constructor.
x = sfi(x0); % Convert to signed fixedpoint w = sfi(w0); % Convert to signed fixedpoint % Rerun the same algorithm with the fixedpoint inputs y = fi_m_radix2fft_algorithm1_6_2(x,w); fidemo.fi_fft_demo_plot(real(x),y,y0,Fs,'Fixedpoint data', ... {'Fixedpoint FFT Algorithm 1.6.2','Builtin'});
Note that the magnitude plot (center) of the fixedpoint FFT does not resemble the plot of the builtin 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 fiobject 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 regenerated, 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 FixedPoint 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 lengthn FFT by 1/n, you obtain a noisetosignal 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 noisetosignal ratio proportional to n [Oppenheim & Schafer 1989, equation 9.105], [Welch 1969].
An efficient way to scale by 1/2 in fixedpoint is to rightshift 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:(r1) for j=0:(L21) temp = w(L21+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 fixedpoint 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, 'Fixedpoint data', ... {'Fixedpoint FFT with scaling','Scaled builtin'});
You can see that the scaled fixedpoint FFT algorithm now matches the builtin FFT to a tolerance that is expected for 16bit fixedpoint data.
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, DiscreteTime Signal Processing, Prentice Hall, 1989.
Peter D. Welch, "A FixedPoint Fast Fourier Transform Error Analysis," IEEE® Transactions on Audio and Electroacoustics, Vol. AU17, No. 2, June 1969, pp. 151157.
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;