# Loren on the Art of MATLAB

## Recent Question about Speed with Subarray Calculations

Recently someone asked me to explain the speed behavior doing a calculation using a loop and array indexing vs. getting the subarray first.

### Contents

#### Example

Suppose I have a function of two inputs, the first input being the column (of a square array), the second, a scalar, and the output, a vector.

myfun = @(x,z) x'*x+z;


And even though this may be calculated in a fully vectorized manner, let's explore what happens when we work on subarrays from the array input.

I am now creating the input array x and the results output arrays for doing the calculation two ways, with an additional intermediate step in one of the methods.

n = 500;
x = randn(n,n);
result1 = zeros(n,1);
result2 = zeros(n,1);


#### First Method

Here we see and time the first method. In this one, we create a temporary array for x(:,k) n times through the outer loop.

tic
for k = 1:n
for z = 1:n
result1(z) = myfun(x(:,k), z);
end
result1 = result1+x(:,k);
end
runtime(1) = toc;


#### Second Method

In this method, we extract the column of interest first in the outer loop, and reuse that temporary array each time through the inner loop. Again we see and time the results.

tic
for k = 1:n
xt = x(:,k);
for z = 1:n
result2(z) = myfun(xt, z);
end
result2 = result2+x(:,k);
end
runtime(2) = toc;


#### Same Results?

First, let's make sure we get the same answer both ways. You can see that we do.

theSame = isequal(result1,result2)

theSame =
1


#### Compare Runtime

Next, let's compare the times. I want to remind you that doing timing from a script generally has more overhead than when the same code is run inside a function. We just want to see the relative behavior so we should get some insight from this exercise.

disp(['Run times are: ',num2str(runtime)])

Run times are: 2.3936      1.9558


#### What's Happening?

Here's what's going on. In the first method, we create a temporary variable n times through the outer loop, even though that array is a constant for a fixed column. In the second method, we extract the relevant column once, and reuse it n times through the inner loop.

Be thoughtful if you do play around with this. Depending on the details of your function, if the calculations you do each time are large compared to the time to extract a column vector, you may not see much difference between the two methods. However, if the calculations are sufficiently short in duration, then the repeated creation of the temporary variable could add a tremendous amount of overhead to the calculation. In general, you should not be worse off always capturing the temporary array the fewest number of times possible.

Have you noticed similar timing "puzzles" when analyzing one of your algorithms? I'd love to hear more here.

Get the MATLAB code

Published with MATLAB® R2013a

## Using Symbolic Math Toolbox to Compute Area Moments of Inertia

Once more I am pleased to introduce guest blogger Kai Gehrs. Kai has been a Software Engineer at MathWorks for the past five years working on the Symbolic Math Toolbox. He has a background in mathematics and computer science. He already contributed to my BLOG in the past writing about Using Symbolic Equations And Symbolic Functions In MATLAB as well as on approaches for Simplifying Symbolic Results.

### Contents

If you are interested in using MATLAB and the Symbolic Math Toolbox in teaching some basics in mechanical engineering, this might be of interest to you.

Computing area moments of inertia is an important task in mechanics. For example, area moments of inertia play a critical role in stress, dynamic, and stability analysis of structures. In this article, we use capabilities of the Symbolic Math Toolbox to compute area moments for cross sections of elliptical tubes.

We start with a basic case involving only numeric parameters, and then make the computations more general by introducing symbolic parameters.

All plots used in this article have been created in the MuPAD Notebook app. You can find the source code for these plots at the end of the article.

#### Basic Example: Cross Section of an Elliptical Tube

The following picture shows the cross section of an elliptical tube.

The following ellipses define outer and inner contours of the section:

• $y = \pm 2 \, \sqrt{1 - \frac{x^2}{9}}$ describe the outer contour line for $x \in [-3,3]$.
• $y = \pm \sqrt{1 - \frac{x^2}{4}}$ describe the inner contour line for $x \in [-2,2]$.

We will compute the area moment of inertia of this section with respect to the $x$-axis.

#### Area Moment Of Inertia

The moment of inertia of an area $A$ with respect to the $x$-axis is defined in terms of the double integral

$$I_x = {\int \int}_A y^2\, \mathrm{d}A.$$

You can find this definition on Wikipedia.

#### Math Behind This Example

In our example, the area $A$ is the hatched area shown in the previous plot. Taking advantage of the symmetry of the section with respect to both the $x$- and $y$-axes, we can restrict our computations to the first quadrant of the $x-y$-plane.

To compute the necessary double integral over the hatched area, we divide this area into two separate areas $A1$ and $A2$.

To compute the final area moment of inertia about the x-axis, we multiply the sum of the double integrals over $A_1$ and $A_2$ by four.

$$I_x = {\int\int}_A y^2 \, \mathrm dA = 4 \cdot \Big({\int\int}_{A_1} y^2 \, \mathrm dA_1 + {\int\int}_{A_2} y^2 \, \mathrm dA_2\Big).$$

Now, we only need to compute these two double integrals. The mathematical theory behind multi-dimensional integration tells us that we can rewrite each of these double integrals in terms of two integrals:

$$I_1 = {\int\int}_{A_1} y^2 \, \mathrm dA_1 = \int_0^2 \int_{\sqrt{1 - \frac{x^2}{4}}}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2 \, \mathrm dy \, \mathrm dx,$$

$$I_2 = {\int\int}_{A_2} y^2 \, \mathrm dA_2 = \int_2^3 \int_{0}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2\, \mathrm dy\, \mathrm dx.$$

The int function from the Symbolic Math Toolbox can compute these integrals.

#### Symbolic Integration

We define $x$ and $y$ as symbolic variables:

syms x y;


We start computing the inner integral $\int_{\sqrt{1 - \frac{x^2}{4}}}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2 \, \mathrm dy$ of $I_1$:

innerI1 = int(y^2, y, sqrt(1-x^2/4), 2*sqrt(1-x^2/9))

innerI1 =
(8*(9 - x^2)^(3/2))/81 - (4 - x^2)^(3/2)/24


Next, we integrate innerI1 with respect to $x$ from $0$ to $2$. This provides $I_1$:

I1 = int(innerI1, x, 0, 2)

I1 =
3*asin(2/3) - pi/8 + (74*5^(1/2))/81


We use the same strategy to compute $I_2$. First, we compute the inner integral $\int_{0}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2\, \mathrm dy$. Then, we integrate the resulting expression with respect to $x$ from $2$ to $3$:

I2 = int(int(y^2, y, 0, 2*sqrt(1-x^2/9)), x, 2, 3);


Hence, the area moment of inertia of the elliptical tube with respect to the $x$-axis is

Ix = 4 * (I1 + I2)

Ix =
(11*pi)/2


We can approximate the result numerically by using the vpa function. For example, we approximate the result using 5 significant (nonzero) digits:

vpa(Ix, 5)

ans =
17.279


#### Advanced Example: Cross Section of an Elliptical Tube Defined by Symbolic Parameters

We can use a numerical integrator, such as MATLAB's integral2, to compute the area moment of inertia in the previous example. A numerical integrator might return slightly less accurate results, but other than that there is not much benefit from using symbolic integration there.

But what if we consider the cross section of a more general elliptical tube whose shape is defined by arbitrary symbolic parameters?

For example, we want to compute the area moment of inertia of this elliptical tube.

Its contour lines are still ellipses, but they are defined by the symbolic parameters $r_1$, $r_2$, $R_1$ and $R_2$:

• $y = \pm R_2 \, \sqrt{1 - \frac{x^2}{R_1^2}}$ describe the outer contour line for $x \in [-R_1,R_1]$.
• $y = \pm r_2 \, \sqrt{1 - \frac{x^2}{r_1^2}}$ describe the inner contour line for $x \in [-r_1,r_1]$.

Now we need to compute these integrals:

$$I_1 = {\int\int}_{A_1} y^2 \, \mathrm dA_1 = \int_0^{r_1} \int_{r_2 \, \sqrt{1 - \frac{x^2}{r_1^2}}}^{R_2 \, \sqrt{1 - \frac{x^2}{R_1^2}}} y^2 \, \mathrm dy \, \mathrm dx,$$

$$I_2 = {\int\int}_{A_2} y^2 \, \mathrm dA_2 = \int_{r_1}^{R_1} \int_{0}^{R_2 \, \sqrt{1 - \frac{x^2}{R_1^2}}} y^2\, \mathrm dy\, \mathrm dx.$$

Although the situation is more complicated now, we can still use the same strategy as above, using symbolic variables instead of numbers. Nevertheless, we need to add one more step: specify relationships between symbolic parameters. This is because the variables $r_1$, $r_2$, $R_1$, and $R_2$ can be any complex number, unless we explicitly restrict their values. For example, the system does not know if these variables are positive or negative, if $r_1 < R_1$ or otherwise. In this example, we want to specify that $r_1>0$, $r_2>0$, $R_1>r_1$, and $R_2>r_2$. In the Symbolic Math Toolbox, such restrictions are called assumptions on variables. They can be set by using the assume and assumeAlso functions:

syms x y r1 r2 R1 R2;
assume(r1 > 0);
assume(r2 > 0);
assumeAlso(r1 < R1);
assumeAlso(r2 < R2);
I1 = int(int(y^2, y, r2*sqrt(1-x^2/r1^2), R2*sqrt(1-x^2/R1^2)), x, 0, r1);
I2 = int(int(y^2, y, 0, R2*sqrt(1-x^2/R1^2)), x, r1, R1);
Ixgeneral = 4 * (I1 + I2);
pretty(Ixgeneral)


/          4     / r1 \ \         /                     4     / r1 \ \
|      3 R1  asin| -- | |         |             4   3 R1  asin| -- | |
3 |                \ R1 / |       3 |      3 pi R1              \ R1 / |
4 R2  | #1 + ---------------- |   4 R2  | #1 - -------- + ---------------- |
\             8         /         \         16             8         /
------------------------------- - ------------------------------------------
3                                      3
3 R1                                   3 R1

3
pi r1 r2
- ---------
4

where

/     2        3 \
| 5 R1  r1   r1  |    2     2 1/2
#1 == | -------- - --- | (R1  - r1 )
\    8        4  /


Now, let's substitute our symbolic parameters with the same values that we had in the first example. Do we get the same result? Let's double-check: when substituting $r_1=2,r_2,=1,R_1=3,R_2=2$ in Ixgeneral, do we get the same result for Ix as obtained for the section initially considered?

As expected, the general solution Ixgeneral reduces to the specific solution Ix when we substitute our original dimensions:

vpa(subs(Ixgeneral, [r1, R1, r2, R2],[2, 3, 1, 2]), 5)

ans =
17.279


Note that it is essential to set the assumptions on $r_1$, $r_2$, $R_1$ and $R_2$. By default, the Symbolic Math Toolbox assumes that all variables including symbolic parameters are arbitrary complex numbers. This makes computing integrals much more complicated and time-consuming. For example, try computing Ixgeneral without setting the assumptions on $r_1$, $r_2$, $R_1$ and $R_2$. In general, it can be impossible to even get the result without any additional assumptions.

To generate the graphics shown in this article, use MuPAD Notebook app. Note that the code in this section does not run in the MATLAB Command Window.

To open the MuPAD Notebook app:

1. On the MATLAB Toolstrip, click the Apps tab.
2. On the Apps tab, click the down arrow at the end of the Apps section.
3. Under Math, Statistics and Optimization, click the MuPAD Notebook button.

Alternatively, type mupad in the MATLAB Command Window.

Paste the following commands to a MuPAD notebook to obtain the four graphics used above:

Ellipse:= (x,a,b) -> sqrt(b^2*(1-x^2/a^2)):
f1 := plot::Function2d(Ellipse(x,2,1), x = -3..3):
f2 := plot::Function2d(Ellipse(x,3,2), x = -3..3):
f3 := plot::Function2d(Ellipse(x,3,2), x = -3..-2):
f4 := plot::Function2d(Ellipse(x,3,2), x = 2..3):
f5 := plot::Function2d(-Ellipse(x,2,1), x = -3..3):
f6 := plot::Function2d(-Ellipse(x,3,2), x = -3..3):
f7 := plot::Function2d(-Ellipse(x,3,2), x = -3..-2):
f8 := plot::Function2d(-Ellipse(x,3,2), x = 2..3):
f9 := plot::Function2d(Ellipse(x,2,1), x = 0..2, VerticalAsymptotesVisible = FALSE):
f10 := plot::Function2d(Ellipse(x,3,2), x = 0..3, VerticalAsymptotesVisible = FALSE):
f11 := plot::Function2d(Ellipse(x,3,2), x = 2..3, VerticalAsymptotesVisible = FALSE):
plot(plot::Hatch(f1,f2),plot::Hatch(f3),plot::Hatch(f4),f1,f2,
plot::Hatch(f5,f6),plot::Hatch(f7),plot::Hatch(f8),f5,f6,
Scaling = Constrained,
GridVisible = TRUE,
VerticalAsymptotesVisible = FALSE,
Height = 120, Width = 200,
Header = "Cross section of elliptical tube (x-y-plane)"):
plot(plot::Hatch(f9,f10),plot::Hatch(f11),f9,f10,
Scaling = Constrained,
GridVisible = TRUE,
Height = 120, Width = 200,
Header = "Restriction to first quadrant (x-y-plane)"):
plot(plot::Hatch(f9,f10,Color=RGB::Magenta),plot::Hatch(f11,Color = RGB::Green),f9,f10,
Scaling=Constrained,
GridVisible = TRUE,
plot::Text2d("A1",[1.0,1.25]),
plot::Text2d("A2",[2.5,0.25]),
Height = 120, Width = 200,
Header = "Integration areas"):
plot(plot::Hatch(f1,f2),plot::Hatch(f3),plot::Hatch(f4),f1,f2,
plot::Hatch(f5,f6),plot::Hatch(f7),plot::Hatch(f8),f5,f6,
plot::Text2d("r1",[1.8,0.05]),
plot::Text2d("r2",[0.05,0.85]),
plot::Text2d("R1",[2.74,0.05]),
plot::Text2d("R2",[0.05,1.84]),
Scaling = Constrained,
GridVisible = TRUE,
VerticalAsymptotesVisible = FALSE,
XTicksLabelsVisible = FALSE,
YTicksLabelsVisible = FALSE,
ViewingBox = [-3.5..3.5,-2.5..2.5],
Height = 120, Width = 200,
Header = "Cross section of elliptical tube (more general situation)"):

#### Have You Tried Symbolic Math Toolbox?

Have you used the Symbolic Math Toolbox in computations related to mechanics? Let me know here.

Get the MATLAB code

Published with MATLAB® R2013a

## MATLAB to FPGA using HDL Coder(TM)

It's my pleasure to introduce guest blogger Kiran Kintali. Kiran is the product development lead for HDL Coder at MathWorks. In this post, Kiran introduces a new capability in HDL Coder™ that generates synthesizable VHDL/Verilog code directly from MATLAB and highlights some of the key features of this new MATLAB based workflow.

### Contents

#### Introduction to HDL Code Generation from MATLAB

If you are using MATLAB to model digital signal processing (DSP) or video and image processing algorithms that eventually end up in FPGAs or ASICs, read on...

FPGAs provide a good compromise between general purpose processors (GPPs) and application specific integrated circuits (ASICs). GPPs are fully programmable but are less efficient in terms of power and performance; ASICs implement dedicated functionality and show the best power and performance characteristics, but require extremely expensive design validation and implementation cycles. FPGAs are also used for prototyping in ASIC workflows for hardware verification and early software development.

Due to the order of magnitude performance improvement when running high-throughput, high-performance applications, algorithm designers are increasingly using FPGAs to prototype and validate their innovations instead of using traditional processors. However, many of the algorithms are implemented in MATLAB due to the simple-to-use programming model and rich analysis and visualization capabilities. When targeting FPGAs or ASICs these MATLAB algorithms have to be manually translated to HDL.

For many algorithm developers who are well-versed with software programming paradigms, mastering the FPGA design workflow is a challenge. Unlike software algorithm development, hardware development requires them to think parallel. Other obstacles include: learning the VHDL or Verilog language, mastering IDEs from FPGA vendors, and understanding esoteric terms like "multi-cycle path" and "delay balancing".

In this post, I describe an easier path from MATLAB to FPGAs. I will show how you can automatically generate HDL code from your MATLAB algorithm, implement the HDL code on an FPGA, and use MATLAB to verify your HDL code.

#### MATLAB to Hardware Workflow

The process of translating MATLAB designs to hardware consists of the following steps:

1. Model your algorithm in MATLAB - use MATLAB to simulate, debug, and iteratively test and optimize the design.
2. Generate HDL code - automatically create HDL code for FPGA prototyping.
3. Verify HDL code - reuse your MATLAB test bench to verify the generated HDL code.
4. Create and verify FPGA prototype - implement and verify your design on FPGAs.

There are some unique challenges in translating MATLAB to hardware. MATLAB code is procedural and can be highly abstract; it can use floating-point data and has no notion of time. Complex loops can be inferred from matrix operations and toolbox functions.

Implementing MATLAB code in hardware involves:

• Converting floating-point MATLAB code to fixed-point MATLAB code with optimized bit widths suitable for efficient hardware generation.
• Identifying and mapping procedural constructs to concurrent area- and speed-optimized hardware operations.
• Introducing the concept of time by adding clocks and clock rates to schedule the operations in hardware.
• Creating resource-shared architectures to implement expensive operators like multipliers and for-loop bodies.
• Mapping large persistent arrays to block RAM in hardware

HDL Coder™ simplifies the above tasks though workflow automation.

#### Example MATLAB Algorithm

Let’s take a MATLAB function implementing histogram equalization and go through this workflow. This algorithm, implemented in MATLAB, enhances image contrast by transforming the values in an intensity image so that the histogram of the output image is approximately flat.

type mlhdlc_heq.m

% Histogram Equalization Algorithm
function [pixel_out] = mlhdlc_heq(x_in, y_in, pixel_in, width, height)

persistent histogram
persistent transferFunc
persistent histInd
persistent cumSum

if isempty(histogram)
histogram = zeros(1, 2^8);
transferFunc = zeros(1, 2^8);
histInd = 0;
cumSum = 0;
end

% Figure out indices based on where we are in the frame
if y_in < height && x_in < width % valid pixel data
histInd = pixel_in + 1;
elseif y_in == height && x_in == 0 % first column of height+1
histInd = 1;
elseif y_in >= height % vertical blanking period
histInd = min(histInd + 1, 2^8);
elseif y_in < height % horizontal blanking - do nothing
histInd = 1;
end

%Read histogram

%Read transfer function

%If valid part of frame add one to pixel bin and keep transfer func val
if y_in < height && x_in < width
transValWrite = transValRead; %Write back same value
cumSum = 0;
elseif y_in >= height %In blanking time index through all bins and reset to zero
histValWrite = 0;
cumSum = transValWrite;
else
end

%Write histogram
histogram(histInd) = histValWrite;

%Write transfer function
transferFunc(histInd) = transValWrite;

pixel_out = transValRead;


#### Example MATLAB Test Bench

Here is the test bench that verifies that the algorithm works with an example image. (Note that this testbench uses Image Processing Toolbox functions for reading the original image and plotting the transformed image after equalization.)

type mlhdlc_heq_tb.m

%% Test bench for Histogram Equalization Algorithm
clear mlhdlc_heq;
testFile = 'office.png';

% Get intensity part of color image
YCBCR = rgb2ycbcr(RGB);
imgOrig = YCBCR(:,:,1);

[height, width] = size(imgOrig);
imgOut = zeros(height,width);
hBlank = 20;
% make sure we have enough vertical blanking to filter the histogram
vBlank = ceil(2^14/(width+hBlank));

for frame = 1:2
disp(['working on frame: ', num2str(frame)]);
for y_in = 0:height+vBlank-1
%disp(['frame: ', num2str(frame), ' of 2, row: ', num2str(y_in)]);
for x_in = 0:width+hBlank-1
if x_in < width && y_in < height
pixel_in = double(imgOrig(y_in+1, x_in+1));
else
pixel_in = 0;
end

             [pixel_out] = mlhdlc_heq(x_in, y_in, pixel_in, width, height);
             if x_in < width && y_in < height
imgOut(y_in+1,x_in+1) = pixel_out;
end
end
end
end
% Make color image from equalized intensity image
% Rescale image
imgOut = double(imgOut);
imgOut(:) = imgOut/max(imgOut(:));
imgOut = uint8(imgOut*255);

YCBCR(:,:,1) = imgOut;
RGBOut = ycbcr2rgb(YCBCR);

figure(1)
subplot(2,2,1); imshow(RGB, []);
title('Original Image');
subplot(2,2,2); imshow(RGBOut, []);
title('Equalized Image');
subplot(2,2,3); hist(double(imgOrig(:)),2^14-1);
title('Histogram of original Image');
subplot(2,2,4); hist(double(imgOut(:)),2^14-1);
title('Histogram of equalized Image');


Let's simulate this algorithm to see the results.

mlhdlc_heq_tb

working on frame: 1
working on frame: 2


The HDL Workflow Advisor (see the snapshot below) helps automate the steps and provides a guided path from MATLAB to hardware. You can see the following key steps of the workflow in the left pane of the workflow advisor:

1. Fixed-Point Conversion
2. HDL Code Generation
3. HDL Verification
4. HDL Synthesis and Analysis

Let's look at each workflow step in detail.

Fixed-Point Conversion

Signal processing applications are typically implemented using floating-point operations in MATLAB. However, for power, cost, and performance reasons, these algorithms need to be converted to use fixed-point operations when targeting hardware. Fixed-point conversion can be very challenging and time-consuming, typically demanding 25 to 50 percent of the total design and implementation time. The automatic floating-point to fixed-point conversion workflow in HDL Coder™ can greatly simplify and accelerate this conversion process.

The floating-point to fixed-point conversion workflow consists of the following steps:

1. Verify that the floating-point design is compatible with code generation.
2. Propose fixed-point types based on computed ranges, either through the simulation of the testbench or through static analysis that propagates design ranges to compute derived ranges for all the variables.
3. Generate fixed-point MATLAB code by applying proposed fixed-point types.
4. Verify the generated fixed-point code and compare the numerical accuracy of the generated fixed-point code with the original floating point code.

Note that this step is optional. You can skip this step if your MATLAB design is already implemented in fixed-point.

HDL Code Generation

The HDL Code Generation step generates HDL code from the fixed-point MATLAB code. You can generate either VHDL or Verilog code that implements your MATLAB design. In addition to generating synthesizable HDL code, HDL Coder™ also generates various reports, including a traceability report that helps you navigate between your MATLAB code and the generated HDL code, and a resource utilization report that shows you, at the algorithm level, approximately what hardware resources are needed to implement the design, in terms of adders, multipliers, and RAMs.

During code generation, you can specify various optimization options to explore the design space without having to modify your algorithm. In the Design Space Exploration and Optimization Options section below, you can see how you can modify code generation options and optimize your design for speed or area.

HDL Verification

Standalone HDL test bench generation:

HDL Coder™ generates VHDL and Verilog test benches from your MATLAB scripts for rapid verification of generated HDL code. You can customize an HDL test bench using a variety of options that apply stimuli to the HDL code. You can also generate script files to automate the process of compiling and simulating your code in HDL simulators. These steps help to ensure the results of MATLAB simulation match the results of HDL simulation.

HDL Coder™ also works with HDL Verifier to automatically generate two types of cosimulation testbenches:

• HDL cosimulation-based verification works with Mentor Graphics® ModelSim® and QuestaSim®, where MATLAB and HDL simulation happen in lockstep.
• FPGA-in-the-Loop simulation allows you to run a MATLAB simulation with an FPGA board in strict synchronization. You can use MATLAB to feed real world data into your design on the FPGA, and ensure that the algorithm will behave as expected when implemented in hardware.

HDL Synthesis

Apart from the language-related challenges, programming for FPGAs requires the use of complex EDA tools. Generating a bitstream from the HDL design and programming the FPGA can be daunting tasks. HDL Coder™ provides automation here, by creating project files for Xilinx® and Altera® that are configured with the generated HDL code. You can use the workflow steps to synthesize the HDL code within the MATLAB environment, see the results of synthesis, and iterate on the MATLAB design to improve synthesis results.

#### Design Space Exploration and Optimization Options

HDL Coder™ provides the following optimizations to help you explore the design space trade-offs between area and speed. You can use these options to explore various architectures and trade-offs without having to manually rewrite your algorithm.

Speed Optimizations

• Pipelining : To improve the design’s clock frequency, HDL Coder enables you to insert pipeline registers in various locations within your design. For example, you can insert registers at the design inputs and outputs, and also at the output of a given MATLAB variable in your algorithm.
• Distributed Pipelining : HDL Coder also provides an optimization based on retiming to automatically move pipeline registers you have inserted to maximize clock frequency, by minimizing the delay through combinational paths in your design.

Area Optimizations

• RAM mapping: HDL Coder™ maps matrices to wires or registers in hardware. If persistent matrix variables are mapped to registers, they can take up a large amount of FPGA area. HDL Coder™ automatically maps persistent matrices to block RAM to improve area efficiency. The challenge in mapping MATLAB matrices to block RAM is that block RAM in hardware typically has a limited set of read and write ports. HDL Coder™ solves this problem by automatically partitioning and scheduling the matrix reads and writes to honor the block RAM’s port constraints, while still honoring the other control- and data-dependencies in the design.
• Resource sharing: This optimization identifies functionally equivalent multiplier operations in MATLAB code and shares them. You can control the amount of multiplier sharing in the design.
• Loop streaming: A MATLAB for-loop creates a FOR_GENERATE loop in VHDL. The body of the loop is replicated as many times in hardware as the number of loop iterations. This results in an inefficient use of area. The loop streaming optimization creates a single hardware instance of the loop body that is time-multiplexed across loop iterations.
• Constant multiplier optimization: This design level optimization converts constant multipliers into shift and add operations using canonical signed digit (CSD) techniques.

#### Best Practices

Now, let's look at few best practices related to writing MATLAB code when targeting FPGAs.

When writing a MATLAB design:

• Use the code generation subset of MATLAB supported for HDL code generation.
• Keep the top-level interface as simple as possible. The top-level function size, types, and complexity determine the interface of the chip implemented in hardware.
• Do not pass in a big chunk of parallel data into the design. Parallel data requires a large number of IO pins on the chip, and would probably not be synthesizable. In a typical image processing design, you should serialize the pixels as inputs and buffer them internally in the algorithm.

When writing a MATLAB test bench:

• Call the design from the testbench function.
• Exercise the design thoroughly. This is particularly important for floating-point to fixed-point conversion, where HDL Coder™ determines the ranges of the variables in the algorithm based on the values the testbench assigns to the variables. You can reuse this testbench to generate an HDL testbench for testing the generated hardware.
• Simulate the design with the testbench prior to code generation to make sure there are no simulation errors, and to make sure all the required files are on the path.

#### Conclusion

HDL Coder™ provides a seamless workflow when you want to implement your algorithm in an FPGA. In this post, I have shown you how to take an image processing algorithm written in MATLAB, convert it to fixed-point, generate HDL code, verify the generated HDL code using the test bench, and finally, synthesize the design and implement it in hardware.

We hope this brief introduction to the HDL Coder™ and MATLAB-to-HDL code generation, verification framework has shown how you can quickly get started on implementing your MATLAB designs and target FPGAs. Please let us know in the comments for this post how you might use this new functionality. Or, if you've already tried using HDL Coder™, let us know about your experiences here.

Get the MATLAB code

Published with MATLAB® R2013a

## New Datatype under Development for Possible MATLAB Release

There is a new datatype we are playing around with that we hope to make available in an upcoming release and we would like your input beforehand.

### Contents

#### New Datatype in Action

Let me show you the new datatype in action so you can first get a feel for it.

inputData = magic(3)

inputData =
8     1     6
3     5     7
4     9     2

outputValues = dis(inputData);


#### Let's Examine the Output

outputValues

Why are you asking?
4     2     3
1     9     6
8     7     5


Well, that's a bit strange, isn't it? I wonder what the relationship between inputData and outputValues is. What can we learn about outputValues?

whos outputValues

  Name              Size            Bytes  Class    Attributes

outputValues      1x1               248  dis



Well, it's a dis array. Let's look at it again.

outputValues

What's it matter to you?
5     6     9
4     3     2
7     8     1


Say what? Let's check it a few more times.

outputValues
outputValues
outputValues

Who wants to know?
6     3     9
5     8     7
4     2     1
Who wants to know?
5     2     9
1     4     8
7     6     3
Who are you to ask me that?
5     3     7
9     8     2
1     4     6


Hoping you get the double meaning here - the dis array not only mixes up the values of the input for display purposes, but also tries to gently *dis*respect you along the way.

Even though this is a silly class, I'll show you the code so you can see how simple it is to make such a class.

type dis

classdef dis
%dis dis is a class.
%   In fact, it's a declasse class.

properties
Data
end
properties (Access=protected)
Original
end

properties (Constant)
'What''s it matter to you?',...
'Who are you to ask me that?',...
'Who wants to know?',...
'What''s the big deal?'}
end

methods
function display(obj)
obj.Data(:) = obj.Data(randperm(numel(obj.Data)));
disp(obj.Data)
end
function obj = dis(in)
obj.Original = in;
obj.Data = reshape(in(randperm(numel(in))),size(in));
end
end

end



#### Should We Invest More Resources?

Of course, I could also add some numeric functions like plus to dis, but I didn't take the time, in case you didn't find this possible new MATLAB addition useful. So please share your thoughts with us here.

Get the MATLAB code

Published with MATLAB® R2013a

## Multiple Y Axes

We were musing here about how common it is to want more than two Y axes on a plot. Checking out the File Exchange, there seem to be several candidates, indicating that this is something at least some people find useful.

### Contents

#### Sample Plot

Here's a sample plot using plotyy that comes with MATLAB.

x = 0:0.01:20;
y1 = 200*exp(-0.05*x).*sin(x);
y2 = 0.8*exp(-0.5*x).*sin(10*x);
plotyy(x,y1,x,y2,'plot');


#### List of Some Possibilities

In addition to plotyy in MATLAB, here's a list of some of the candidates from the File Exchange.

#### What are You Plotting with More Y Axes?

I am curious to know what kind of data or results you are plotting so that having multiple y-axes makes a compelling presentation. Let us know here.

Get the MATLAB code

Published with MATLAB® R2013a

## Using the MATLAB Unit Testing Infrastructure for Grading Assignments

Steven Lord, Andy Campbell, and David Hruska are members of the Quality Engineering group at MathWorks who are guest blogging today to introduce a new feature in R2013a, the MATLAB unit testing infrastructure. There are several submissions on the MATLAB Central File Exchange related to unit testing of MATLAB code. Blogger Steve Eddins wrote one highly rated example back in 2009. In release R2013a, MathWorks included in MATLAB itself a MATLAB implementation of the industry-standard xUnit testing framework.

If you're not a software developer, you may be wondering if this feature will be of any use to you. In this post, we will describe one way someone who may not consider themselves a software developer may be able to take advantage of this framework using the example of a professor grading students' homework submissions. That's not to say that the developers in the audience should move on to the next post; you can use these tools to test your own code just like a professor can use them to test code written by his or her students.

There is a great deal of functionality in this feature that we will not show here. For more information we refer you to the MATLAB Unit Testing Framework documentation.

### Contents

#### Background

In order to use this feature, you should be aware of how to define simple MATLAB classes in classdef files, how to define a class that inherits from another, and how to specify attributes for methods and properties of those classes. The object-oriented programming documentation describes these capabilities.

#### Problem Statement

As a professor in an introductory programming class, you want your students to write a program to compute Fibonacci numbers. The exact problem statement you give the students is:

Create a function "fib" that accepts a nonnegative integer n and returns
the nth Fibonacci number. The Fibonacci numbers are generated by this
relationship:
F(0) = 1
F(1) = 1
F(n) = F(n-1) + F(n-2) for integer n > 1
Your function should throw an error if n is not a nonnegative integer.

#### Basic Unit Test

The most basic MATLAB unit test is a MATLAB classdef class file that inherits from the matlab.unittest.TestCase class. Throughout the rest of this post we will add additional pieces to this basic framework to increase the capability of this test and will change its name to reflect its increased functionality.

dbtype basicTest.m

1     classdef basicTest < matlab.unittest.TestCase
2
3     end

test = basicTest

test =
basicTest with no properties.


#### Running a Test

To run the test, we can simply pass test to the run function. There are more advanced ways that make it easier to run a group of tests, but for our purposes (checking one student's answer at a time) this will be sufficient. When you move to checking multiple students' answers at a time, you can use run inside a for loop.

Since basicTest doesn't actually validate the output from the student's function, it doesn't take very long to execute.

results = run(test)

results =
0x0 TestResult array with properties:

Name
Passed
Failed
Incomplete
Duration
Totals:
0 Passed, 0 Failed, 0 Incomplete.
0 seconds testing time.


Let's say that a student named Thomas submitted a function fib.m as his solution to this assignment. Thomas's code is stored in a sub-folder named thomas. To set up our test to check Thomas's answer, we add the folder holding his code to the path.

addpath('thomas');
dbtype fib.m

1     function y = fib(n)
2     if n <= 1
3         y = 1;
4     else
5         y = fib(n-1)+fib(n-2);
6     end



#### Test that F(0) Equals 1

The basicTest is a valid test class, and we can run it, but it doesn't actually perform any validation of the student's test file. The methods that will perform that validation need to be written in a methods block that has the attribute Test specified.

The matlab.unittest.TestCase class includes qualification methods that you can use to test various qualities of the results returned by the student files. The qualification method that you will likely use most frequently is the verifyEqual method, which passes if the two values you pass into it are equal and reports a test failure if they are not.

The documentation for the matlab.unittest.TestCase class lists many other qualification methods that you can use to perform other types of validation, including testing the data type and size of the results; matching a string result to an expected string; testing that a given section of code throws a specific errors or issues a specific warning; and many more.

This simple test builds upon generalTest by adding a test method that checks that the student's function returns the value 1 when called with the input 0.

dbtype simpleTest.m

1     classdef simpleTest < matlab.unittest.TestCase
2         methods(Test)
3             function fibonacciOfZeroShouldBeOne(testCase)
4                 % Evaluate the student's function for n = 0
5                 result = fib(0);
6                 testCase.verifyEqual(result, 1);
7             end
8         end
9     end


Thomas's solution to the assignment satisfies this basic check. We can use the results returned from run to display the percentage of the tests that pass.

results = run(simpleTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running simpleTest
.
Done simpleTest
__________

results =
TestResult with properties:

Name: 'simpleTest/fibonacciOfZeroShouldBeOne'
Passed: 1
Failed: 0
Incomplete: 0
Duration: 0.0112
Totals:
1 Passed, 0 Failed, 0 Incomplete.
0.011168 seconds testing time.
100% Passed.


#### Test that F(pi) Throws an Error

Now that we have a basic positive test in place we can add in a test that checks the behavior of the student's function when passed a non-integer value (like n = pi) as input. The assignment stated that when called with a non-integer value, the student's function should error. Since the assignment doesn't require a specific error to be thrown, the test passes as long as fib(pi) throws any exception.

dbtype errorCaseTest.m

1     classdef errorCaseTest < matlab.unittest.TestCase
2         methods(Test)
3             function fibonacciOfZeroShouldBeOne(testCase)
4                 % Evaluate the student's function for n = 0
5                 result = fib(0);
6                 testCase.verifyEqual(result, 1);
7             end
8             function fibonacciOfNonintegerShouldError(testCase)
9                 testCase.verifyError(@()fib(pi), ?MException);
10            end
11        end
12    end


Thomas forgot to include a check for a non-integer valued input in his function, so our test should indicate that by reporting a failure.

results = run(errorCaseTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running errorCaseTest
.
================================================================================
Verification failed in errorCaseTest/fibonacciOfNonintegerShouldError.

---------------------
Framework Diagnostic:
---------------------
verifyError failed.
--> The function did not throw any exception.

Expected Exception Type:
MException

Evaluated Function:
@()fib(pi)

------------------
Stack Information:
------------------
In C:\Program Files\MATLAB\R2013a\toolbox\matlab\testframework\+matlab\+unittest\+qualifications\Verifiable.m (Verifiable.verifyError) at 637
In H:\Documents\LOREN\MyJob\Art of MATLAB\errorCaseTest.m (errorCaseTest.fibonacciOfNonintegerShouldError) at 9
================================================================================
.
Done errorCaseTest
__________

Failure Summary:

Name                                            Failed  Incomplete  Reason(s)
=============================================================================================
errorCaseTest/fibonacciOfNonintegerShouldError    X                 Failed by verification.

results =
1x2 TestResult array with properties:

Name
Passed
Failed
Incomplete
Duration
Totals:
1 Passed, 1 Failed, 0 Incomplete.
0.026224 seconds testing time.
50% Passed.


Another student, Benjamin, checked for a non-integer value in his code as you can see on line 2.

rmpath('thomas');
dbtype fib.m

1     function y = fib(n)
2     if (n ~= round(n)) || n < 0
3         error('N is not an integer!');
4     elseif n == 0 || n == 1
5         y = 1;
6     else
7         y = fib(n-1)+fib(n-2);
8     end


Benjamin's code passed both the test implemented in the fibonacciOfZeroShouldBeOne method (which we copied into errorCaseTest from simpleTest) and the new test case implemented in the fibonacciOfNonintegerShouldError method.

results = run(errorCaseTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running errorCaseTest
..
Done errorCaseTest
__________

results =
1x2 TestResult array with properties:

Name
Passed
Failed
Incomplete
Duration
Totals:
2 Passed, 0 Failed, 0 Incomplete.
0.010132 seconds testing time.
100% Passed.


#### Basic Test for Students, Advanced Tests for Instructor

The problem statement given earlier in this post is a plain text description of the homework assignment we assigned to the students. We can also state the problem for the students in code (if they're using release R2013a or later) by giving them a test file they can run just like simpleTest or errorCaseTest. They can directly use this "requirement test" to ensure their functions satisfy the requirements of the assignment.

dbtype studentTest.m

1     classdef studentTest < matlab.unittest.TestCase
2         methods(Test)
3             function fibonacciOfZeroShouldBeOne(testCase)
4                 % Evaluate the student's function for n = 0
5                 result = fib(0);
6                 testCase.verifyEqual(result, 1);
7             end
8             function fibonacciOfNonintegerShouldError(testCase)
9                 testCase.verifyError(@()fib(pi), ?MException);
10            end
11        end
12    end


In order for the student's code to pass the assignment, it will need to pass the test cases given in the studentTest unit test. However, we don't want to use studentTest as the only check of the student's code. If we did, the student could write their function to cover only the test cases in the student test file.

We could solve this problem by having two separate test files, one containing the student test cases and one containing additional test cases the instructor uses in the grading process. Can we avoid having to run both test files manually or duplicating the code from the student test cases in the instructor test? Yes!

To do so, we write an instructor test file to incorporate, through inheritance, the student test file. We can then add additional test cases to the instructor test file. When we run this test it should run three test cases; two inherited from studentTest, fibonacciOfZeroShouldBeOne and fibonacciOfNonintegerShouldError, and one from instructorTest itself, fibonacciOf5.

dbtype instructorTest.m

1     classdef instructorTest < studentTest
2         % Because the student test file is a matlab.unittest.TestCase and
3         % instructorTest inherits from it, instructorTest is also a
4         % matlab.unittest.TestCase.
5
6         methods(Test)
7             function fibonacciOf5(testCase)
8                 % Evaluate the student's function for n = 5
9                 result = fib(5);
10                testCase.verifyEqual(result, 8, 'Fibonacci(5) should be 8');
11            end
12        end
13    end


Let's look at Eric's test file that passes the studentTestFile test, but in which he completely forgot to implement the F(n) = F(n-1)+F(n-2) recursion step.

rmpath('benjamin');
dbtype fib.m

1     function y = fib(n)
2     if (n ~= round(n)) || n < 0
3         error('N is not an integer!');
4     end
5     y = 1;



It should pass the student unit test.

results = run(studentTest);
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running studentTest
..
Done studentTest
__________

100% Passed.


It does NOT pass the instructor unit test because it fails one of the test cases.

results = run(instructorTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

Running instructorTest
..
================================================================================
Verification failed in instructorTest/fibonacciOf5.

----------------
Test Diagnostic:
----------------
Fibonacci(5) should be 8

---------------------
Framework Diagnostic:
---------------------
verifyEqual failed.
--> NumericComparator failed.
--> The values are not equal using "isequaln".

Actual Value:
1
Expected Value:
8

------------------
Stack Information:
------------------
In C:\Program Files\MATLAB\R2013a\toolbox\matlab\testframework\+matlab\+unittest\+qualifications\Verifiable.m (Verifiable.verifyEqual) at 411
In H:\Documents\LOREN\MyJob\Art of MATLAB\instructorTest.m (instructorTest.fibonacciOf5) at 10
================================================================================
.
Done instructorTest
__________

Failure Summary:

Name                         Failed  Incomplete  Reason(s)
==========================================================================
instructorTest/fibonacciOf5    X                 Failed by verification.

results =
1x3 TestResult array with properties:

Name
Passed
Failed
Incomplete
Duration
Totals:
2 Passed, 1 Failed, 0 Incomplete.
0.028906 seconds testing time.
66.6667% Passed.


Benjamin, whose code we tested above, wrote a correct solution to the homework problem.

rmpath('eric');

results = run(instructorTest)
percentPassed = 100 * nnz([results.Passed]) / numel(results);
disp([num2str(percentPassed), '% Passed.']);

rmpath('benjamin');

Running instructorTest
...
Done instructorTest
__________

results =
1x3 TestResult array with properties:

Name
Passed
Failed
Incomplete
Duration
Totals:
3 Passed, 0 Failed, 0 Incomplete.
0.015946 seconds testing time.
100% Passed.


#### Conclusion

In this post, we showed you the basics of using the new MATLAB unit testing infrastructure using homework grading as a use case.

We checked that the student's code worked (by returning the correct answer) for one valid value and worked (by throwing an error) for one invalid value. We also showed how you can use this infrastructure to provide an aid/check for the students that you can also use as part of your grading.

We hope this brief introduction to the unit testing framework has shown you how you can make use of this feature even if you don't consider yourself a software developer. Let us know in the comments for this post how you might use this new functionality. Or, if you've already tried using matlab.unittest, let us know about your experiences here.

Get the MATLAB code

Published with MATLAB® R2013a

## Major New Release at Chebfun – Two Dimensional Capabilities

Perhaps you've read about Chebfun before here or on Cleve''s blog. Chebfun allows you to do numerical computing with functions and not simply numbers. If not, you might want to investigate. If so, you may be very interested in new capabilities available.

### Contents

#### Major New Capability

Just recently, the Chebfun team released a major new capability on their site. While Chebfun is extremely useful, it was restricted to solving problems in one dimension. Now you may solve problems with two independent variables using Chebfun2! As before, their website includes code you can freely download and many interesting examples. The examples are well documented and easy to learn from. I will soon have time to explore Chebfun2 because of an irrational amount of airline time.

I also love how easy it is to install - you can do it from MATLAB. They kindly post the code right on the website so you can paste it into your running session.

#### Many Examples

They have many examples included in the software, and posted on their web site, using the publish command in MATLAB. The high-level examples are organized into these categories:

1. Approximation
2. Complex Analysis
3. Geometry
4. Optimization
5. Statistics
6. Rootfinding
7. Vector Calculus
8. Fun !

I hope with this breadth of examples, you are able to find ones interesting or relevant (or both) to you. Once you've gotten a chance to try out Chebfun2, I'd love to hear about your experiences. Post them here to share with us.

Get the MATLAB code

Published with MATLAB® R2012b

## Logical Indexing – Multiple Conditions

I've talked about logical indexing before in some of the linked posts, but recent work makes me want to show it off again. One of the nice things about logical indexing is that it is very easy and natural to combine the results of different conditions to select items based on multiple criteria.

### Contents

#### What is Logical Indexing?

Suppose I have an array of integers, not sorted, and want to find the ones that are less than a certain number. Here's how I can do it using the function find.

X = randperm(20)
target = 5;

X =
Columns 1 through 13
3     7     1    16    20    15    13     9    14     8    11    19    18
Columns 14 through 20
4    10     5     6    17    12     2

ind = find(X < target)

ind =
1     3    14    20


You can see that find returns the indices into the array X that have values less than the target. And we can use these to extract the values.

Xtarget = X(ind)

Xtarget =
3     1     4     2


Another way to accomplish the same outcome is to use the logical expression to directly perform the indexing operation. Here's what I mean.

logInd = X < target

logInd =
Columns 1 through 13
1     0     1     0     0     0     0     0     0     0     0     0     0
Columns 14 through 20
1     0     0     0     0     0     1


MATLAB returns an array that matches the elements of the array X, element-by-element holding 1s where the matching values in X are the desired values, and 0s otherwise. The array logInd is not an array of double numbers, but have the class logical.

whos logInd

  Name        Size            Bytes  Class      Attributes



I can now use this array to extract the desired values from X.

XtargetLogical = X(logInd)

XtargetLogical =
3     1     4     2


Both methods return the results.

isequal(Xtarget, XtargetLogical)

ans =
1


#### Compound conditions.

Let me create an anonymous function that returns true (logical(1)) for values that are even integers.

iseven = @(x) ~logical(rem(x,2))

iseven =
@(x)~logical(rem(x,2))


Test iseven.

iseven(1:5)

ans =
0     1     0     1     0


#### Find Values Meeting More Than One Condition

Now I would like to find the values in X that are less than target and are even. This is very natural to do with logical indexing. We have the pieces of code we need already.

compoundCondInd = (X < target) & iseven(X)

compoundCondInd =
Columns 1 through 13
0     0     0     0     0     0     0     0     0     0     0     0     0
Columns 14 through 20
1     0     0     0     0     0     1


We can see we found suitable values at locations 3 and 19. And we can extract those values next.

X(compoundCondInd)

ans =
4     2


#### Did You Notice?

Did you see how easy it is to combine multiple conditions? I simply look for each of condition, getting back logical arrays, and then compute a logical array where the two input arrays are both true (via &). I could, of course, calculate a compound condition where only either one or the other condition needs to be true using logical or (via |).

#### A Recent Application

I recently used this in the context of finding suspect data values. I had 2 arrays, hourly temperature and speed. The problem is that when the temperature gets near or below freezing, the speed sensor might freeze. But I didn't want to delete ALL the values below freezing. So I looked for data where the temperature was sufficiently low AND the speed was very low (which could potentially mean the sensor was frozen). That way, I did not need to discard all data at low temperatures.

#### Have You Used Compound Indexing?

Did you do it like I did, using logical expressions? Or did you use some other techniques? What were you trying to achieve with your compound indexing? Let me know here.

Get the MATLAB code

Published with MATLAB® R2012b

## Introduction to Functional Programming with Anonymous Functions, Part 3

Tucker McClure is an Application Engineer with The MathWorks. He spends his time helping our customers accelerate their work with the right tools and problem-solving techniques. Today, he'll be discussing how "functional programming" can help create brief and powerful MATLAB code.

### Contents

#### Recap

When we left off, we had implemented conditional statements, recursion, and multi-line statements in anonymous functions, so today we'll tackle loops.

Before we get started, let's implement the functions that we'll need again.

iif     = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();
recur   = @(f, varargin) f(f, varargin{:});
curly   = @(x, varargin) x{varargin{:}};


#### Loops

Note that the recursive sequences we created in the last part could also have been implemented with for loops. For instance, here's factorial of n:

   factorial = 1;
for k = 1:n
factorial = k * factorial;
end

Many times, recursive functions can be written iteratively in loops. However, we can't use for or while in an anonymous function, so instead of asking how we can unwrap recursive functions into iterative loops, let's ask the reverse: how can we implement loops with recursive functions?

#### Loops via Recursion

To loop properly, one must know:

• What to do each iteration
• If the process should continue to the next iteration
• What's available when the loop begins

Allowing the "what to do" to be a function (fcn) of some state (x), the "if it should continue" to be another function (cont) of the state, and "what's available when the loop begins" to be the initial state (x0), we can write a loop function. This is a big step, so bear with me for some explanation!

On each step, the loop function will call the cont function, passing in all elements of the state, x, as in cont(x{:}). If that returns false (meaning we shouldn't continue), the current state, x, is returned. Otherwise, if we should continue, it calls fcn with all elements of the current state, as in fcn(x{:}), and passes the output from that to the next iteration. Letting this single iteration be denoted as f, we can build the anonymous function loop using our recur function.

loop = @(x0, cont, fcn) ...                                     % Header
recur(@(f, x) iif(~cont(x{:}), x, ...                    % Continue?
true,        @() f(f, fcn(x{:}))), ... %   Iterate
x0);                                               % from x0.


For this trivial example, the state is simply the iteration count. We'll increase the count every iteration until the count >= n and return the final count. All this does therefore is count from 0 to the input n. Not very interesting, but it demonstrates the loop.

count = @(n) loop({0}, ...          % Initialize state, k, to 0
@(k) k < n, ...   % While k < n
@(k) {k + 1});    %   k = k + 1 (returned as cell array)

arrayfun(count, 1:10)

ans =
[1]    [2]    [3]    [4]    [5]    [6]    [7]    [8]    [9]    [10]


I suppose that worked, but why are we using cell arrays to store the state, such as {0} and {k+1}? There are two reasons. First, if x is a cell array, then when we dump all elements of x into fcn, they become multiple arguments! That is, fcn(x{:}) is the same as fcn(x{1}, x{2}, ...). So instead of our function taking a big cell array for an input, it can take named arguments, which we'll use below. Second, we do this because it allows a function to return multiple elements that will be used by the next iteration, so if a function needed to return y and z, which would be arguments to the next iteration, it can simply return one cell array, {y, z}. It makes it easy to use. Here's a factorial example demonstrating this. The state is two different things: the iteration count, k, and factorial of the previous number, x. Note that both values of the state, k and x, are inputs to all of the functions. Note here how we're using @(k, x) for our functions. By allowing x to be a cell array, each element of the array becomes an argument such as k or x!

factorial = @(n) loop({1, 1}, ...          % Start with k = 1 and x = 1
@(k, x) k <= n, ...  % While k <= n
@(k, x) {k + 1, ...  %   k = k + 1;
k * x});    %   x = k * x;


Call it:

factorial(5)

ans =
[6]    [120]


Wait, we wanted 120 (the fifth number of the factorial sequence), so what's that 6 doing there?

#### A Better Loop

Remember how we return the full state? That's not very useful for this factorial example, as we get both k and the number we want as outputs in that cell array. Because the whole state isn't generally useful, let's add a cleanup function to our loop. We'll execute this when the loop is done (right after ~cont(...) returns false). Our cleanup function will take the full state and return only the important parts.

loop = @(x0, cont, fcn, cleanup) ...                            % Header
recur(@(f, x) iif(~cont(x{:}), cleanup(x{:}), ...        % Continue?
true,        @() f(f, fcn(x{:}))), ... %   Iterate
x0);                                               % from x0.


Now here's factorial, with clean output.

factorial = @(n) loop({1, 1}, ...         % Start with k = 1 and x = 1
@(k,x) k <= n, ...  % While k <= n
@(k,x) {k + 1, ...  %   k = k + 1;
k * x}, ... %   x = k * x;
@(k,x) x);          % End, returning x


The result:

factorial(5)

ans =
120


First seven numbers of factorial:

arrayfun(factorial, 1:7)

ans =
Columns 1 through 6
1           2           6          24         120         720
Column 7
5040


That's better.

I'll be the first to admit that the loop is a bit longer and much more rigid than a normal MATLAB loop. On the other hand, it can be used in anonymous functions, and its syntax has a certain cleanliness to it in that it doesn't modify any variables that live outside the loop; it has its own scope. This is one nice feature of loop being a function that takes code (functions) as arguments.

#### Doing More in a Loop

Let's say we want to do something else in the loop, but don't want its output passed to the next iteration, like printing something out. Remember the do_three_things example from last time? We executed numerous statements by putting them in a cell array and used curly to access the output we cared about. We can do that here, in a loop. For example, let's write out a function to print n digits of the factorial sequence. We'll use an array to store two things. The first will be the number that fprintf returns, which we don't care about. The second will be the cell array we want to return, k and x. We'll access that cell array with curly, as in curly({..., {k, x}}, 2), which just returns {k, x}.

say_it = @(k, x) fprintf('Factorial(%d): %d\n', k, x);
@(k,x) k <= n, ...              % While k <= n
{k + 1, ...       %   k = k + 1;
k * x}}, ...     %   x = k * x;
2), ...            %   Return {k+1,k*x}.
@(k,x) x);                      % End, returning x

print_factorial(7);

Factorial(1): 1
Factorial(2): 2
Factorial(3): 6
Factorial(4): 24
Factorial(5): 120
Factorial(6): 720
Factorial(7): 5040


Now we're executing multiple things and only returning what we want while inside a loop built built on recursion and anonymous conditionals! We've come a long way since Part 1.

As a practical note, recall that because these loops use recursion, there's a limit to the number of times they can loop (MATLAB has a recursion limit, which is a setting in Preferences). Also, a recursive implementation of a loop isn't the most efficient. For this reason, it's best to implement loop itself in a file that can then be used in the same way. If it's in a file, it can also be kept on the MATLAB path so that it can be used anywhere.

   function x = loop(x, cont, f, cleanup)
while cont(x{:})
x = f(x{:});
end
if nargin == 4
x = cleanup(x{:});
end
end

#### Final Example

This brings us to our final example. Below, we'll simulate a simple harmonic oscillator over time, using a structure to store dissimilar states, including a complete time history of the oscillator. This might simulate, for example, the sway of a lamp that's hanging from the ceiling after an earthquake.

% First, calculate a state transition matrix that represents a harmonic
% oscillator with damping. Multiplying this by |x| produces |x| at a
% slightly later time. The math here isn't important to the example.
Phi = expm(0.5*[0 1; -1 -0.2]);

% Now create the loop.
x   = loop({[1; 0], 1}, ...                  % Initial state, x = [1; 0]
@(x,k) k <= 100, ...              % While k <= 100
@(x,k) {[x, Phi * x(:, end)], ... %   Update x
k + 1}, ...               %   Update k
@(x,k) x);                        % End, return x

% Create a plot function.
plot_it = @(n, x, y, t) {subplot(2, 1, n), ...            % Select subplot.
plot(x(n, :)), ...               % Plot the data.
iif(nargin==4, @() title(t), ... % If there's a
true,      []), ...          % title, add it.
ylabel(y), ...                   % Label y
xlabel('Time (s)')};             % and x axes.

% Plot the result.
plot_it(1, x, 'Position (m)', 'Harmonic Oscillator');
plot_it(2, x, 'Velocity (m/s)');


#### Summary

That's it for loops via recursion!

Let's look back at what we did over these three parts. First, we started with a simple map utility function to demonstrate the function-of-functions idea. Then we created our ubiquitous inline if, which further enabled recursion (a conditional is necessary to make recursion stop!). We also showed using multiple statements by storing their outputs in a cell array. Finally, we created a loop construct on top of our recursion functions.

At this point, we've done more than just scratch the surface of functional programming. We've used MATLAB's interesting constructs, such as function handles, cell arrays, and varargin to implement a functional programming framework, allowing a new syntax within MATLAB, where code can be arguments to flow control functions. Here's a roundup of what we created.

map   = @(val, fcns) cellfun(@(f) f(val{:}), fcns);
mapc  = @(val, fcns) cellfun(@(f) f(val{:}), fcns, 'UniformOutput', 0);
iif   = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();
recur = @(f, varargin) f(f, varargin{:});
paren = @(x, varargin) x(varargin{:});
curly = @(x, varargin) x{varargin{:}};
loop  = @(x0,c,f,r)recur(@(g,x)iif(c(x{:}),@()g(g,f(x{:})),1,r(x{:})),x0);


These have also been programmed as "normal" MATLAB functions so that they can be kept on the path and used whenever they're needed. These can be found under "Functional Programming Constructs" in File Exchange, here.

Thanks for reading. I hope this has both enabled a new level of detail in anonymous functions in MATLAB and helped demonstrate the wide range of possibilities available within the MATLAB language.

Do you have other functional programming patterns you use in your code? For instance, a do-while loop is just like our loop above except that it always runs at least one iteration. Any ideas how to program this or other interesting constructs in anonymous functions? Please let us know here!

Get the MATLAB code

Published with MATLAB® R2012b

## Introduction to Functional Programming with Anonymous Functions, Part 2

Tucker McClure is an Application Engineer with The MathWorks. He spends his time helping our customers accelerate their work with the right tools and problem-solving techniques. Today, he'll be discussing how "functional programming" can help create brief and powerful MATLAB code.

### Contents

#### Recap

Last time, we said that functional programming was marked by storing functions as variables (function handles) and working with functions that act on other functions. We put these ideas together to implement our own version of a map function for handling multiple inputs and outputs from multiple functions simultaneously, and we created iif, an "inline if", to allow the use of conditional statements inside of anonymous functions. So how might we work with recursive functions -- functions of themselves? We'll see how a functional programming style allows us to implement recursive functionality inside anonymous functions, and this will pave the way for the final part, in which we'll implement loops, without ever using for or while (which we can't use in anonymous functions).

Before we get started, let's implement iif again; we're going to need it frequently.

iif = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();


#### Anonymous Function Recursion

Recall that a recursive function is a function that calls itself. It therefore needs some way to refer to itself. When we write an anonymous function, it isn't "named" (hence, "anonymous"), so it can't call itself by name. How can we get around this?

Let's start with a Fibonacci sequence example. Recall that the nth number of the Fibonacci sequence is the sum of the previous two numbers, starting with 1 and 1, yielding 1, 1, 2, 3, 5, 8, 13, 21, etc. This is easy to implement recursively.

   fib = @(n) iif(n <= 2, 1, ...                    % First two numbers
true,   @() fib(n-1) + fib(n-2)); % All later numbers

But hey, that can't work! We haven't defined fib yet, so how could this anonymous function call it? In fact, the anonymous function will never "know" we're referring to it as fib, so this won't work at all. Therefore, instead of trying to call fib directly, let's provide another input: the handle of a function to call, f.

fib = @(f, n) iif(n <= 2, 1, ...                      % First two numbers
true,   @() f(f, n-1) + f(f, n-2)); % All later numbers


Getting closer. Now, if we pass fib to fib along with the number we want, it will call fib, passing in fib as the first argument, recursively until we get our answer.

fib(fib, 6)

ans =
8


Ok, that's right. The sixth number of the sequence is 8. On the other hand, the syntax we've created is terrible. We have to provide the function to itself? I'd rather not. Instead, let's just write a new function that hands fib to fib along with the input n.

fib2 = @(n) fib(fib, n);

fib2(4)
fib2(5)
fib2(6)

ans =
3
ans =
5
ans =
8


That's a lot closer to what we want, but there's one more step. Let's write a function called recur to hand a function handle to itself, along with any other arguments. This makes recursion less cumbersome.

recur = @(f, varargin) f(f, varargin{:});


That was simple, so now let's re-write fib. The first argument to recur is the function, which we'll define inline. The second is n. That's all there is to it. It now reads as "Recursively call a function that, if k <= 2, returns one, and otherwise returns the recursive function of k-1 plus that of k-2, starting with the user's input n." (If it doesn't read quite this clearly at first, that's ok. It takes some getting used to. Comment liberally if necessary!)

fib = @(n) recur(@(f, k) iif(k <= 2, 1, ...
true,   @() f(f, k-1) + f(f, k-2)), ...
n);


And we can find the first ten numbers of the sequence via arrayfun.

arrayfun(fib, 1:10)

ans =
1     1     2     3     5     8    13    21    34    55


Factorial (f(n) = 1 * 2 * 3 * ... n) is another easy operation to represent recursively.

factorial = @(n) recur(@(f, k) iif(k == 0, 1, ...
true,   @() k * f(f, k-1)), n);
arrayfun(factorial, 1:7)

ans =
Columns 1 through 6
1           2           6          24         120         720
Column 7
5040


A number to an integer power has a nearly identical form. Here's 4.^(0:5).

pow = @(x, n) recur(@(f, k) iif(k == 0, 1, ...
true,   @() x * f(f, k-1)), n);
arrayfun(@(n) pow(4, n), 0:5)

ans =
1           4          16          64         256        1024


That was a big step for anonymous functions, using both recursion and an inline conditional together with ease. Like map and iif, recur, looks strange at first, but once it's been seen, it's hard to forget how it works (just make one of the inputs a function handle and pass it to itself). And recursion doesn't have to stop at interesting mathematical sequences of numbers. For instance, in the next part, we'll use this to implement loops in, but first, we'll need a some helper functions and a good way to execute multiple statements in an anonymous function.

#### Helpers

These little functions are useful in many circumstances, and we're going to need curly frequently.

paren = @(x, varargin) x(varargin{:});
curly = @(x, varargin) x{varargin{:}};


They allow us to write x(3, 4) as paren(x, 3, 4) and similarly for curly braces. That is, now we can think of parentheses and curly braces as functions! At first this might not seem useful. However, imagine writing a function to return the width and height of the screen. The data we need is available from this call:

get(0, 'ScreenSize')

ans =
1           1        1920        1200


However, we don't need those preceeding ones. We could save the output to a variable, say x, and then access x(3:4), but if we need this in an anonymous function, we can't save to a variable. How do we access just elements 3 and 4? There are numerous ways, but paren and curly are similar to constructs found in other languages and are easy to use, so we'll use those here.

Now we can write our screen_size function to return just the data we want.

screen_size = @() paren(get(0, 'ScreenSize'), 3:4);

screen_size()

ans =
1920        1200


While on the subject, note that we can actually use any number of indices or even ':'.

magic(3)
paren(magic(3), 1:2, 2:3)
paren(magic(3), 1:2, :)

ans =
8     1     6
3     5     7
4     9     2
ans =
1     6
5     7
ans =
8     1     6
3     5     7


We do the same with the curly braces. Here, the regular expression pattern will match both 'rain' and 'Spain', but we'll only select the second match.

spain = curly(regexp('The rain in Spain....', '\s(\S+ain)', 'tokens'), 2)

spain =
'Spain'


(Click for Regexp help.)

It also works with ':' (note that the single quotes are required).

[a, b] = curly({'the_letter_a', 'the_letter_b'}, ':')

a =
the_letter_a
b =
the_letter_b


#### Executing Multiple Statements

With curly in place, let's examine something a little different. Consider the following:

do_three_things = @() {fprintf('This is the first thing.\n'), ...
fprintf('This is the second thing.\n'), ...
max(eig(magic(3)))};

do_three_things()

This is the first thing.
This is the second thing.
ans =
[25]    [26]    [15]


We've executed three statements on a single line. All of the outputs are stored in the cell array, so we have three elements in the cell array. The first two outputs are actually garbage as far as we're concerned (they're just the outputs from fprintf, which is the number of bytes written, which we don't care about at all). The last output is from max(eig(magic(3))); That is, the biggest eigenvalue of magic(3) is exactly 15. Let's say we just wanted that final value, the eigenvalue. It's the third element of the cell array, so we can grab it with curly.

do_three_things = @() curly({fprintf('This is the first thing.\n'), ...
fprintf('This is the second thing.\n'), ...
max(eig(magic(3)))}, 3);

do_three_things()

This is the first thing.
This is the second thing.
ans =
15


For a more complex example, let's say we want to write a function to:

1. Create a small figure in the middle of the screen
2. Plot some random points
3. Return the handles of the figure and the plot

Then by storing all of the outputs in a cell array and using curly to access the outputs we care about, we can make a multi-line function with multiple outputs, all in a simple anonymous function.

dots = @() curly({...
figure('Position', [0.5*screen_size() - [100 50], 200, 100], ...
'MenuBar',  'none'), ...                % Position the figure
plot(randn(1, 100), randn(1, 100), '.')}, ...  % Plot random points
':');                                          % Return everything

[h_figure, h_dots] = dots()

h_figure =
3
h_dots =
187


(As a quick aside, note that if a statement doesn't return anything, we can't put it in a cell array, and so we can't use it this way. There are ways around this, discussed here.)

#### To Be Continued

Today, we've come a long way, from a simple condition through recursion and executing multiple statements. Here's a roundup of the functions so far.

map     = @(val, fcns) cellfun(@(f) f(val{:}), fcns);
mapc    = @(val, fcns) cellfun(@(f) f(val{:}), fcns, 'UniformOutput', 0);
iif     = @(varargin) varargin{2*find([varargin{1:2:end}], 1, 'first')}();
recur   = @(f, varargin) f(f, varargin{:});
paren   = @(x, varargin) x(varargin{:});
curly   = @(x, varargin) x{varargin{:}};


These can also be found here, implemented as regular MATLAB functions that can be kept on the path.

Next time, we'll look at loops. Until then, have you worked with functions such as paren or curly? How else are people implementing these or similar operations? Let us know here.

Get the MATLAB code

Published with MATLAB® R2012b

Loren Shure works on design of the MATLAB language at MathWorks. She writes here about once a week on MATLAB programming and related topics.

These postings are the author's and don't necessarily represent the opinions of MathWorks.