Develop Fixed-Point Algorithms

Open Script

This example shows how to develop and verify a simple fixed-point algorithm. This example follows these steps for algorithm development:

1) Implement a second-order filter algorithm and simulate in double-precision floating-point.

2) Instrument the code to visualize the dynamic range of the output and state.

3) Convert the algorithm to fixed point by changing the data type of the variables. The algorithm itself does not change.

4) Compare and plot the fixed-point and floating-point results.

Define Double-Precision Floating-Point Variables

Develop the algorithm in double-precision floating-point. The algorithm used in this example is a second-order lowpass filter that removes the high frequencies in the input signal.

Define numerator coefficients, a, and denominator coefficients, b.

b = [ 0.25 0.5      0.25    ];
a = [ 1    0.09375  0.28125 ];

Generate a random input that has both high and low frequencies.

s = rng; rng(0,'v5uniform');
x = randn(1000,1);
rng(s); % restore |rng| state

Pre-allocate the output, y, and state, z, for speed.

y = zeros(size(x));
z = [0;0];

Implement Data Type Independent Algorithm

This is a second-order filter that implements the standard difference equation:

|y(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) - a(2)*y(n-1) - a(3)*y(n-2)|

for k=1:length(x)
    y(k) =  b(1)*x(k) + z(1);
    z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k);
    z(2) =  b(3)*x(k)         - a(3)*y(k);
end

Save the floating-point result.

ydouble = y;

Instrument Floating-Point Code to Visualize Dynamic Range

To convert to fixed point, we need to know the range of the variables. Depending on the complexity of an algorithm, this task can be simple or quite challenging. In this example, the range of the input value is known, so selecting an appropriate fixed-point data type is simple. Concentrate on the output (y) and states (z) since their range is unknown. To view the dynamic range of the output and states, modify the code slightly to instrument it.

Create a NumericTypeScope object and view the dynamic range of states (z).

z = [0;0]; % Reset states

hscope1 = NumericTypeScope;

for k=1:length(x)
    y(k) =  b(1)*x(k) + z(1);
    z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k);
    z(2) =  b(3)*x(k)         - a(3)*y(k);

    % Process the data and update the visual
    step(hscope1,z);
end

% clear the information stored in the object hscope1
reset(hscope1);

% Create a |NumericTypeScope| object and view the dynamic range of the
% output (|y|).
hscope2 = NumericTypeScope;
step(hscope2,y);

First, analyze the information displayed for variable z (state). From the histogram, observe that the dynamic range lies within (2^(1) 2^(-12)].

By default, the scope uses a word length of 16 bits with zero tolerable overflows. This results in a data type of numerictype(true,16,14) since at least 2 integer bits are needed to avoid overflows. You can get more information on the statistical data from the Input Data and Resulting Type panels. The Input Data panel shows that the data has both positive and negative values and hence a signed quantity which is reflected in the suggested numerictype. Also, the maximum data value is 1.51 which can be represented by the suggested type.

Next, look at variable y (output). From the histogram plot we see that the dynamic range lies within (2^(1) 2^(-13)].

Define Fixed-Point Variables

Convert variables to use fixed-point data types and run the algorithm again.

Enable logging to see the overflows and underflows introduced by the selected data types.

FIPREF_STATE = get(fipref);
reset(fipref)
fp = fipref;
default_loggingmode = fp.LoggingMode;
fp.LoggingMode = 'On';

Capture the present state of and reset the global fimath to the factory settings.

globalFimathAtStart = fimath;
resetglobalfimath;

Define the fixed-point types for the variables in the below format: fi(Data, Signed, WordLength, FractionLength)

b = fi(b, 1, 8, 6);
a = fi(a, 1, 8, 6);

x = fi(x, 1, 16, 13);
y = fi(zeros(size(x)), 1, 16, 13);
z = fi([0;0], 1, 16, 14);

Implement the Same Data Type Independent Algorithm

for k=1:length(x)
    y(k) =  b(1)*x(k) + z(1);
    z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k);
    z(2) =  b(3)*x(k)         - a(3)*y(k);
end
% Reset the logging mode.
fp.LoggingMode = default_loggingmode;

In this example, we have redefined the fixed-point variables with the same names as the floating-point so that we could inline the algorithm code for clarity. However, it is a better practice to enclose the algorithm code in a MATLAB® file function that could be called with either floating-point or fixed-point variables.

Compare and Plot the Floating-Point and Fixed-Point Results

Plot the magnitude response of the floating-point and fixed-point results and the response of the filter to see if the filter behaves as expected when it is converted to fixed-point.

n = length(x);
f = linspace(0,0.5,n/2);
x_response = 20*log10(abs(fft(double(x))));
ydouble_response = 20*log10(abs(fft(ydouble)));
y_response = 20*log10(abs(fft(double(y))));
plot(f,x_response(1:n/2),'c-',...
    f,ydouble_response(1:n/2),'bo-',...
    f,y_response(1:n/2),'gs-');
ylabel('Magnitude in dB');
xlabel('Normalized Frequency');
legend('Input','Floating point output','Fixed point output','Location','Best');
title('Magnitude response of Floating-point and Fixed-point results');

h = fft(double(b),n)./fft(double(a),n);
h = h(1:end/2);
clf
hax = axes;
plot(hax,f,20*log10(abs(h)));
set(hax,'YLim',[-40 0]);
title('Magnitude response of the filter');
ylabel('Magnitude in dB')
xlabel('Frequency');

Notice that the high frequencies in the input signal are attenuated by the low-pass filter which is the expected behavior.

Plot the Error

clf
n = (0:length(y)-1)';
e = double(lsb(y));
plot(n,double(y)-ydouble,'.-r', ...
     [n(1) n(end)],[e/2 e/2],'c', ...
     [n(1) n(end)],[-e/2 -e/2],'c')
text(n(end),e/2,'+1/2 LSB','HorizontalAlignment','right','VerticalAlignment','bottom')
text(n(end),-e/2,'-1/2 LSB','HorizontalAlignment','right','VerticalAlignment','top')
xlabel('n (samples)'); ylabel('error')

Implement the Algorithm in Simulink

If you have Simulink® and Fixed-Point Designer™, you can run this model, which is the equivalent of the algorithm above. The output, y_sim is a fixed-point variable equal to the variable y calculated above in MATLAB code.

As in the MATLAB code, the fixed-point parameters in the blocks can be modified to match an actual system; these have been set to match the MATLAB code in the example above. Double-click the blocks to see the settings.

% Set up the From Workspace variable
x_sim.time = n;
x_sim.signals.values = x;
x_sim.signals.dimensions = 1;

% Run the simulation
out_sim = sim('fitdf2filter_demo', 'SaveOutput', 'on', ...
    'SrcWorkspace', 'current');

% Open the model
open_system('fitdf2filter_demo')

% Verify that the Simulink results are the same as the MATLAB file
isequal(y, out_sim.get('y_sim'))

ans =

  logical

   1

Assumptions Made for This Example

In order to simplify the example, we have taken the default math parameters: round-to-nearest, saturate on overflow, full precision products and sums. You can modify all of these parameters to match an actual system.

The settings were chosen as a starting point in algorithm development. Save a copy of this MATLAB file, start playing with the parameters, and see what effects they have on the output. How does the algorithm behave with a different input? See the help for fi, fimath, and numerictype for information on how to set other parameters, such as rounding mode, and overflow mode.

% Reset the global fimath
globalfimath(globalFimathAtStart);
fipref(FIPREF_STATE);