Jiro's pick this week is `roundsd` by François Beauducel.

There are a number of entries, including this, this, this, and this, that perform this useful task of rounding numbers to specific significant digits. François provides a simple, well-implemented function with good documentation. Here's how it works.

- rounded to 3 significant digits

```
format long
roundsd(pi,3)
```

ans = 3.140000000000000

- /100 rounded to 3 significant digits:

roundsd(pi/100,3)

ans = 0.031400000000000

It works with arrays as well.

roundsd([pi pi/100],3)

ans = 3.140000000000000 0.031400000000000

Did you know that in R2014b, the `round` function gives you the option to do the same thing? It provides the option of rounding to significant digits or to digits
respect to the decimal point.

- Round respect to the decimal point

round([pi pi/100],3)

ans = 3.142000000000000 0.031000000000000

- Round to significant digits

round([pi pi/100],3,'significant') % Reset format format short

ans = 3.140000000000000 0.031400000000000

Thanks, François (and others), for providing this functionality to folks using MATLAB versions prior to R2014b!

**Comments**

Give it a try and let us know what you think here or leave a comment for François.

Get
the MATLAB code

Published with MATLAB® R2014b

An email thread about a date-related design question had me thinking about terminology this week. *(Note: There's nothing about image processing in today's post. Sorry!)*

Specifically, I was thinking about the term *time zone*. In the U.S. you might see a date written this way:

November 21, 2004 1:02 PM (EST)

EST means *Eastern Standard Time* in the U.S. Or you might see a date written like this:

11/21/2004 8:30 (-04:00)

Many people refer to EST or -04:00 as the "time zone." However, because I work across the hall from the chief designer and implementer (Hi, P^2) of the new datetime suite of functions in MATLAB R2014b, I know that these are more properly referred to as *time zone offsets*.

A time zone is a much more complex concept than a simple offset from UTC. For example, a time zone encompasses both current and historical rules for using Daylight Saving Time.

The time zone offset for a given time depends not only on the time zone but also both the year and the day-of-year.

The dependence on the day-of-year is pretty obvious. In the fall season where I live, the time zone offset depends on whether the day is before or after the first Sunday in November.

d1 = datetime('October 15, 2014','TimeZone','America/New_York')

d1 = 15-Oct-2014

tzoffset(d1)

ans = -4:00

d2 = datetime('November 15, 2014','TimeZone','America/New_York')

d2 = 15-Nov-2014

tzoffset(d2)

ans = -5:00

Most of the state of Arizona does not observe Daylight Saving Time.

d3 = datetime('October 15, 2014','TimeZone','America/Phoenix')

d3 = 15-Oct-2014

tzoffset(d3)

ans = -7:00

d4 = datetime('November 15, 2014','TimeZone','America/Phoenix')

d4 = 15-Nov-2014

tzoffset(d4)

ans = -7:00

The dependence of time zone offset on the year is somewhat less obvious. It turns out that Daylight Saving Time rules change occasionally. In the U.S. the rules were last changed in 2007. Let's compare the time zone offset for noon on Halloween for 2007 and 2006.

d5 = datetime('October 31, 2007 12:00 PM','TimeZone','America/New_York')

d5 = 31-Oct-2007 12:00:00

tzoffset(d5)

ans = -4:00

d6 = datetime('October 31, 2006 12:00 PM','TimeZone','America/New_York')

d6 = 31-Oct-2006 12:00:00

tzoffset(d6)

ans = -5:00

Now for just a little fun with MATLAB R2014b date calculations. Most of the U.S. changes from Daylight Saving time to Standard Time at 2:00 AM on the first Sunday of November. Given the year, how can we compute that?

Here's just one way. Start with 2:00 AM on November 1 of the specified year.

year = 2014; d = datetime(year,11,1,2,0,0)

d = 01-Nov-2014 02:00:00

Now use `dateshift` to move to the first Sunday.

dateshift(d,'dayofweek','Sunday') % You can also use 1 instead of 'Sunday'.

ans = 02-Nov-2014 02:00:00

My horoscope for today said "Your next blog will include an example of an anonymous function for no good reason." So here it is. Let's make an anonymous function that returns the exact time for the transition from Daylight Saving Time to Standard time in the U.S.

f = @(year) dateshift(datetime(year,11,1,2,0,0),'dayofweek','Sunday')

f = @(year)dateshift(datetime(year,11,1,2,0,0),'dayofweek','Sunday')

f(2014)

ans = 02-Nov-2014 02:00:00

f(2015)

ans = 01-Nov-2015 02:00:00

f(2016)

ans = 06-Nov-2016 02:00:00

Of course, this function only works for years between 2007 and whenever the U.S. Congress changes its mind again!

Next week, I promise, I'll get back to writing about image processing. Or maybe colormaps. Unless I write about something else.

Get
the MATLAB code

Published with MATLAB® R2014b

Brett's Pick this week includes nods to two readers, and a submission called Fast Noise Estimation in Images by Tolga Birdal.

In my last post, back in September, I discussed several files for quantifying noise in images, without reference to a gold-standard image. I Picked and compared three submissions that simplified the quantification of image noise in isolated images. All of my tests, however, were conducted against image files that I had artificially blurred. Eric and Marco correctly took me to task for that, pointing out that 'blur' is not synonymous with 'noise' (though maybe we can agree that it can be a *type* of noise.)

Marco further pointed me to Tolga's submission, a fiendishly simple function that convolves the test image with a reference kernel to quantify noise.

In today's post, I would like to revisit noise quantification using my previous Picks, *and* Tolga's, and to look at their application to images containing different types of noise. I will again use the "rice.png" image that ships with the Image Processing Toolbox to consider how each does in quantifying the five types of noise-- "contant Gaussian," "locally varying Gaussian," Poisson, salt-and-pepper, and speckle--created by the `imnoise` function. To start, I'm going to use default input parameters (where possible) in the calls to `imnoise`.

img = cell(6,1); img{1} = imread('rice.png'); subplot(2,3,1) imshow(img{1});title('No Noise'); img{2} = imnoise(img{1},'gaussian'); subplot(2,3,2) imshow(img{2});title('Gaussian'); img{3} = imnoise(img{1},'localvar',0.05*rand(size(img{1}))); subplot(2,3,3) imshow(img{3});title('Local Gaussian'); img{4} = imnoise(img{1},'poisson'); subplot(2,3,4) imshow(img{4});title('Poisson'); img{5} = imnoise(img{1},'salt & pepper'); subplot(2,3,5) imshow(img{5});title('Salt & Pepper'); img{6} = imnoise(img{1},'speckle'); subplot(2,3,6) imshow(img{6});title('Speckle');

Now consider how each of the tools we've discussed:

quantifies noise in each of these images:

for ii = 1:6 imageBlur(ii) = blurMetric(img{ii}); blindImageQuality(ii) = blindimagequality(img{ii},2); noiseLevel(ii) = NLEstimate(img{ii}); fastNoiseEstimate(ii) = estimate_noise(img{ii}); end

NoiseTypes = {'NoNoise';'Gaussian';'LocalGaussian';'Poisson';'SaltAndPepper';'Speckle'}; Results = table(imageBlur,blindImageQuality,noiseLevelEstimage,fastNoise,... 'RowNames',NoiseTypes');

Interpreting this table is difficult and subjective. We can't necessarily compare metrics (which are on different scales), and we can't assess how each metric reflects different *levels* of noise. So let's do another series of analyses, increasing the noise successively, and recalculating the metrics. We can then plot some trends.

Let's increment four types of noise (omitting Poisson noise, for which `imnoise` is not adjustable) five times and recompute the metrics:

noNoise = imread('rice.png'); index = 0; for noiseType = 1:4 % Loop over noise types for noiseLevel = 1:5 index = index + 1; level = 0.01*noiseLevel; % Create noisy images and loop over metrics switch noiseType case 1 str = ['Gaussian (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'gaussian',0,level); case 2 str = [' Local Gaussian (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'localvar',level*rand(size(noNoise))); case 3 str = ['Salt & Pepper (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'salt & pepper',level); case 4 str = ['Speckle (' num2str(level) ')']; img{noiseType} = imnoise(noNoise,'speckle',level); end subplot(4,5,index); imshow(img{noiseType});title(str); imageBlur(noiseType,noiseLevel) = blurMetric(img{noiseType}); blindImageQuality(noiseType,noiseLevel) = blindimagequality(img{noiseType},2); noiseLevelEstimage(noiseType,noiseLevel) = NLEstimate(img{noiseType}); fastNoise(noiseType,noiseLevel) = estimate_noise(img{noiseType}); end end

f = figure; set(f,'DefaultAxesFontsize',8) NoiseTypes = {'Gaussian';'LocalGaussian';'SaltAndPepper';'Speckle'}; for noiseType = 1:4 ax = subplot(4,1,noiseType); noise = [Results{1,1},imageBlur(noiseType,:)]; plot(noise/max(noise)); if noiseType == 1 title('Normalized Noise Metrics','fontsize',12); end hold on noise = [Results{1,2},blindImageQuality(noiseType,:)]; plot(noise/max(noise)); noise = [Results{1,3},noiseLevelEstimage(noiseType,:)]; plot(noise/max(noise)); noise = [Results{1,4},fastNoise(noiseType,:)]; plot(noise/max(noise)); ax.XTick = 1:6; ax.XTickLabel = []; ylabel(NoiseTypes{noiseType}) end ax.XTickLabel = 0:5; xlabel('NOISE LEVEL') legend(MetricNames)

Conclusions are easier to draw now. For the most part, we can see metric values monotonically increasing or decreasing with the amount of noise. (The trend is less important than the monotonicity; we can easily think of the metric values as quantifying image quality.)

One final comment: it's very possible--maybe even likely--that these metrics can be tweaked to give better results, and that they will perform differently for different image and noise types. I encourage you to play around with the inputs using your own images; hopefully this discussion will give you a framework for evaluating the different noise metrics.

Swag, and thanks, are due to Eric, Marco, and Tolga!

As always, I welcome your thoughts and comments.

Get
the MATLAB code

Published with MATLAB® R2014b

**Reuse: Or how I learned to love the Bat Man Lego set…**

Like many an engineer of my generation I grew up with the basic Lego® starter set: the 1x1, 2x2, 2x3, 1x8 and the massive 2x8. We considered the wheel a high tech bit of plastic. Years later when I first saw a custom set with its detailed instructions, I was outraged by what I assumed was the corporate-sponsored crushing of my nephew’s creative possibilities. I scoffed, I scorned, I worried about the decline of western civilization if we fed these kids ships fully built!

My nephew had a whole raft of “complete” kits. From these kits he learned how to build houses and starships, bridges and mountains. What’s more he learned how to build in a way that was rugged and would survive the hard-play action of space ships crashing into erupting volcanos. One day I stopped and ** saw** what he had built — a hybrid submarine from bits of the Bat Man street cruiser and the Millennium Falcon brought together with good old basic Legos. He christened it the

*Figure 1 Not actual Subbatfalc*

**What's the link with Simulink?**

So what does this have to do with MATLAB and Simulink? Well, everything. As engineers creating something new, we need to figure out where we start. Do we use the primitive blocks and create everything from first principals or do we find existing components and build from there?

MATLAB and Simulink gives you a wonderful mixture of basic and advanced blocks. Further, with object-oriented programing in MATLAB and Simulink Libraries and Model blocks in Simulink, we give you the ability to create your own advanced blocks (heck, if you use TLC we give you the basic plastic to craft anything).

However as a consultant in the Controls Design Area I have lost track of the number of times I have worked with a client who has built up a set of custom functions such as transfer functions, integrators, or table look-up algorithms. They’ve recreated basic blocks. Now, truth be told, in many cases these custom blocks had a slight performance edge over the built-in blocks. But customers also have to maintain and validate these blocks, taking time away from building something more important.

So how do you determine when to build from scratch and when to use built-in functionality? I would suggest the following 6 State chart.

*Figure 2: Note remember to thank Sally!*

My nephew took less than an hour to build the Subbatfalc. Having worked with the kits, he had the skills to build it from basic blocks but it would have taken him ten times longer. He made the rational choice to reuse; and so am I. So, here I am recanting my harangue against kits. It isn’t the end of civilization — it’s the start of more interesting time. Now, pardon me, I am going to go help him build the Mount Doom/Death Star mash-up (Doom Star??) so we can see what happens when the Subbatfalc attacks!

**Now it's your turn**

How do you decide between using an existing complex block, versus building it from basic blocks? Let us know by leaving a comment here.

]]>The Method of Particular Solutions computes a highly accurate approximation to the eigenvalue we have been seeking, and guaranteed bounds on the accuracy. It also provides flexibility involving the boundary conditions that leads to the MathWorks logo.

In 1965-66, after finishing grad school, I spent a postdoctoral year at the ETH, the Swiss Federal Institute of Technology in Zurich. At the time, the faculty included three world renowned numerical analysts, Eduard Stiefel, Heinz Rutishauser, and Peter Henrici. Moreover, IBM had a branch of their research division in Ruschlikon, a Zurich suburb, where they held occasional seminars.

In October, Leslie Fox, from Oxford, visited the IBM center and gave a talk about the work two of his students, Joan Walsh and John Reid, were doing on the L-shaped membrane. Expansions involving Bessel functions that matched the singularity of the eigenfunction near the reentrant corner were coupled to finite difference grids in the rest of the domain. Shortly after Fox's talk, it occurred to me to use the Bessel function expansions over the entire domain and avoid finite differences altogether.

Within a matter of weeks I had done computations on the CDC 1604 mainframe at the ETH that had me very excited. To celebrate, I made an L-shaped ornament out of aluminum foil for the Christmas tree in our apartment.

Henrici contributed some theoretical background relating what I was doing to work of Stephan Bergman and a Russian mathematician I. N. Vekua. Together with Fox, we published a paper in the SIAM Journal on Numerical Analysis in 1967.

You don't need to know much about Bessel functions to appreciate the method of particular solutions. Using polar coordinates, a particular solution is any function of the form

$$ v_\alpha(\lambda,r,\theta) = B_\alpha(\sqrt{\lambda} r) \sin{(\alpha \theta}) $$

First of all, notice that the $\sin{(\alpha \theta)}$ portion insures that this solution satisfies the boundary condition

$$ v_\alpha(\lambda,r,\theta) = 0 $$

along the two rays with $\theta = 0$ and $\theta = {\pi}/{\alpha}$. For our L-shaped domain the reentrant corner has an angle of $\frac{3}{2}\pi$, so we will want to have

$$ \alpha = \frac{2}{3}j $$

where $j$ is an integer.

The Bessel function $B_\alpha$ is defined by an ordinary differential equation that implies any $v_\alpha$ satisfies the partial differential equation

$$ \Delta v_\alpha + \lambda v_\alpha = 0 $$

If $\alpha$ is not an integer, then $B_\alpha(r)$ has another important property. Its derivatives blow up as $r$ goes to 0. This asymptotic singularity is just what we need to match the singularity in the eigenfunction.

We can also take advantage of symmetries. The first eigenfunction is symmetric about the centerline of the L. We can also eliminate components that come from the eigenfunctions of the unit squares that are subdomains. This implies that for the first eigenfunction $\alpha = \frac{2}{3}j$ where $j$ is an odd integer that is not a multiple of 3, that is

$$ \alpha_j = \frac{2}{3}j, \ \ j = 1,5,7,11,13, ... $$

Form a linear combination of particular solutions.

$$ u(\lambda,r,\theta) = \sum_\alpha c_\alpha v_\alpha(\lambda,r,\theta) $$

This function satisfies our differential equation throughout the domain, satisfies the boundary condition on the two edges that meet at the reentrant corner, and has the desired singular behavior there. All we have to do is pick $\lambda$ and the coefficients $c_\alpha$ so that $u$ satisfies the boundary condition, $u = 0$, on the remaining edges of the L.

The general idea is to pick a number of points on the outer boundary and investigate the matrix obtained by evaluating the particular solutions in the linear combination at these points. Vary $\lambda$ until that matrix is nearly rank deficient. The resulting null vector provides the desired coefficients.

The details of how we carry out this general plan have varied over the years. The best approach is implemented in the function `membranetx`, available here, described in the PDE chapter of *Numerical Computing with MATLAB*. It uses the matrix computation workhorse, the SVD, the singular value decomposition.

The number of rows in the matrix involved is the number of points on the boundary and the number of columns is the number of terms in the sum. The figures below come from matrices that are 46-by-9. Let $(r_i,\theta_i)$ be the points chosen on the boundary. The matrix elements are

$$ a_{i,j}(\lambda) = B_{\alpha_j}(\sqrt{\lambda} r_i) \sin{(\alpha_j \theta_i)} $$

The function `fmintx` is used to vary $\lambda$ and find a minimum of the SVD of this matrix. The resulting matrix is essentially rank deficient and the resulting right singular vector is a null vector that provides the coefficients.

The minimum is found at

$$ \lambda = 9.639723844 $$

Here are the first four particular solutions. The vertical scale reflects the magnitude of the corresponding coefficient.

mps_figs

None of these four terms is remotely close to zero on the outer boundary, but when they are added together the sum is zero to graphical accuracy.

mps_figs(4)

So this is actually the first eigenfunction of the L-shaped membrane.

Now we can have some fun. Here is a `surf` plot of the sum of just the first two particular solutions.

mps_figs(2)

This function is an exact solution of the partial differential equation

$$ \Delta u + \lambda u = 0 $$

with the value of $\lambda$ reported above. It satisfies the boundary condition $u = 0$ on the two edges meeting at the reentrant corner, but it fails significantly on the outer edges. So it is only with some artistic license that we can refer to it as the first eigenfunction of L-shaped membrane.

When I first saw this `surf` plot, it didn't have this gorgeous color map, but I found the shape so attractive that I kept it around. The evolution of the shading, lighting and colors that turned this into our logo is the subject of my next blog.

The Fox, Henrici and Moler paper was not just about the method of particular solutions. We were also concerned about finding guaranteed upper and lower bounds for the exact eigenvalue of the continuous differential operator. A theorem proved in the paper relates the discrepancy in meeting the boundary condition to the width of a small interval containing the eigenvalue. In 1966 we proved

$$ 9.6397238_{05} < \lambda < 9.6397238^{84} $$

Today with `membranetx` on MATLAB, I am finding a value of $\lambda$ in the center of this interval.

```
[L,lambda] = membranetx(1);
format long
lambda
```

lambda = 9.639723844021997

In 2005 Timo Betcke and Nick Trefethen published a paper in *SIAM Review* that pointed out a serious shortcoming in the approach we described in the Fox, Henrici and Moler paper. When the method is attempted on other domains, especially domains with more than one singular corner, the matrix has a tendency to appear nearly rank deficient for all values of $\lambda$. It is difficult to identify local minima. The recommended cure is to add a number of interior points and require the approximate eigenfunction to be nonzero there.

Betcke and Trefethen also prove an approximation theoretic result that, when combined with their computations, adds four more digits to our rigorous bounds on the first eigenvalue.

$$ 9.63972384402_{16} < \lambda < 9.63972384402^{22} $$

Moreover, they are so confident in their numerical results that they are sure the exact eigenvalue lies in the center of this interval, and they are pretty sure they even know the next digit.

$$ \lambda \approx 9.63972384402195 $$

This agrees with what I get out of `membranetx` to almost full double precision accuracy. I think we finally have this fellow pinned down.

I recommend that anyone reading this blog who wants to know more about the method of particular solutions and computing the eigenvalues of the L-shaped membrane read the Betcke and Trefethen paper.

Leslie Fox, Peter Henrici, and Cleve Moler, Approximations and Bounds for Eigenvalues of Elliptic Operators, SIAM Journal Numerical Analysis 4, 89-102, 1967, <http://blogs.mathworks.com/images/cleve/Fox_Henrici_Moler.pdf>

Timo Betcke and Lloyd N. Trefethen, Reviving the Method of Particular Solutions SIAM Review 47, 469-491, 2005, https://people.maths.ox.ac.uk/trefethen/publication/PDF/2005_112.pdf

Cleve Moler, *Numerical Computing with MATLAB*, Chapter 11, Partial Differential Equations, 2004. <http://www.mathworks.com/moler/pdes.pdf>

Get
the MATLAB code

Published with MATLAB® R2014b

Sean's pick this week is votebar by paul koch.

Paul's `votebar` is a nice utility for viewing cylindrical three dimensional bar plots like you would see after an election. For example:

```
n = 4;
pct = randfixedsum(n,1,100,0,100); % File Exchange
votebar(100, pct, parula(n), 0);
```

Or in the case of the "Is my wife right?" CODY problem:

votebar(100,[100 0],[0 1 0; 0 0 1],0); text(1.5,-1,'Yes','FontSize',20) text(1.5,1.5,'No','FontSize',20) title('Is My Wife Right?')

Let's take a look at this with the election data from the 2014 Massachusetts gubernatorial race.

% Data and setup votes = [1041640 1001279 71144 19192 16125]; pct = votes./sum(votes)*100; names = {'Charlie Baker', 'Martha Coakley', 'Evan Falchuk', 'Scott Lively', 'Jeff McCormick'}; colors = [1 0 0; 0 0 1; 0.541 0.706 0.486]; % RGBcolors idx = [1 2 3 3 3].'; % Index into colors % Votebar hp = votebar(100, pct, colors(idx,:), 0); % Decorating legend(hp(3,:),names) lighting gouraud material shiny title('Massachusetts Gubernatorial Results') view(84,14) text(-6.5,-2,'%','FontSize',20)

Source: *http://www.politico.com/2014-election/results/map/governor/massachusetts/#.VGIlaE1OXcw*

Give it a try and let us know what you think here or leave a comment for Paul.

Get
the MATLAB code

Published with MATLAB® R2014b

I hesitated for some time today about where to go next in my discussion about the new default colormap in MATLAB. I finally decided to spend a little time on the reactions so far. Those reactions have ranged from "Great, what took you so long?" to "I hate it; give me back jet!" In the Twitterverse, one scientist proclaimed "End of an era!" Another simply said "yuck." (Jet is still in MATLAB, by the way.)

Today I want to dig into some of those reactions a little bit. Along the way I'll show some experiments with a test pattern to explain some of my comments about jet.

To remind you, last time I summarized 15 years or so of criticisms about rainbow colormaps:

- Rainbow colormaps confuse viewers because there is no natural perceptual ordering of the spectral colors.
- Rainbow colormaps obscure small details in the data.
- Rainbow colormaps mislead viewers by suggesting data features that are not really there. These "phantom features" often take the form of false boundaries. This effect can falsely segment the data.
- Rainbow colormaps lose critical information about high and low data values when displayed or printed on a gray-scale device.
- Rainbow colormaps can be difficult to interpret for some color-impaired viewers.

I have seen a claim that that rainbow incorporates a universally understood order based on associating blue with cold and red with hot. There's no doubt that this order makes sense for some people (although blue flame is hotter than red flame?). I quibble with the "universally understood" part, though. The dominant evidence is that individuals do not universally order the spectral hues the same way unless they refer to spectral order mnemonic, such as "Roy G. Biv." It is well-established, however, that individuals can easily and consistent order colors by relative brightness.

Another commenter said that the rainbow colormap has higher contrast than parula. Yes, this is true. Rainbow goes all the way from dark to bright in the bottom third of the data range, and it goes all the way from bright to dark in the top third of the data range. For some parts of the data range, therefore, there is higher visual contrast than with parula, which uses the entire data range to go from dark to bright. A price paid for this increased contrast, however, is the visual illusions and confusion caused by the reversed roles of dark and bright in different parts of the colormap.

Matteo Niccoli gives the following nice example of a rainbow-based visual illusion. The arrow in the right-hand image below points to something that gives the strong but false impression of being a depression in the bone structure.

*Image credit: Matteo Niccoli, "The Rainbow Is Dead ? Long Live the Rainbow!" Copyright 2014 MYCARTA. Used with permission.*

And I showed another visual confusion example last time with these visualizations of two different sinusoidal patterns:

These visualizations give the impression that the two sinusoidal patterns have opposite phase, but in reality they have exactly the same phase.

Another commenter said that, with sufficiently familiarity, the bright "stripes" that are characteristic of the rainbow colormap can be an advantage for certain kinds of data. That might be true in some particular cases, but I don't think it's true in general. Let me show an example, using a test pattern, that illustrates how much difficulty the cyan and yellow stripes can cause in interpreting what should be easily recognizable data.

```
url = 'http://blogs.mathworks.com/images/steve/2014/eia-2000-cropped.png';
I = imread(url);
imshow(I)
colormap(jet)
```

Let me label some portions of the visualization and then ask some quetions about it.

text(35,90,'A','FontSize',25,'HorizontalAlignment','center',... 'BackgroundColor','white') text(125,90,'A','FontSize',25,'HorizontalAlignment','center',... 'BackgroundColor','white') text(215,90,'A','FontSize',25,'HorizontalAlignment','center',... 'BackgroundColor','white') text(190,320,'B','FontSize',25,'HorizontalAlignment','center') line([77 160],[225 225],'Color','white','LineWidth',5) text(65,225,'C','FontSize',25,'HorizontalAlignment','right',... 'BackgroundColor','white')

At periodic intervals, the dark blue vertical bars have a bright outline (see the points marked "A"). What is the data significance of these bright outlines? Especially when compared with the dark blue vertical bars that do not have a bright outline?

The vertical bars in between the "A" marks don't have bright outlines. Those bars alternate between dark blue and dark red, a low-contrast combination that is a little difficult to focus clearly on (a red-blue visual illusion called Chromostereopsis).

Above the "B", I see very thin vertical dark bars. When the bars are this thin, it's difficult to distinguish between the dark red and dark blue, which is unfortunate because those colors represent opposite ends of the data range. The bright outlines between the bars is broken up by a diagonal pattern of blue stripes.

The numbers around the border are oddly broken up and hard to decipher.

Let's investigate why some vertical bars are outlined and some not by looking at the data cross-section along the white path marked "C".

```
c = improfile(I,[77 160],[225 225]);
figure
plot(c,'*-')
```

I'll superimpose two lines that show the data values corresponding to the bright cyan and bright yellow in the jet colormap.

hold on plot([1 84],[1 1]*0.375*255,'c','LineWidth',5) plot([1 84],[1 1]*0.625*255,'y','LineWidth',5) hold off

It is now apparent why some vertical bars are outlined but not others. In each transition from low to high and from high to low, there is one data value that is in between. When that one transitional data value happens to be close the data values mapping to cyan or yellow, the vertical bar is shown with a bright outline. Otherwise, there is no outline. Although the effect is visually dramatic and makes some regions of the data look quite different from others, the visual difference has no real data significance.

OK, now let's show this data using the parula colormap.

imshow(I) colormap(parula)

This visualization is trivially easy to interpret. The vertical bars in the middle are all the same. The digits are easy to read. The visual difference between low data values and high data values has high contrast and is easy on the eyes. And the data above the point marked "B" is seen to have a completely different form than the way it appeared with jet. The stripes have a lower spatial frequency and are seen to be slightly fanning out in direction instead of appearing vertical.

A valid complaint about this example is that it is completely artificial, with data rather unlike what people are likely to be looking at with jet or parula. I'd like to suggesting thinking about it this way, though: If jet can take familiar, easily interpretable data and render it almost unrecognizable, what can jet do to your own data?

Finally, I'll mention that a blog reader commented that using jet made it easy to distinguish positive and negative values from each other in a data set whose values are centered around zero. And yes, jet can be better than parula for that. Jet is bright and relatively unsaturated in the middle, and it is dark and saturated at both ends. These are characteristics often shared by colormaps designed specifically for highlighting a central, specific value and distinguishing the values above and below that. These are called *divergent* colormaps. In a future post I'll show where to find some good divergent colormaps and how to use them in MATLAB.

Get
the MATLAB code

Published with MATLAB® R2014b

**Enabling Interface Display**

In R2014b, you will notice a new entry at the top of the Display menu:

For a model that originally looks like this:

Enabling the interface display makes it look like this:

As you can see, this feature puts the emphasis on the inputs and outputs of the system.

**Advantages of the Interface Display**

For someone like me who receives new large Simulink models every day, this feature is very useful. When I explore new models, I can use the highlighting to trace sources and destinations of signals through multiple levels of nested subsystems. As you can see in the following animation, by zooming in, I can select one element of a bus signal and follow the highlighted signals to quickly find out where it is used.

Once the model has been simulated or updated once, the interface view displays for each port its data type and sample time, making it easier to understand the model.

One more thing to note. When in interface Display, you cannot edit the model. This is convenient when your goal is just to explore and understand the model without modifying it.

**Now it's your turn**

Give the new interface view a try and let us know what you think by leaving a comment here.

]]>Today, David Garrison, our guest blogger, will continue his series on the new graphics system in R2014b.

- Part 1: Features of the New Graphics System
- Part 2: Using Graphics Objects
**Part 3: Compatibility Considerations in the New Graphics System**

Here is Part 3 of the series.

- What have we learned so far?
- Section 1: Compatibility Considerations for R2014b graphics
- Why did we make these changes?
- How do the changes affect me?
- How do I find out more about the changes in R2014b graphics?
- Section 2: Visual Differences
- The New Default Colormap
- Line Colors in Plots
- Plot Titles and Labels
- Section 3: Changes that Affect Advanced Graphics Users
- Graphics Functions Return Objects, not Numeric Handles
- Colorbars and Legends are No Longer Axes
- Objects Returned by Certain Charts Have Changed
- UI Controls May Not Appear in a GUI
- What do you do now?
- Have you encountered any incompatibilities?
- Conclusions

In Part 1 of this series, I provided an introduction to the new MATLAB graphics system in R2014b with descriptions of a number of new features. In Part 2, I talked about one of the big changes in R2014b -- *graphics functions return MATLAB objects, not numeric handles*. In this post, I will talk about compatibility considerations in using R2014b graphics.

This post is longer than my two previous posts. In the first section, I'll talk about the kind of changes we've made in R2014b graphics, why we made those changes, how they affect you, and where to go for more information. The second section describes visual changes in the new graphics system including a new colormap, new line colors in plots, and new sizes for plot titles and axis labels, and these affect everyone. The last section is about changes we've made that affect advanced graphics users and user interface builders. Most of you can skip the more advanced material in the last section.

When we talk about compatibility considerations we are talking about changes in R2014b that may not be compatible with code written before R2014b. These *incompatibilities* can affect your code in one of three ways.

- A visual change that affects how a plot looks
- A change that affects the behavior of your code
- A change that causes an error that will require you to modify your code

First of all, let me say that here at MathWorks we take compatibility very seriously. We have a rigorous process for vetting potential changes in MATLAB that may affect existing MATLAB code. We evaluate them very thoroughly to determine the extent of their impact. We document those changes in our release notes and provide alternatives whenever possible. Those changes are necessary to allow us to continue to develop MATLAB.

One of our goals for the new graphics system was to introduce a new architecture that would allow us to build new graphics capabilities for the future. We couldn't continue with the old graphics architecture because it was too limiting. The new architecture has also allowed us to fix a large number of outstanding graphics bugs that have accumulated over the years.

Our second goal was to support existing graphics functionality and to minimize the disruption to MATLAB users. We did that by evaluating every change that would introduce an incompatibility in the graphics system. We tested those changes in house with over 100,000 MATLAB files and in person with a large number of advanced graphics users to determine the following:

- How often will a particular incompatibility affect user code?
- Which users will be impacted?
- How hard will it be to identify the effect of a specific incompatibility?
- How much work will be required if code changes have to be made?

The results of all this testing allowed us to do two things. First, it allowed to make changes in the new graphics system before the release to minimize the effect of certain incompatibilities. Second, it told us what kind of documentation and tools would be required to help people make the transition.

Anyone using graphics will notice a visual difference in R2014b. People who've seen R2014b tell us the new visual appearance is a big improvement over previous versions of MATLAB. Differences include changes to colormaps, line colors, font sizes, grids, and other visual properties. **For most people that's it. That's all you need to know.**

Those of you who do advanced graphics programming and build user interfaces may be affected by incompatibilities that will require changes to your code. I'll talk about the most common ones in the last section of this post. We've also provided plenty of documentation and some tools to help you make the transition.

In preparation for the release of the new graphics system, we created a lot of material to help people understand the changes and make any modifications to their code. Here is a list of resources for information on R2014b graphics:

- MATLAB Graphics in R2014b: Learn about the new graphics system with links to information about new features, compatibility, help and support, and R2014b graphics examples.
- MATLAB Release Notes: Read detailed information about changes in the
**Graphics**and**GUI Building**sections. - Graphics Changes in R2014b: Read about changes to R2014b graphics and how to troubleshoot your code.
- Tools for transitioning to R2014b Graphics: Download tools to check for incompatibilities and find missing graphics in GUIs.
- Graphics Performance: Find techniques to maximize graphics performance.
- System Requirements for Graphics: Troubleshoot issues associated with graphics cards and drivers.
- MATLAB Answers: Find answers to common questions about using R2014b graphics.

Anyone who uses MATLAB graphics in R2014b will see that plots look quite different. In this section, I will show you some of the big visual differences. As we go through these examples, I think you'll agree that the changes in R2014b are a big improvement over previous MATLAB versions. Of course, there may be times when you may want to change some of these visual properties. I will either tell you how to do that or I will provide a link to the documentation for more information.

In MATLAB, the `colormap` is a matrix of RGB color values used to set the colors for images and surfaces. For example, when you call the `surf` function you get a plot where the surface color is proportional to the height of the plot. Colors are based on the values from the colormap.

In R2014b, we introduced a new colormap called *parula* which is now the default. In previous versions of MATLAB, the default colormap was *jet*. Here's a comparison of the colormaps between R2014a and R2014b.

We changed the default colormap to address problems with rainbow colormaps like *jet*. The new *parula* colormap is ordered from dark to light and is perceptually uniform. For more information about the issues with rainbow colormaps like *jet*, see the post on Steve Eddins blog called A New Colormap for MATLAB – Part 2 – Troubles with Rainbows. You can use the `colormap` function to change the colormap to whatever you want it to be. For example, to change the colormap to *jet* use the following command after you create your plot.

colormap(jet)

We've also changed the lines colors used for plots. These line colors were selected to make lines easier to distinguish and to help people with certain types of color blindness.

The line colors used in plots are controlled by the `ColorOrder` property of the `Axes` object. One other note regarding line colors in plots. In R2014b, when plot is called with `hold on`, MATLAB will use the next color in the color order for each call to `plot`. That gives plots with different line colors. In previous versions, each call to plot would start over in the color order so that each line would be the same color.

For information on controlling the colors in the ColorOrder see the documentation section Why Are Plot Lines Different Colors?

One last change regarding appearance of plots. In R2014b, plot titles and axis labels use a larger font size than in previous versions of MATLAB and plot titles are bold by default.

I intentionally used a very long title in this plot to make a point. In this case, the title is so long it extends beyond the boundaries of the axes. You can control your plot title and axis labels using the `TitleFontWeight`, `TitleFontSizeMultiplier`, and `LabelFontSizeMultiplier` properties. For more information see How Do I Make the Graph Title Smaller?

That covers the major visual changes for plots in R2014b. The next section covers changes that only affect advanced graphics users and GUI builders.

*Most of you can stop reading here and get on with the rest of your day.*

If you are still reading, I assume you do some advanced graphics programming in MATLAB or maybe you build MATLAB user interfaces. In this section I'm going to cover the four most common non-visual changes in R2014b graphics that you are likely to encounter. If you do encounter one of these, you will probably have to make changes to your code.

I mentioned this briefly in Part 2 of this series. Graphics functions now return objects, not numeric handles. The R2014b documentation has detailed information about this subject in the section called Graphics Handles Are Now Objects, Not Doubles. I will give a couple of simple examples to illustrate what can happen with code written before R2014b.

Prior to R2014b, you could store a set of handles to graphics objects in an array and then add some numeric data to that array. In R2014b, that will cause an error.

```
x = -pi:0.1:pi ;
y1 = sin(x);
y2 = cos(x);
myLines = plot(x,y1,x,y2) % plot returns an array of two Line objects
```

If you then try to set `myLines(3) = 1.2`, you get the following error.

Cannot convert double value 1.2 to a handle

MATLAB won't let you add numeric values to an array of graphics objects. A similar problem occurs if you try to use an object handle in a function where MATLAB expects a numeric value. A simple example of this happens with the `sprintf` function.

```
a = sprintf('You clicked on figure %d\n', gcf);
```

The `%d` specification in the `sprintf` format string expects an integer value. However, since gcf is a `figure` object, you get the following error.

Error using sprintf Function is not defined for 'matlab.ui.Figure' inputs.

Here is one final example. Because graphics handles used to be numbers, you could use them in logical expressions.

if (get(0, 'CurrentFigure')) disp(['Figure ' get(gcf, 'Name')']) % display the figure name for gcf else disp('No open figures') % there is no open figure end

This worked in earlier versions of MATLAB because `get(0,'CurrentFigure')` would return either an empty array or a numeric figure handle. Both of these values are valid in the logical test of the `if` statement above. In R2014b, this will cause an error.

Conversion to logical from matlab.ui.Figure is not possible.

We have tried to maintain compatibility with previous releases in some cases. For example, you can still use `0` to refer to the graphics root in functions like `get` and `set`. As a best practice, however, we now recommend using the `groot` function to get the graphics root. Similarly, we still support the use of literal integer values to refer to figures in functions like `set`, `get`, and `figure`. Again, the best practice is to use a variable which contains the object when using these functions.

If you find yourself really stuck, it is possible to cast object handles to numeric handles using the `double` function. You can then cast the number back to an object handle using the `handle` function. We don't recommend this as a long term solution. Be aware that we may choose to remove this feature in a future version of MATLAB. If we do, we'll let you know in advance.

In earlier versions of MATLAB, the `legend` and `colorbar` functions created objects whose type was 'axes'. Technically, legends and colorbars were subclasses of axes. In R2014b, the `legend` function creates a Legend object and the `colorbar` function creates a ColorBar object.

So what does that mean? First, the results of a findobj call in R2014b may be different than in previous versions. In R2014b, the command

findobj('Type', 'Axes')

will not return any legend or colorbar objects. To get those objects you will need to do the following:

findobj('Type', 'Legend') findobj('Type', 'ColorBar')

Second, legends and colorbars cannot become the current axes. Prior to R2014b, you could write code that looks like this:

plot(1:10) cb = colorbar; axes(cb) % Make the colorbar the current axes title('My Colorbar') % set the title of the colorbar

If you try to run this code in R2014b you will see an error.

Error using axes Handles of type ColorBar cannot be made the current axes

In R2014b, you can still use the `title` function. Just give it the colorbar handle as the first argument.

```
title(cb, 'My Colorbar')
```

Alternatively you can use the colorbar's `Label` property to get the text object for the label and set its string property.

```
cb.Label.String = 'My Colorbar'
```

Charting functions like `bar`, `contour`, `stem` and others return one or more graphics objects. In previous versions of MATLAB, those objects had names like *barseries*, *contourgroup*, and *stemseries*. These objects were each a special type of object called an *hggroup*. Each of these hggroup objects had a set of children. The children were low level objects like lines or patches. Here is an example from R2014a.

[~,h] = contour(peaks, 'LineLevels', -6:1:8) ; get(h, 'Type') % an hggroup object is returned by the contour function

ans =

hggroup

ch = get(h, 'Children') ; get(ch(1), 'Type') % the children of the hggroup are patch objects

ans =

patch

Prior to R2014b, the `contour` command returned an hggroup object that had children that were patch objects (one for each contour line). Using the patch objects, it was then possible to do interesting things with the contour lines. You could, for example, make the even numbered contour lines solid and make the odd numbered contour lines dashed as shown below.

As you've probably guessed, things are different in R2014b. First, the types of the objects returned by these functions are different. For example, the `bar` function returns a *Bar* object and the `contour` function returns a *Contour* object. Second, the objects returned by these functions no longer have any children. That means you can't get to the low-level objects.

So what do you do? You have to take a different approach. In the code below, I've create two contours -- one with solid lines and one with dashed lines. This code will create a plot like the one shown above. It works in R2014b *and* in previous versions of MATLAB.

major = -6:2:8; minor = -5:2:7; [~,hmajor] = contour(peaks,'LevelList',major); % contour with even-numbered levels

hold on [~,hminor] = contour(peaks,'LevelList',minor); % contour with odd-numbered levels set(hminor,'LineStyle',':') % make the odd-numbered levels dotted hold off

There is one last change that I want to talk about. It can affect MATLAB GUIs created before R2014b. Suppose I have a simple GUI with a panel and a button. I might have code that looks something like this:

figure('NumberTitle', 'off', 'Name', version('-release'), ... 'Position', [100 100 350 260], 'MenuBar', 'none', 'Toolbar', 'none') ; button = uicontrol('Style', 'pushbutton', 'String', 'My Button', ... 'Units', 'normalized', 'Position', [0.4 0.5 0.25 0.15]) ; panel = uipanel('Position', [0.10 0.10 0.8 0.8], 'Title', 'My UIPanel') ;

Seems OK right? Here's what it looks like in R2014a and R2014b.

So what happened here? I don't see my button in R2014b. Where did it go? It turns out that the button is still there but in R2014b it is drawn under the panel.

In previous versions of MATLAB, uicontrols were always drawn on top on uipanels regardless of the order in which they were created. In R2014b, the two components are drawn in creation order. Since the panel was created after the button, it is drawn on top. What I really want is to have the button be a child of the panel. In order to do that, I can change my code to look like this:

figure('NumberTitle', 'off', 'Name', version('-release'), ... 'Position', [100 100 350 260], 'MenuBar', 'none', 'Toolbar', 'none') ; panel = uipanel('Position', [0.10 0.10 0.8 0.8], 'Title', 'My UIPanel') ; button = uicontrol('Parent', panel, 'Style', 'pushbutton', 'String', 'My Button', ... 'Units', 'normalized', 'Position', [0.4 0.5 0.25 0.15]) ;

You can also fix this problem in GUIDE by moving the button out of the panel and then moving it back in. This will automatically make the button a child of the panel.

You can check your existing GUIs by downloading a tool from Tools for transitioning to R2014b Graphics. For more information, see Why Are Some Components Missing or Partially Obscured?

Read the materials I've listed above to learn more about the changes to graphics in R2014b and how they might affect your MATLAB code. Download the tools and run them on your code and your GUIs. If you aren't able to find the information you need to work around one of the changes in R2014b, contact our Technical Support Team . They are very familiar with the changes in R2014b and can help you find a solution.

Have you encountered any of these changes in R2014b? Have you been able to find the information you need? Have you been able to use this information to make necessary changes to your code? Have you tried the transition tools? We'd love to hear your thoughts here.

I am hopeful that this series of posts has given you a good introduction to the R2014b graphics system. There's a lot of new stuff to digest in this release. Try these things out and let me know what you think.

Thanks to Loren for letting me take over her space for the last few weeks.

Get
the MATLAB code

Published with MATLAB® R2014b

The Partial Differential Equation Toolbox contains tools for the analysis of PDEs in two space dimensions and time. It is perhaps not surprising that one of the primary examples involves the L-shaped membrane.

If you have the PDE toolbox installed, bring up `pdetool`. Click on the first triangle in the toolbar. This initializes a tool for creating a planar domain with an unstructured mesh. The default domain is our L. The default mesh has equally spaced points on the boundary connected by an irregular grid in the interior with 150 nodes and 258 triangles.

The PDE tab allows the specification of a general elliptic operator with variable coefficients, but if we accept the defaults we have our desired Laplacian. Here is the first eigenfunction, obtained with the coarse grid, plotted with the default `cool` colormap. The first eigenvalue, reported in the title, is 9.9707.

Refine the grid three times to one with 8417 nodes and 16512 triangles, change to our new `parula` colormap, and add contour lines. The reported eigenvalue is 9.6525.

The `pdetool` has an option to plot the absolute value of the gradient. We see that the singularity at the origin acts like a black hole, sucking in all the color.

The mesh refinement is accomplished by adding grid points at the midpoints of the triangles, thereby quadrupling the number of triangles. It is possible to do this five times, producing a fine mesh with `258 * 4^5 = 264192` triangles. (Actually the mesh can be refined six times to one with over a million triangles, but the subsequent eigenvalue calculation runs out of memory.)

Here are the results for the first eigenvalue of the L-shaped membrane obtained with successive bisections of the unstructured triangular mesh.

format long load pdetool_results L

L = 9.970747465677787 9.745769128365856 9.675647712787020 9.652453140975609 9.644395207950412 9.641482807142362

These values are candidates for acceleration by Aitken's delta-squared process.

$$ x_n - \frac {(\Delta x_n)^2}{\Delta^2 x_n} $$

d1 = diff(L,1); L2 = L(3:end) - d1(2:end).^2./diff(L,2)

L2 = 9.643895738979518 9.640988739828559 9.640105597452624 9.639834371546435

We can even use delta-squared a second time, although this is not often justified.

t1 = diff(L2,1); L2(3:end) - t1(2:end).^2./diff(L2,2)

ans = 9.639720224107141 9.639714153353552

Compare this with the value I reported in my previous post, obtained by extrapolating from a regular square grid.

lambda1 = 9.639723753239963

We will see in my next post a much better way to compute the eigenvalue that will provide a result to nearly full double precision accuracy.

Get
the MATLAB code

Published with MATLAB® R2014b

*by Wendy Fullam*

Have you ever found yourself thinking, “I wish there was a place where I could go to see a bunch of MATLAB examples…”?

MathWorks has just what you’re looking for. We just launched a big new feature on the Support page: MATLAB Examples. On it, you can discover thousands of code examples for MATLAB and Simulink: data manipulation, signal processing, machine learning, statistics, finance, computational biology, finance, you name it. Maybe you’ll even discover a few things that you didn’t know were possible.

MATLAB Examples is a single destination to find high-quality code examples, including those authored by both MathWorks staff and contributors to the File Exchange. Here’s what the main page looks like.

What can you do here?

You can browse by topic area:

You can search by keyword, and filter the results:

After finding an example that you want to use, you can quickly see what product contains that example:

or access the code:

You can also find related examples:

The Help Team at MathWorks hopes you enjoy the new interface and welcomes your thoughts and reactions. Check out the new area and then let us know what you think in the comments below!

]]>How often do you find a File Exchange submission accompanied by a video tutorial? I've only explored a small fraction of all the content on the File Exchange, but in my limited sampling, I have to say that this is the first contribution with such an achievement. The model is simple enough: a prediction of the water level of two tanks as they are filled and drained. But as John notes in his video, the model doesn't jump up and scream "I'm clearly two draining tanks!" when look you at it.

As I reviewed the model, other thoughts came to mind on how to improve it. The tank levels are calculated using a Level- 1 MATLAB S-Function. This technique is so old that MathWorks doesn't even provide a documentation page for it. The only thing available is a discussion on how to maintain them. I recommend avoiding the S-Function and instead rely on a MATLAB Function block for the derivative calculation combined with an Integrator block. Or better yet, the whole thing could be done with Simuink blocks without MATLAB code. I'd also suggest eliminating all the Muxes and Demuxes in favor of separate signal lines or buses. So there's room for improvement, but a solid model nevertheless.

In my mind, it's material like this that the File Exchange was meant for. John combines MATLAB, Simulink, Excel, and YouTube together into an interactive educational experience. It's a lot of fun, and I hope you check it out.

Let us know what you think here or leave a comment for John.

]]>

Sean's pick this week is Display progress, pause or stop a time-consuming loop by Rafael Oliveira.

It is not uncommon to see questions such asking "How do I stop a `for` or `while` loop when a user pushes a button?" in MATLAB Answers. There are a few ways to do this, but usually the results involve checking the state of a toggle button or some flag that
has been set by a user. Alternatively, one could consider replacing the `for` or `while` loop with a `timer` which is what I typically recommend.

Rafael's `stop` makes stopping a loop easy and elegant to do.

for ii = 1:100 if ~stop.requested disp(['Still Running, iteration: ' num2str(ii)]) % Computationally expensive operation. fft2(rand(5000)); else % I clicked said stop disp('Stopped') break end drawnow end

Still Running, iteration: 1 Still Running, iteration: 2 Still Running, iteration: 3 Still Running, iteration: 4 Still Running, iteration: 5 Still Running, iteration: 6 Still Running, iteration: 7 Stopped

This utility also allows you to pause and resume elegantly.

I typically find File Exchange submissions either by randomly browsing, needing some piece of functionality, or trying to help someone not reinvent the wheel. However, I came by this one a little differently. Out of curiosity, I wanted to know which File Exchange submissions contain the string "xkcd". This one came up with a nice little easter egg related to a specific comic.

Turning that mode on:

for ii = 0:0.3:1 stop.setStatus('xkcd',ii) drawnow snapnow end

What approaches do you usually take to allow for user interactions with a loop or timer?

Give it a try and let us know what you think here or leave a comment for Rafael.

Get
the MATLAB code

Published with MATLAB® R2014b

If you're going to be there, come see my *MATLAB Techniques for Image Processing* talk on Wednesday, 12:45 - 13:45. Or stop by the MathWorks booth.

Now you'll have to excuse me while I finish preparing my talk.

]]>**Scrolling**

The Scope block now has a new **scroll** option to control the behavior when the duration of the simulation is longer than the Scope time range. The best way to describe that is to show it:

**Signal Viewer and Floating Scope Look and Feel**

In R2014b, the Signal Viewer and Floating Scope interfaces are now identical to the Simulink Scope block. Whether you prefer having a block or not in the model, you can experience the same visualization:

This means that the Viewer line properties, axes and figure colors can now be customized. Here is an example of how I like to configure it:

**Signal Viewer supports multiple new features**

For a few releases, the Signal Viewer and Floating Scope have lagged behind with many recent Simulink features. Not anymore! In R2014b, you will notice that the Signal Viewer and Floating Scope now offer features similar to the Scope block and support:

- Logging signal in dataset format
- Simulation Stepping
- Fast Restart

**Now it's your turn**

Do you prefer to use the Scope, Floating Scope, or Viewer? Why? Let us know by leaving a comment here.

]]>The first one I saw was Cleve's follow-up post on the eigenvalues and eigenfunctions of the L-shaped membrane, the basis of the MathWorks logo.

The last image from Cleve's post is a lovely contour plot of one of the eigenfunctions of the shape, rendered using the new parula colormap I've been discussing this month.

The second was the hilariously over-the-top high-tech Friday afternoon cookie detector created by Tom Lowell and posted on Loren's blog. Tom uses an IP camera on the ceiling and Image Processing Toolbox functions and the new webread function in R2014b and some File Exchange contributions to solve a difficult problem faced by MathWorkers in Natick, Massachusetts at the end of every week.

]]>

Today, I'd like to introduce a guest blogger, Tom Lowell, who is a program manager here at MathWorks. He works in our hardware connectivity group, and plays with Arduinos and Raspberry Pis in his spare time. He's relatively new to both MathWorks and MATLAB and decided to write a cookie detector as his first MATLAB project.

When I joined MathWorks, one of the first traditions I learned was that a large plate of home made cookies is delivered every Friday to every floor of every building here in our corporate headquarters in Natick. The problem is, they go fast! And you don't know when they'll be delivered. So if you're in a meeting or hunkered down at your desk, it's easy to miss out.

This seemed like a great opportunity to get some hands-on experience with MATLAB's image processing capabilities and solve an important problem. I'll describe the cookie detector and how it's built. And hopefully this will whet your appetite to go off and explore MATLAB's image processing capabilities on your own.

The cookie detector does one thing. It detects when a large circular plate of circular cookies is placed anywhere on the counter in the kitchen area on my floor. The system is very simple with only two components; an IP camera in the kitchen and MATLAB running on my computer.

IP cameras are similar to web cams except they run independently from any particular computer. They're attached to your network just like your laptop, your printer, or your smart phone. Most IP cameras are accessible via their IP address. Many, including mine, have a tiny web server running on their chipset. You just type their IP address into your web browser and you get back a JPG image or even video. The IP camera that I happen to be using is a Foscam FI8910W. But you can use any IP camera that supports URL addressing. Your MATLAB code will need to know the IP address of your IP camera. Therefore, it's desirable to configure your IP camera with a static IP address so it doesn't change over time.

You need a nearby outlet to power the IP camera, and a network connection. Most IP cameras support both WIFI or hardwired connections. We're using a WIFI connection because the computer that's running MATLAB is on the other side of the building. The ceiling in our kitchen happens to be metal. This made mounting the camera a breeze (using magnets).

All code for this project is contained within MATLAB. You will you also need the Image Processing Toolbox. The code is very simple. The main function loops forever, until a plate of cookies is detected. It downloads a new image from the camera every 30 seconds. MATLAB analyses the image using its built-in function `imfindcircles`. Once the cookies are detected, it sends a text message to several members on our team and then exits.

Getting the image from the camera into MATLAB was is the easiest part of the project. The camera has a tiny, built-in web server running on it that will return an image when requested via a URL. And MATLAB has several functions that will fetch data from web servers. I used `webread` which is new in R2014b. It reads content from a RESTful web service and returns it as Internet media types such as JSON, XML, images, or text.

Create a string with the camera's IP address plus specific API commands to pass to the camera. The format of this URL is camera specific. - snapshot.cgi tells the camera that we want a single image This camera also requires user and pwd as RESTful name-value pairs.

IPcamURL = 'http://172.31.28.24:1026/snapshot.cgi'; user = 'cookies'; password = 'lovecookies'

Use a timeout so that network issues don't cause it to hang.

options = weboptions('Timeout', 10);

The IP camera's web interface returns content as a JPG image stream. webread converts that for you into an RGB image that is immediately usable by MATLAB functions.

RGB = webread(IPcamURL, 'user', user, 'pwd', password, options); imshow(RGB);

The analysis for cookies is also pretty straightforward. The code first checks the image for a large circle (the round cookie plate) and then checks the image for smaller circles (cookies) contained within the large circle. My code calls the Image Processing Toolbox function `imfindcircles`, passing it the image from the camera, along with a range of valid radii to check for. Any results are returned in a vector of circles.

Check for a plate. If one or more plates are found, check for cookies... `detectPlateImfind()` and `detectCookiesImfind()` are just customized wrappers around `imfindcircles`

[plateCenter, plateRadius] = detectPlateImfind(tempJpgFile, smallPlateRad, bigPlateRad, 'g', 1);

nCookies = 0; nPlates = size(plateRadius, 1); if nPlates > 0 nCookies = detectCookiesImfind(tempJpgFile, smallCookieRad, bigCookieRad, 'r', 0, plateCenter, plateRadius); end

The analysis takes about 1 second per image on my MacBook Air.

Writing the basic code was easy because MATLAB insulates the programmer from the complexities of the various image-processing algorithms. However, you can (and need to) tune your code by passing various parameters to these image processing algorithms. I'll briefly describe the type of tuning I did here. For more information, there's a great example of tuning `imfindcircles` in the MathWorks documentation center.

First, there are few things to notice about the real life image above. The plate is a perfect circle with well-defined edges, a known radius, and high contrast relative to the counter top. Also, the cookies are generally circular but they're a bit irregular. They may be different colors and sizes each week.

The eventual call to the built-in function `imfindcircles` takes advantage of all of these. Here's the call to find the plate. The range of acceptable plate radii is passed in (smallest, biggest). The "object polarity" of the silver plate relative to the dark blue counter is passed in as "Bright". And some other parameters are set that were determined by tuning with a positive test case.

Find the plate(s)

[cirCen, cirRad, cirMetric] = imfindcircles(grayimg,[smallest, biggest],... 'ObjectPolarity','Bright', 'Method', 'Twostage','Sensitivity',0.96);

Notice that the parameters to the second call to imfindcircles to find the cookies are different. The object polarity is "dark" because the cookies are dark relative to the silver plate. A different (default) detection method was used and an edge threshold was used.

Find any cookies (on the plates)

[cirCen, cirRad] = imfindcircles(grayimg,[smallest, biggest],... 'ObjectPolarity','dark','Sensitivity',0.96, 'EdgeThreshold',0.05);

There are trade-offs to be made when you adjust these parameters. For example, if you want to detect every cookie on the plate, you might have to relax the parameters so much that the algorithm also incorrectly identifies other artifacts in your image as circles (other shapes, arcs, flashes of light). But if you tighten up your parameters to filter out those incorrect detections, you may lose detection of overlapping cookies, irregularly shaped cookies, etc.

The same logic applies to the plate or anything else you might want to detect. We don't want to miss the cookie plate, but we also don't want the tuning to be so permissive that we accidentally send 20 people rushing to the kitchen when someone sets their lunch plate on the counter.

The take-way here is that there is no one set of parameters to `imfindcircles` or any other image processing algorithm that will guarantee perfect results. You need to tune the parameters based on the real world characteristics of your image (such a the shape and color of the object, the lighting, angle, and other factors). This tuning is iterative and can take a while.

The output is a display of the image with the plate circled in green and the cookies in red. This is very easy to produce with MATLAB's basic plotting functions because imfindcircles returns a vector of circles (including the center points on the image and radii).

Although this is a pretty simple image, there are a few things to notice. It correctly detected plates of two different sizes. However, there are some false detects and it didn't get the exact number of cookies right. Lastly, those red crosses seen outside of the plate boundaries are detections that the loosely tuned algorithm thought were circles. However, some post-processing code was able to filter those out by determining they were not on the plate. This code can detect the plates with very good accuracy because of their well-defined shape and contrast. The code does not do as well with the cookies, which are irregular and offer less contrast relative to their background.

Remember, the object here is not a perfect cookie detector. The goal is correct and timely notification of the cookies' arrival. And it's pretty good at that. These cookies were detected before they were unwrapped!

Once the code has determined that there is a plate of cookies, it sends a text message to interested team members. I'm not going to cover the text messaging in this blog. It works by sending an email through the unique mail gateway of the cellular carriers. For example, to send a text to a Verizon user, MATLAB can send an email to 'xxxxxxxxxx@vtext.com' (where 'xxxxxxxxxx' is the cell phone number). I used Ke Feng's send text message example that's posted on MATLAB Central. The cool thing is that it's all done in MATLAB.

Do you have a project that uses an IP camera? We'd love to hear how you might use MATLAB with your IP camera project at school, work, or home. Tell us about it here.

Get
the MATLAB code

Published with MATLAB® R2014b

After reviewing the state of affairs fifty years ago, I use classic finite difference methods, followed by extrapolation, to find the first eigenvalue of the region underlying the MathWorks logo.

My Ph. D. dissertation, submitted to the Stanford Math Department in 1965, was titled "Finite Difference Methods for Eigenvalues of Laplace's Operator". The thesis was not just about the L-shaped region, but the L was my primary example.

My thesis advisor was George Forsythe, who had investigated the L in some of his own research. At the time, many people were interested in methods for obtaining guaranteed upper and lower bounds for eigenvalues of differential operators. Forsythe had shown that certain kinds of finite difference methods could provide lower bounds for convex regions. But the L is not convex, so that made it a region to investigate.

Nobody knew the first eigenvalue very precisely. Forsythe and Wolfgang Wasow published a book in 1960, "Finite Difference Methods for Partial Differential Equations" (John Wiley & Sons), that features the L in several different sections. In section 24.3 Forsythe references some calculations he had done in 1954 that led him to estimate that

$$ \lambda_1 = 9.636 $$

He goes on to report an informal conjecture made in the late '50s from his colleague R. De Vogelaere that

$$ 9.63968 < \lambda_1 < 9.63972 $$

I repeat these numbers here to give you some idea of the accuracy we were working with back then. We will see in a later post that the right answer is just outside De Vogelaere's upper bound.

Here's a page from my thesis. The contour plot at the bottom of the page would eventually evolve into the MathWorks logo. The plot was made on a Calcomp plotter, which involved a computer controlled solenoid holding a ball-point or liquid ink pen moving across a paper rolled on a rotating drum. This was a very effective device and was used in computer centers worldwide for a number of years.

The calculations shown on this page are the first eigenvalue of the finite difference Laplacian for ten values of the mesh size $h$. The finest mesh, $h$ = 1/100, has 29,601 interior points. As I recall, the calculation of the first eigenvalue, which I did with a kind of inverse power method, took about half an hour on Stanford's campus main frame, an IBM 7090.

Here is a plot of these values I computed fifty years ago. I've added a black line at what we now know to be the limiting value, the first eigenvalue of the continuous differential operator. You can see that the finite difference values are converging quite slowly. Furthermore, they are converging from above. Forsythe was not going to get his lower bound for this nonconvex region.

thesis_graph

The show rate of convergence of the finite difference eigenvalue with decreasing mesh size is caused by the fact that the corresponding eigenfunction is singular at the reentrant corner in the region. The gradient becomes infinite at that corner. You can see the contour lines crowding together near the corner. If you tried to make a drum head or tambourine by stretching a membrane over an L-shaped frame, it would rip.

To be more precise, use polar coordinates, $r$ and $\theta$, centered at the reentrant corner. Note that to cover the interior of the region $\theta$ will run from $0$ to $\frac{3}{2}\pi$. It turns out that approaching the corner the first eigenfunction has the asymptotic behavior

$$ u(r,\theta) \approx r^\alpha \sin{\alpha \theta} $$

with $\alpha = \frac{2}{3}$. This implies that the $k$ -th derivatives of the eigenfunction have singularities that behave like

$$ \partial^k u \approx r^{\alpha-k} $$

When this fact is combined with the Taylor series analysis of the discretization error of the five-point finite difference operator we find that the first eigenvalue converges at the rate

$$ \lambda_{1,h} - \lambda_1 \approx O(h^{2 \alpha}) = O(h^\frac{4}{3}) $$

This is significantly slower than the $O(h^2)$ convergence rate obtained when the eigenfunction is not singular.

Now let's use my laptop and the sparse capabilities in MATLAB. The functions `numgrid`, `delsq`, `spy`, and `eigs` have been included in the `sparfun` directory since its introduction. We begin by generating a small L-shaped numbering scheme.

```
n = 10;
L = numgrid('L',n+1)
```

L = 0 0 0 0 0 0 0 0 0 0 0 0 1 5 9 13 17 21 30 39 48 0 0 2 6 10 14 18 22 31 40 49 0 0 3 7 11 15 19 23 32 41 50 0 0 4 8 12 16 20 24 33 42 51 0 0 0 0 0 0 0 25 34 43 52 0 0 0 0 0 0 0 26 35 44 53 0 0 0 0 0 0 0 27 36 45 54 0 0 0 0 0 0 0 28 37 46 55 0 0 0 0 0 0 0 29 38 47 56 0 0 0 0 0 0 0 0 0 0 0 0

To illustrate how this scheme works, superimpose it on a finite difference grid.

Lgrid(L)

The statement

A = delsq(L);

constructs matrix representation of the five-point discrete Laplacian for the grid `L`.

```
whos A
issym = isequal(A,A')
nnzA = nnz(A)
```

Name Size Bytes Class Attributes A 56x56 3156 double sparse issym = 1 nnzA = 244

So for this grid, the discrete Laplacian `A` is a 56-by-56 sparse double precision symmetric matrix with 244 nonzero elements. That's an average of

ratio = nnz(A)/size(A,1)

ratio = 4.3571

nonzero elements per row. Each interior grid point is connected to its four nearest neighbors. For example, point number 44 is connected to points 35, 43, 45, and 53. So the nonzero elements in column 44 are

A(:,44)

ans = (35,1) -1 (43,1) -1 (44,1) 4 (45,1) -1 (53,1) -1

There are 4's on the diagonal and -1's is positions on the off-diagonals corresponding to connections between neighbors. They can be seen in the `spy` plot which shows the location of the nonzeros.

spy(A)

Scale `A` by the square of the mesh size and then use `eigs` to find its five smallest eigenvalues.

```
h = 2/n
A = A/h^2;
eigvals = eigs(A,5,'sm')
```

h = 0.2000 eigvals = 30.3140 27.7177 19.0983 14.6708 9.6786

Run for decreasing mesh size, up to the memory capacity of my laptop.

```
type L_finite_diffs
```

% L_finite_diffs % Script for MathWorks Logo, Part Two % Finite Differences opts.tol = 1.0e-12; eigvals = zeros(13,5); tic for n = 50:50:650 L = numgrid('L',n+1); h = 2/n; A = delsq(L)/h^2; e = eigs(A,5,'sm',opts); eigvals(n/50,:) = fliplr(e'); fprintf('%4.0f %10.6f %9.0f %10.0f %16.12f\n', ... n,h,size(A,1),nnz(A),e(5)) end toc save eigvals eigvals

```
type L_finite_diffs_results
```

n h size(A) nnz(A) lambda_1,h 50 0.040000 1776 8684 9.661331286788 100 0.020000 7301 36109 9.649547111146 150 0.013333 16576 82284 9.645736333124 200 0.010000 29601 147209 9.643932372382 250 0.008000 46376 230884 9.642903412412 300 0.006667 66901 333309 9.642247498956 350 0.005714 91176 454484 9.641797238078 400 0.005000 119201 594409 9.641471331746 450 0.004444 150976 753084 9.641225852818 500 0.004000 186501 930509 9.641035124947 550 0.003636 225776 1126684 9.640883203502 600 0.003333 268801 1341609 9.640759699387 650 0.003077 315576 1575284 9.640657572983 Elapsed time is 29.729023 seconds.

I know the Taylor series analysis of the discretization error involves terms with both $h^{4/3}$ and $h^2$, so let's use these as the basis for extrapolation. Use backslash to do a least squares fit.

load eigvals eigvals = eigvals(:,1); n = (50:50:650)'; h = 2./n; e = ones(size(h)); X = [e h.^(4/3) h.^2]; c = X(4:end,:)\eigvals(4:end); lambda1 = c(1); n = (50:10:650)'; h = 2./n; fit = c(1) + c(2)*h.^(4/3) + c(3)*h.^2; plot(n(1:5:end),eigvals,'o',n,fit,'-',[0 700],[lambda1 lambda1],'k-') set(gca,'xtick',100:100:600,'xticklabel',... num2str(sprintf('2/%3.0f\n',100:100:600))) xlabel('h')

Frankly, I was surprised when the extrapolation of the finite difference eigenvalues produced this result.

```
format long
lambda1
```

lambda1 = 9.639723753239963

I know from techniques that I will describe in later posts that this agrees with the exact eigenvalue of the continuous differential operator to seven decimal places. Before I wrote this blog I never tried this computation on computers with enough memory to handle meshes this fine and was never able to obtain this kind of accuracy from extrapolated finite difference eigenvalues.

Let's do a contour plot of the eigenfunction. The tricky part is using the grid numbering in `L` to insert the eigenvector back onto the grid.

close n = 100; h = 2/n; L = numgrid('L',n+1); A = delsq(L)/h^2; [V,E] = eigs(A,5,'sm'); L(L>0) = -V(:,5); contourf(fliplr(L),12); axis square axis off

The plot shows off `parula`, the MATLAB default colormap with Handle Graphics II and Release 2014b. The colormap derives its name from the parula, which is a small warbler found throughout eastern North America that exhibits these colors. As our contour plot demonstrates, the color intensity varies smooth and uniformly as the contour levels increase. This is not true of `jet`, the old default. Steve Eddins has more about `parula` in his blog.

Get
the MATLAB code

Published with MATLAB® R2014b

Last week I showed you the new MATLAB colormap, parula. Parula has replaced jet as the default colormap in R2014b, which was released earlier this month.

This week I want to explain some of the motivations for replacing jet. Jet is an example of a *rainbow colormap*. A rainbow colormap is based on the order of colors in the spectrum of visible light. Some of you might know the mnemonic "Roy G. Biv" for remembering the order of colors: red, orange, yellow, green, blue, indigo, and violet.

Multiple methods exists for constructing rainbow colormaps, so they don't all look alike. Here is jet, the MATLAB version:

showColormap('jet','bar')

In my introductory post last week, I showed you a few visualization examples using jet, and I asked you some questions about them. Let's look at the questions again and figure out the answers.

**Question 1:** In the chart below, as you move from left to right along the line shown in yellow, how does the data change? Does it trend higher? Or lower? Or does it trend higher in some places and lower in others?

fluidJet hold on plot([100 160],[100 135],'y','LineWidth',3) hold off colormap(gca,'jet')

By my count there are six very distinct color changes as you move along the line from left to right. So what does that mean about the data?

I posed this question to someone who is very familiar with the jet colormap, and I was surprised by the response: "The data starts low, goes high, and then goes back low again."

"But what about all those other color stripes," I asked? "What do they mean?"

"That's just what jet does. I've learned to ignore all those stripes."

Oh! Well, it turns out that some of the color stripes are indeed mostly meaningless. But not all of them!

Here's what the data actually does along that path.

improfile(s.X,[100 160],[100 135])

Yes, the data starts low, goes high, and then goes back low. But if you ignore all the color stripes in between the blue and the red, then you miss the presence of the smaller-magnitude peak on the left. Alternatively, if you believe that all the color stripes are meaningful, then the yellow stripes seem to falsely suggest distinct data regions to the left and the right of the main peak.

Now let's compare with parula.

fluidJet colormap(gca,parula) hold on plot([100 160],[100 135],'y','LineWidth',3) hold off

You can see that now there are only four distinct color changes as you follow the yellow line. These changes correspond directly to real features in the data. There are no extraneous color stripes that do not correspond to real data features.

**Question 2:** In the filled contour plot below, which regions are high and which regions are low?

filledContour15 colormap(gca,rgb2gray(jet))

In the plot above, I have done a crude simulation of what a jet-based visualization looks like if printed on a grayscale printer. Here's the full color version:

colormap(gca,jet)

**If** you are familiar with jet and know therefore that blue is low and red is high, then you can probably look at the full-color version and give a reasonable interpretation of it. When rendered as grayscale, however, all hope is lost. The red and blue colors of jet are just about equally dark.

Here are the full-color and grayscale results with parula.

filledContour15 colormap(gca,parula)

filledContour15 colormap(gca,rgb2gray(parula))

With parula, dark consistently means low, and bright consistently means high.

The last questions relate to the three plots below (A, B, and C) showing different horizontal oscillations.

**Question 3:** Which horizontal oscillation (A, B, or C) has the highest amplitude?

**Question 4:** Which horizontal oscillation (A, B, or C) is closest to a pure sinusoid?

**Question 5:** In comparing plots A and C, which one starts high and goes low, and which one starts low and goes high?

subplot(2,2,1) horizontalOscillation1 title('A') subplot(2,2,2) horizontalOscillation2 title('B') subplot(2,2,3) horizontalOscillation3 title('C') colormap(jet)

The answers are subjective and might depend on your monitor and ambient lighting. But here's what I see.

The data pattern in image A has significantly higher visual contrast than either image B or image C. That suggests that the data oscillation in A has the highest amplitude.

The data patterns in images B and C both exhibit vertical stripes within each period of the oscillation. That strongly suggests that horizontal oscillations shown in B and C have some constant or almost-constant regions and that they go up and down in something like a stair-step pattern. To my eye, image A looks the most like a pure sinusoid.

And finally, the pattern in image A starts with a bright stripe on the left, whereas the image B pattern starts with a dark stripe. That suggests that oscillation A starts high, while oscillation B starts low (or maybe vice versa).

You probably won't be surprised to find out that these are all trick questions. Actually, oscillations A, B, and C are all sinusoids with exactly the same amplitude and shape. They differ only in their constant offsets, which place them in different regions of the jet colormap. Note especially that oscillations A and C are actually in phase with each other; they both start high and go low. The differences we see are all visual artifacts caused by the jet colormap.

Compare the result using parula:

colormap(parula)

With these examples, I wanted to show several ways in which the jet colormap, and rainbow colormaps in general, can mislead the viewer about properties of the data.

It turns out that some people have been writing about these problems with rainbow colormaps for years. Here's a summary of the main criticisms:

- Rainbow colormaps confuse viewers because there is no natural perceptual ordering of the spectral colors. In addition to causing visual confusion (such as whether oscillations A and C above are in-phase or out-of-phase), the lack of perceptual ordering can slow down tasks because viewers have to refer to the color key more often in order to interpret the data.

- Rainbow colormaps obscure small details in the data. The primary reason is that the green and especially the cyan sections of the rainbow colormap are perceptually indistinct, which makes the data in the corresponding ranges appear to be uniform or flat. (I haven't shown examples of this in the blog, yet, but there are examples in the technical paper mentioned below.)

- Rainbow colormaps mislead viewers by suggesting data features that are not really there. These "phantom features" often take the form of false boundaries. This effect, in combination with perceptually indistinct green or cyan regions, can falsely segment the data.

- Rainbow colormaps lose critical information about high and low data values when displayed or printed on a gray-scale device.

- Rainbow colormaps can be difficult to interpret for some color-impaired viewers.

For a much more detailed summary of the literature and online sources regarding rainbow colormaps, see the paper "Rainbow Color Map Critiques: An Overview and Annotated Bibliography" on mathworks.com.

Next time I'll starting getting into colormap construction principles.

Get
the MATLAB code

Published with MATLAB® R2014b

by Toshi Takeuchi

One of the new R2014b features that deserves your attention is Git integration. Git is a source control system (also known as version control or source code management system) that enables collaborative software development. Why does that matter to you? Programming is an essential skill in many technical fields even outside computer science, and some universities now offer software carpentry workshops to enhance coding skills for researchers. Source control is one of those essential skills in software carpentry.

Until now, you may have tinkered alone with code you needed for your project. However, there are other people who may be working on similar problems and they may be writing similar programs. Source control enables you to work with other people so that you don’t have to do it all alone. Collaboration lets you be more productive in other aspects of your project.

Even if you don’t care about such collaboration, wouldn’t it be cool to share your personal project and see other people using it? They may even fix bugs and improve your code for you!

GitHub is one of the most popular websites that host Git repositories. The best place to share your MATLAB projects is File Exchange because of its popularity with the MATLAB user community. And guess what – File Exchange is integrated with GitHub! Now you see the connection?

What is a Git repository? A repo (repository) is a directory that holds your source code and any associated files. Local repos are on your local drive, and the remote repos are on GitHub or other hosts, and you sync the local repos to remote repos as you write your code. You can start with either the local or remote repos, but in this example I am going to start with a remote repo.

The process looks like this for a single developer:

- Create or fork a repo on GitHub
- Clone the repo to your local drive – this is your local repo
- Add your files to the local repo
- Sync your local repo to remote repo
- Repeat this process to keep your source code in sync as you write more code
- Share your GitHub repo on File Exchange

When you talk about Git, you cannot go without mentioning “forking”. In the simplest terms forking means copying someone else’s public repo on the remote server, rather than starting a repo from scratch. In practice forking is used as a way to contribute to the existing projects or to start a new project using an existing project as a starting point. Once you make changes to your forked project, you can send a merge request to the original developer, and your changes may be accepted and incorporated into the main project.

Forking enables a flexible distributed style of collaboration and number of forks you have on your project acts as a measure of popularity – similar to the count of likes or followers on Facebook or Twitter. The social aspect of forking is an interesting topic on its own, but we need to skip it for this post.

Signing up with GitHub is easy – just click the *sign up* button on the homepage and follow instructions. A free account would do for now. You also need to download and install Git. Even though GitHub has GUI apps for Windows and Mac, you need to set up the command-line tool for use with MATLAB. You also want to follow these instructions for registering binary files with Git.

Creating a repo on GitHub is very easy – just follow these instructions. From this point on I will assume you named your repo *Hello-World* and initialized it with a README file. Please note that you can only create a public repo with a free account.

Until recently, you needed to use the command line tool for this step, but starting with R2014b we can just use MATLAB’s Current Folder window. No more Git commands like `git init`, `git status`, `git add`, or `git commit`!

Open your copy of MATLAB and create an empty folder. Right-clicking the empty space in the Current Folder window to bring up a contextual menu, and select *Source Control > Manage Files*.

This will open a new dialog box: *Manage Files Using Source Control*.

- For
*Select control integration*, choose*Git* - For
*Repository path*, click*Change*

You now see a new dialog box: *Select a Repository*. Copy and paste the URL of the remote repo you just created. You can find the URL in the right sidebar of your new repo on GitHub.

You choose either SSH or HTTPS depending on how you setup your authentication on GitHub.

Click *Validate*. You may be asked for your login password for authentication. You can close the dialog box when your path is validated.

Back in *Manage Files* dialog box, the *sandbox* should be already set to your current folder. All you need to do now is hit *Retrieve*.

You have now successfully cloned the remote repo to your local drive. Check your Current Folder window. You should see just one file – README.md, but with a green circle next to it. This is just a text file but you can apply wiki-like syntax called markdown to make it appear like a regular web page when viewed on GitHub. README serves as the front page of your repo on GitHub.

Let’s add a new MATLAB script file *helloworld.m*. It will appear with a blank circle – it means it is not added to Git source control yet. To add it to Git, right-click on the file and select *Add to Git*. The empty circle changes to “+” symbol. When you modify a file already under source control, the symbol becomes a blue square.

You can continue editing files as you like, but at some point, you want to take a snapshot of the edits you made. That’s when you do a commit. You can select any empty space in the Current Folder window to bring up the contextual menu and select *Commit to Git Repository*. This will bring up a dialog box where you can add your comment about the changes you made since the last commit. Comments will be helpful to keep to track of your changes and revert back to earlier commits if necessary.

When you commit, the snapshot is saved in the local repo, but it is also a good idea to mirror the changes to the remote repo as well. To do so, bring up the contextual menu by right-clicking an empty space in the Current Folder window and select *Push*. That will push your changes to the remote repo. You may need to enter your password.

The real power of source control comes from the ability to create multiple branches from your project. By default, you have a branch called “master” in your repo. You can create a new branch from the master branch, makes changes, and then merge those changes back to the master. This mechanism is used for working on new experimental features without affecting the working code on the master. You can branch and merge in MATLAB but the details are beyond the scope of this post.

If you have been curious about Git but put off by its complicated command line interface, Git integration in R2014b makes Git much more pleasant and approachable. I hope this quick introduction motivates you to take advantage of this new feature. When you do, please don’t forget to post your project to the File Exchange. To learn more about Git, it is actually helpful to start with reviewing the underlying concept about how Git works.

]]>Jiro's pick this week is MATLAB Plot Gallery by Plot Gallery.

This gallery, which hosts a collection of File Exchange entries showing various plotting examples, has been highlighted once before. I want to bring this up once again, since the examples are now updated with the new graphics system in R2014b. The example code is backward compatible, so you can still run them in older releases, but running them in the current release (R2014b) will give you the new look.

So what's changed with the new graphics system? ** A lot!!** Rather than, talking about them here, let me point you to a few great resources. The first is a series of blog posts in Loren's blog:

These posts will give you a nice introduction to the new graphics system.

Second, there are a couple of blog posts about the new colormap in Steve's blog:

And finally, there is the newly introduced blog, Mike on MATLAB Graphics. Mike is a MathWorks developer working on MATLAB's graphics system, and he will be writing about computer graphics concepts and MATLAB graphics techniques.

**Comments**

Take a look at these resources and let us know what you think here.

Get
the MATLAB code

Published with MATLAB® R2014b

Coursera is a technology platform that kickstarted the current MOOCs boom. Even though there are more MOOCs players now, it still remains one of the leading companies in this space. But how are they doing these days for delivering higher education to the masses online?

Today's guest blogger, Toshi Takeuchi, would like to share an analysis using Courera's data.

I am a big fan of MOOCs and I benefited a lot from free online courses on Coursera, such as Stanford's Machine Learning course. Like many websites these days, Coursera offers its data through REST APIs. Coursera offers a number of APIs, but Catalog APIs are available without OAuth authentication. We can find out the details of courses offered by Coursera with these APIs.

We can try to answer questions like *"how do STEM and non-STEM courses break down among universities?*"

JSON is a very common data format for REST APIs, and Coursera's APIs also returns results in JSON format. MATLAB now supports JSON out of the box in R2014b. You could always use JSON from within MATLAB by taking advantage of user contributed MATLAB programs on File Exchange, but built-in JSON support makes it easy for us to share scripts that use JSON, because we don't have to worry about dependencies.

Let's try the new feature using Coursera APIs. Calling a REST API is very simple with `webread`.

restApi='https://api.coursera.org/api/catalog.v1/courses'; params = 'sessions,universities,categories'; resp=webread(restApi,'includes',params,weboptions('Timeout',60));

`webread` returns the JSON response as a structure array. The data is further processed in a separate script `processData.m` - check out the details if interested.

We need to decide which categories represent STEM subjects. When there are multiple categories assigned to a given course, we treat it as a STEM course as long as one of them is included in STEM categories.

processData

STEM categories 'Computer Science: Theory' 'Economics & Finance' 'Medicine' 'Mathematics' 'Physical & Earth Sciences' 'Biology & Life Sciences' 'Computer Science: Systems & Security' 'Computer Science: Software Engineering' 'Engineering' 'Statistics and Data Analysis' 'Computer Science: Artificial Intelligence' 'Physics' 'Chemistry' 'Energy & Earth Sciences'

As a sanity check, let's plot the number of courses vs. number of sessions by university. A single course can be offered repeatedly in multiple sessions. Therefore you can determine the longevity or age of a given course by the count of sessions.

If it is a new course, or it was not repeated, then you only have one session per course. We can use this as the baseline, and check how universities scaled up their courses relative to this baseline.

R2014b comes with new MATLAB Graphics System, but you can still use the familiar commands for plotting.

% group by number of courses grouping = ones(height(universities),1)*2; grouping(universities.courses > 25) = 1; grouping(universities.courses <= 10) = 3; % plot figure gscatter(universities.courses,universities.sessions,grouping) h = refline(1,0); set(h,'Color','m','LineStyle',':') h = refline(2,0); set(h,'Color','m','LineStyle',':') h = refline(3,0); set(h,'Color','m','LineStyle',':') h = refline(6,0); set(h,'Color','m','LineStyle',':') xlabel('Number of Courses'); ylabel('Number of Sessions'); title('\fontsize{14} Courses by Sessions by University'); legend('Universities with 25+ courses','Universities with 10+ courses',... 'Universities with 1-10 courses','Ref line: 1 session per course',... 'Ref line: 2 sessions per course','Ref line: 3 sessions per course',... 'Ref line: 6 sessions per course','Location','NorthWest') % add university names for i = 1:height(universities) if universities.courses(i) > 10 && universities.sessions(i) > 20 text(universities.courses(i),universities.sessions(i),... universities.shortName{i},'FontSize',12) end end

You can see that Stanford, Penn (University of Pennsylvania), JHU (Johns Hopkins), and Duke are leading the pack. They are the early adopters, based on the number of sessions. It is interesting to see PKU (Peking University) leading international institutions. They offer a number of courses in Chinese. Coursera didn't start international partnership until recently, so it is quite remarkable the PKU has broadened their online content in relatively short time. More recent entrants are on the left with fewer courses and sessions.

Established players are trying to scale up by repeating the sessions. JHU seems to be particularly aggressive in terms of the number of courses they offer and how they are repeated as sessions.

Let's plot the number of courses by ratio of STEM courses by university. This will tell us which schools are making investments in online education content, and whether they focus on STEM or non-STEM subjects. The size of the marker indicates the total number of sessions they are associated with, so it also gives us how long they have been involved in Coursera. Notice the `parula` colormap used in the colorbar, the new default colormap in R2014b.

% use sesion count for setting marker sizes markerSize = universities.sessions; % we need to scale the marker size markerSize = (markerSize - min(markerSize))/(max(markerSize)-min(markerSize)); markerSize = markerSize * 1000; markerSize(markerSize == 0) = 1; % change the tick labels to reflect the original values barticks = num2cell(20:20:200); % create a scatter plot figure scatter(universities.courses,universities.stem_ratio,markerSize,markerSize,'fill') xlim([0 40]) h = colorbar('TickLabels',barticks); h.Label.String = '\fontsize{11}Number of Sessions'; title('\fontsize{14} Ratio of STEM courses by University on Coursera') xlabel('\fontsize{11}Number of Courses'); ylabel('\fontsize{11}Ratio of STEM Courses'); % add university names for i = 1:height(universities) if universities.stem_ratio(i) ~= 0 && universities.stem_ratio(i) ~= 1 && universities.courses(i) >= 5 text(universities.courses(i),universities.stem_ratio(i),universities.shortName{i},'FontSize',12) end end % add reference lines line([25 25],[0 1],'LineStyle',':') line([10 10],[0 1],'LineStyle',':') line([0 40],[0.5 0.5],'LineStyle',':')

Stanford is very heavy on STEM subjects, while others are more balanced. More recent entrants on the left have a wider variance in how STEM heavy their courses are. Perhaps rate of adaption is different among different academic disciplines?

We can plot the ratio of courses per category in order to see the relative representation of academic disciplines on Coursera. A course can belong to multiple categories, and in such cases a count is split equally across the included categories. Note that you can now rotate axis tick labels in R2014b.

% get the count of categories by university catByUniv = zeros(height(universities),height(categories)); for i = 1:length(T.categories) row = ismember(universities.id,T.universities(i)); col = ismember(categories.id,T.categories{i}); catByUniv(row,col) = catByUniv(row,col) + 1/length(T.categories{i}); end % segment the universities by number of courses catByTiers = [sum(catByUniv(grouping == 1,:));... sum(catByUniv(grouping == 2,:)); sum(catByUniv(grouping == 3,:))]; % get the ranking of categories by number of courses [~,ranking] = sort(sum(catByUniv(universities.courses > 25,:)),'descend'); % get the ratio of courses by category catByTiers = bsxfun(@rdivide,catByTiers,sum(catByTiers,2)); % plot a bar graph figure xticks = [{''};categories.name(ranking);{''}]; h = bar(catByTiers(:,ranking)'); xlim([0 26]); ax = gca; set(ax,'XTick',0:26);set(ax,'XTickLabel',xticks);set(ax,'XTickLabelRotation',270); title('\fontsize{14} Ratio of Courses Per Category') legend('Universites with 25+ courses','Universites with 10+ courses','Universites with 1-10 courses','Location','Best')

It looks like there was more STEM bias among the early adopters (universities with a lot of courses) but new entrants (universities with fewer courses) tend to have more non-STEM courses. Categories like Social Sciences, Humanities, Business and Management, Education, Teacher Professioal Development, Music, Film and Audio are on the rise.

Why do we see this non-STEM shift? There are a number of possible explanations.

- In the beginning, Coursera courses relied on autograders. They were well suited for quantitative STEM subjects, but not for non-STEM subjects.
- Autograders were custom built for respective courses and they are in fact full SaaS applications. It was difficult to scale the number of courses if you needed to build for each course a custom SaaS app that can withstand substantial peak traffic near the deadline -
*this human behavior is pretty universal* - Later, Coursera introduced a crowd sourced essay grading system that can be used across multiple courses. This freed universities from the burden of creating custom SaaS apps.
- This led to rapid expansion of course offerings and made non-STEM subjects viable. In fact, I took a number of STEM courses from JHU, and they tend to use essay grading system rather than autograders.

There are questions we cannot answer with the data at hand. For example, is the shift driven by the convenience of supply side (universities) or by the demand for non-STEM subjects by the public?

There are no strict prerequisites for Coursera courses, but the bar is still high for STEM courses. Therefore it is quite possible that the potential market size is larger for non-STEM subjects.

You also saw how easy it is to use REST API with JSON response within R2014b, and got a quick look at some of the new features of updated MATLAB Graphics System. Download the new release, try those new features yourself and share what you find here!

Get
the MATLAB code

Published with MATLAB® R2014b

**Differential Equation Editor**

Try typing `dee` in MATLAB. This will launch an example model that looks like:

If you open one of the demo and double-click on the block, you will see a nice little user interface:

In this interface, you can type any equation you want, using the format of the Fcn block.

When I first saw this model and user interface, I thought *"This thing looks really old, has it always been there, unnoticed?"*. So I looked into our source code history and found out this example ships with Simulink and has been untouched since the middle of the 90's!

**How does that work?**

I strongly encourage you to look under the mask and observe how this block is implemented. It is a pretty good example of how to create models programmatically. Once you click the rebuild button in the user interface, functions like `delete_block`, `add_block` and `add_line` draw the differential equation in a systematic manner.

For any set of equation, the resulting model will look like this, where one Fcn block is used for each state equation and for each output equation:

To help you visualize it better, I formatted the above model to make it easier to read:

**Now it's your turn**

Give this cute little example a try and let us know what you think by leaving a comment here.

]]>Today, David Garrison, our guest blogger, will continue his series on the new graphics system in R2014b.

- Part 1: Features of the New Graphics System
**Part 2: Using Graphics Objects**- Part 3: Compatibility Considerations in the New Graphics System

Here is Part 2 of the series.

In Part 1 of this series, I provided an introduction to the new MATLAB graphics system in R2014b. I described a number of new features and alluded to one of the big changes in R2014b -- *graphics functions return MATLAB objects, not numeric handles*.

When I use MATLAB, I usually don't think about the internal graphics system. Usually, I'm just trying to visualize some data. I might call one of the MATLAB charting functions, like `plot` or `bar` to understand something important about my data. I also create user interfaces using functions like `uicontrol` or `uipanel`. These are functions that are part of the MATLAB graphics system.

Let's suppose that I am creating a plot and I want the line in the plot to be red. I can do that in one of two ways. I can call the `plot` function with some extra arguments like this.

x = 0:0.1:10; y = sin(x); plot(x,y,'Color','red')

I can also use the output argument of the `plot` command to modify my plot after I've created it. Here I created the plot first, then used the `set` command to change the color of the line.

p = plot(x,y) ; set(p, 'Color', 'red')

Prior to R2014b, a call to a graphics function would return a number like this

p = plot(x,y)

p = 174.346

The value of the variable `p` was a special kind of number called a *handle*. The handle was a reference to a graphics object. The object itself was not available to the user but you could use the handle to retrieve or change properties of the object using the `set` and `get` commands. In the example above, we used the handle of the line object to set its color. In fact, when this capability first became available in MATLAB, this approach was so new and useful it was given a special name -- *Handle Graphics*.

All graphics handles were just numeric values. Those values were integers for figures and fractional values for everything else. You could use a handle anywhere you could use a number. I've seen MATLAB code where people stored their handles with other data in a double array or used handles in functions where a number is expected as the argument (e.g. math functions).

As you probably know, graphics objects are stored as a tree structure with parents and children. There is a special object at the top of the tree called the *graphics root*. Prior to R2014b, the handle to the graphics root was always `0`. For example to get the list of all open figures, you could type:

```
myFigures = get(0, 'Children');
```

In R2014b, graphics functions do not return numeric handles any more. They now return MATLAB objects. Now when you call a graphics function with an output argument, you see something like this:

p = plot(x,y)

p = Line with properties: Color: [0 0.4470 0.7410] LineStyle: '-' LineWidth: 0.5000 Marker: 'none' MarkerSize: 6 MarkerFaceColor: 'none' XData: [1x101 double] YData: [1x101 double] ZData: [1x0 double] Use GET to show all properties

The `plot` function now returns a `Line` object. That line object has properties. When MATLAB displays the line object in the command window, it shows the object's most common properties. The `all properties` text is a hyperlink that you can click to see all the properties of the object.

If you use the `whos` function to see more information about `p`, you see something like this:

```
whos p
```

Name Size Bytes Class Attributes p 1x1 112 matlab.graphics.chart.primitive.Line

The *Class* column in the `whos` output indicates what kind of object if is. What you see here is the full class name of the object. Don't worry too much about the details of the class name. Most of you will never need to use it.

In R2014b, use the `groot` function to refer to the graphics root.

groot

ans = Graphics Root with properties: CurrentFigure: [1x1 Figure] ScreenPixelsPerInch: 96 ScreenSize: [1 1 1920 1200] MonitorPositions: [1 1 1920 1200] Units: 'pixels' Use GET to show all properties

One of the benefits of having graphic objects is that you can now get and set their values using dot notation. Dot notation has the form `object.Property`. It is the same syntax that you use when you want to refer to the field of a structure. For example, in R2014b, you can get the color of the line above like this:

p.Color

ans = 0 0.4470 0.7410

Similarly, you can set the line width of the line using dot notation.

p.LineWidth = 2;

Dot notation is also useful if you want to set properties of two objects to be the same. Prior to R2014b, you might write code that looks like this:

set(p2, 'Color', get(p1, 'Color'))

With dot notation, you can do the same thing with

p2.Color = p1.Color;

With dot notation you can also use tab completion when you are trying to access a graphics object property.

You can still use the `set` and `get` functions but I find dot notation a much more convenient way to refer to object properties. One note of caution, however. When you use dot notation you must use the correct capitalization of the property name. For example,

p.LineWidth

will return the value of the line width but

p.linewidth

will cause an error.

As I mentioned above, you can still use the `set` and `get` functions and there are cases where `set` and `get` are useful.

You can use the `get` command if you want to know all the properties of a graphics object. This is the same list of properties you see if you click the `all properties` link in the example above.

get(p)

AlignVertexCenters: 'off' Annotation: [1x1 matlab.graphics.eventdata.Annotation] BeingDeleted: 'off' BusyAction: 'queue' ButtonDownFcn: '' Children: [] Clipping: 'on' Color: [0 0.4470 0.7410] CreateFcn: '' DeleteFcn: '' DisplayName: '' HandleVisibility: 'on' HitTest: 'on' Interruptible: 'on' LineStyle: '-' LineWidth: 2 Marker: 'none' MarkerEdgeColor: 'auto' MarkerFaceColor: 'none' MarkerSize: 6 Parent: [1x1 Axes] PickableParts: 'visible' Selected: 'off' SelectionHighlight: 'on' Tag: '' Type: 'line' UIContextMenu: [] UserData: [] Visible: 'on' XData: [1x101 double] XDataMode: 'manual' XDataSource: '' YData: [1x101 double] YDataSource: '' ZData: [1x0 double] ZDataSource: ''

The `set` function is useful if you want to know what options are available for a given property. In this case you can call `set` with the property name but no value for the property.

```
set(p,'LineStyle')
```

'-' '--' ':' '-.' 'none'

Another use of `set` occurs when you need to reference a property from an array of graphics objects. Suppose you have an array of line objects created by the `plot` command:

x = 0:0.1:10; y1 = sin(x); y2 = cos(x); lineArray = plot(x,y1,x,y2)

lineArray = 2x1 Line array: Line Line

You can use the `set` function to change the colors of all the lines in the array in a single command.

set(lineArray, 'Color', 'blue')

What do you think about the change to graphics objects in R2014b? Have you tried using dot notation for getting and setting graphics object properties? We'd love to hear your thoughts here.

There are some other changes in R2014b that will impact some existing MATLAB code. These changes mostly impact advanced graphics users and people who build more complicated user interfaces. Next time, I will describe the important changes and what steps you should take to make your code compatible with R2014b. I'll also give you some guidance on how to write code that works in multiple releases.

Get
the MATLAB code

Published with MATLAB® R2014b

I believe it was almost four years ago that we started kicking around the idea of changing the default colormap in MATLAB. Now, with the major update of the MATLAB graphics system in R2014b, the colormap change has finally happened. Today I'd like to introduce you to parula, the new default MATLAB colormap:

showColormap('parula','bar')

(Note: I will make `showColormap` and other functions used below available on MATLAB Central File Exchange soon.)

I'm going spend the next several weeks writing about this change. In particular, I plan to discuss why we made this change and how we settled on parula as the new default.

To get started, I want to show you some visualizations using jet, the previous colormap, and ask you some questions about them.

**Question 1:** In the chart below, as you move from left to right along the line shown in yellow, how does the data change? Does it trend higher? Or lower? Or does it trend higher in some places and lower in others?

fluidJet hold on plot([100 160],[100 135],'y','LineWidth',3) hold off colormap(gca,'jet')

**Question 2:** In the filled contour plot below, which regions are high and which regions are low?

filledContour15 colormap(gca,rgb2gray(jet))

The next questions relate to the three plots below (A, B, and C) showing different horizontal oscillations.

**Question 3:** Which horizontal oscillation (A, B, or C) has the highest amplitude?

**Question 4:** Which horizontal oscillation (A, B, or C) is closest to a pure sinusoid?

**Question 5:** In comparing plots A and C, which one starts high and goes low, and which one starts low and goes high?

subplot(2,2,1) horizontalOscillation1 title('A') subplot(2,2,2) horizontalOscillation2 title('B') subplot(2,2,3) horizontalOscillation3 title('C') colormap(jet)

Next time I'll answer the questions above as a way to launch into consideration of the strengths and weaknesses of jet, the previous default colormap. Then, over the next few weeks, I'll explore issues in using color for data visualization, colormap construction principles, use of the L*a*b* color space, and quantitative and qualitative comparisons of parula and jet. Toward the end, I'll even discuss how the unusual name for the new colormap came about. I've created a new blog category (colormap) to gather the posts together in a series.

If you want a preview of some of the issues I'll be discussing, take a look at the technical paper "Rainbow Color Map Critiques: An Overview and Annotated Bibliography." It was published on mathworks.com just last week.

Get
the MATLAB code

Published with MATLAB® R2014b

MathWorks is the only company in the world whose logo satisfies a partial differential equation. Why is the region for this equation shaped like a capital letter L?

The wave equation describes how a disturbance travels through matter. If the units are chosen so that the wave propagation speed is equal to one, the amplitude of a wave satisfies

$$ {{\partial^2 u} \over {\partial t^2}} = \triangle u $$

The $\triangle$ denotes Laplace's operator

$$ \triangle = {\partial^2 \over {\partial x^2}} + {\partial^2 \over {\partial y^2}} $$

Geometry plays a crucial role here. Initial values of the amplitude and velocity of the wave are prescribed on a certain region. Values of the amplitude or its normal derivative are also prescribed on the boundary of the region. If the wave vanishes outside the region, these boundary values are zero.

Separating out periodic time behavior leads to solutions of the form

$$ u(x,y,t) = \cos{(\sqrt{\lambda}\,t)} v(x,y) $$

The functions $v(x,y)$ also depend upon $\lambda$ and the region. They satisfy the differential equation

$$ \triangle v + \lambda v = 0 $$

and are zero on the boundary of the region. The quantities $\lambda$ that lead to nonzero solutions are the *eigenvalues*, and the corresponding functions $v(x,y)$ are the *eigenfunctions* or *modes*. They are determined by the physical properties of the medium and the geometry of the region. The square roots of the eigenvalues are resonant frequencies. A periodic external driving force at one of these frequencies generates an unboundedly strong response in the medium.

Any solution of the wave equation can be expressed as a linear combination of these eigenfunctions. The coefficients in the linear combination are obtained from the initial conditions.

In one dimension, the eigenvalues and eigenfunctions are easily determined. The simplest example is a violin string, held fixed at the ends of an interval. The eigenfunctions are trig functions.

$$ v_k(x) = \sin{(k x)} $$

If the length of the interval is $\pi$, the eigenvalues are determined by the boundary condition, $v_k(k \pi) = 0$. Hence, $k$ must be an integer and

$$ \lambda_k = k^2 $$

If the initial condition is expanded in a Fourier sine series,

$$ u(x,0) = \sum_k a_k \sin{(k x)} $$

(And the initial velocity is zero), then the solution to the wave equation is

$$ u(x,t) = \sum_k a_k \cos{(\sqrt{\lambda_k}\,t)} v_k(x) $$

Here are graphs of the first nine eigenfunctions in one dimension. The corresponding eigenvalues are the squares of integers.

```
eigenvals = (1:9).^2
plot_modes('1d')
```

eigenvals = 1 4 9 16 25 36 49 64 81

The simplest region in two dimensions is a square. The eigenfunctions are again trig functions.

$$ v_{k,j}(x,y) = \sin{(k x)}\,\sin{(j y)} $$

If the sides have length $\pi$, the boundary conditions imply that $k$ and $j$ must be integers. Here are the first nine eigenvalues and eigenfunctions.

```
[k,j] = meshgrid(1:3); e = k.^2+j.^2; eigenvals = e(:)'
plot_modes('square')
```

eigenvals = 2 5 10 5 8 13 10 13 18

If the region is a circular disc, we switch to polar coordinates, $r$ and $\theta$. Trig functions are replaced by Bessel functions. The eigenfunctions become

$$ v_{k,j}(r,\theta) = B_j(\mu_k r) \,\sin{(j \theta)} $$

where $B_j$ is the $j$ -th order Bessel function and $\mu_k = \sqrt{\lambda_k}$. To find the eigenvalues we need to have the eigenfunctions vanish on the boundary of the disc. If the radius is one, we require

$$ B_j(\mu_k) = 0 $$

In other words, we need to compute zeros of Bessel functions. Here are the first nine eigenvalues and eigenfunctions of the circular disc. The violin string has become a tambourine.

```
eigenvals = [bjzeros(0,3) bjzeros(1,3) bjzeros(2,3)].^2
plot_modes('disc')
```

eigenvals = Columns 1 through 7 5.7832 30.4713 74.8870 14.6820 49.2185 103.4995 26.3746 Columns 8 through 9 70.8500 135.0207

Replace the full circular disc by a three-quarter circular sector. The angle at the origin is $270^\circ$ or $\frac{3}{2}\pi$ radians. We can make our eigenfunctions adapt to this angle. Take

$$ v_{k,j}(r,\theta) = B_{\alpha_j}(\mu_k r) \,\sin{(\alpha_j \theta)} $$

with

$$ \alpha_j = \frac{2}{3} j $$

and fractional order Bessel functions. The eigenfunctions satisfy their differential equation and also satisfy the boundary conditions on both sides of the angle.

$$ v_{k,j}(r,\theta) = 0 $$

at $\theta = 0$ and at $\theta = \frac{3}{2}\pi$.

By finding the zeros of the Bessel functions we can also have the eigenfunctions satisfy the boundary conditions on the outer circular portion of the boundary. Here are the first nine eigenvalues and eigenfunctions of the three-quarter circular sector.

```
eigenvals = [bjzeros(2/3,3) bjzeros(4/3,3) bjzeros(6/3,3)].^2
plot_modes('sector')
```

eigenvals = Columns 1 through 7 11.3947 42.6442 93.6362 18.2785 56.1131 113.6860 26.3746 Columns 8 through 9 70.8500 135.0207

These eigenfunctions have another important property. Most of them are singular; the derivatives of the fractional order Bessel functions are unbounded at the origin. You can see the black concentration of grid lines in the plots. If you tried to make a tambourine with this sector shape, it would rip at the sharp corner. This singular behavior is needed to model the solution to the wave equation on this region.

For all the regions we have discussed so far it is possible to express the eigenvalues as zeros of analytic functions. For the interval and the square, the eigenvalues are integers, which are the zeros of $\sin{\pi x}$. For the circular disc and sector, the eigenvalues are zeros of Bessel functions. Once we have the eigenvalues, it is easy to compute the eigenfunctions using sines and Bessel functions.

So, an L-shaped region formed from three unit squares is interesting for at least two reasons. It is the simplest geometry for which solutions to the wave equation cannot be expressed analytically; numerical computation is necessary. Furthermore, the 270 degree nonconvex corner causes a singularity in the solution. Mathematically, the gradient of the first eigenfunction is unbounded near the corner. Physically, a membrane stretched over such a region would rip at the corner. This singularity limits the accuracy of finite difference methods with uniform grids.

I used the L-shaped region as the primary example in my doctoral thesis fifty years ago. MathWorks has adopted a modified surface plot of the first eigenfunction as the company logo. I am going to devote a series of blog posts to the L.

Here are the first nine eigenvalues and eigenfunctions, computed by a function from Numerical Computing with MATLAB, which I will discuss in a future posting. Compare these eigenfunctions with the ones from the circular sector, which has the same reentrant corner and resulting singularity.

for k = 1:9 [~,eigenvals(k)] = membranetx(k); end eigenvals plot_modes('L')

eigenvals = Columns 1 through 7 9.6397 15.1973 19.7392 29.5215 31.9126 41.4745 44.9485 Columns 8 through 9 49.3480 49.3480

Simple model problems involving waves on an L-shaped region include an L-shaped membrane, or L-shaped tambourine, and a beach towel blowing in the wind, constrained by a picnic basket on one fourth of the towel.

A more practical example involves ridged microwave waveguides. One such device is a waveguide-to-coax adapter. The active region is the channel with the H-shaped cross section visible at the end of the adapter. The ridges increase the bandwidth of the guide at the expense of higher attenuation and lower power-handling capability. Symmetry of the H about the dotted lines shown in the contour plot of the electric field implies that only one quarter of the domain needs to be considered and that the resulting geometry is our L-shaped region. The boundary conditions are different than our membrane problem, but the differential equation and the solution techniques are the same.

The photo is courtesy of Advanced Technical Materials, Inc. See their website, <http://atmmicrowave.com/>, for lots of devices like this.

waveguide

Get
the MATLAB code

Published with MATLAB® R2014b

Greg's pick this week is Quad-Sim by David.

There have been a bunch of models for multiple-rotor flying machines floating around the File Exchange, but this entry is far and above the most comprehensive entry I have seen.

It has neat, sophisticated Simulink models of the dynamic system behavior; graphical user interfaces for data analysis, data formatting, and visualization; and finally, it includes the notion of model development and verification.

To me this entry represents the power MATLAB and Simulink can provide to the design of complex dynamic systems such as flying machines.

This video summarizes it very nicely.

Okay, maybe not everything. At some point you need to have substance in order to be useful. And documentation helps create the context in which to operate the tools. The Quad-Sim entry has all of these. The content of this entry is well organized in a way which helps detail what exactly is available within the content of the Quad-Sim tools.

In addition, the Simulink simulation models are well organized in terms of functional hierarchy. In other words, elements contained within the Simulink model are placed in well-named subsystems and libraries according to their functionality. There are clear partitions between inputs, controllers, plant dynamics, sensors, etc.

It sounds simple, but it makes all the difference. The organization of the content and the models alone made it very easy for me to understand, and start using the models.

I am definitely a proponent of Simulink as a platform for the development of dynamic systems and their control. Because I believe in Simulink as a design tool, it also means I spend a great deal of time and effort thinking about how to create and share models and supporting materials.

I have a self-imposed bias that makes it difficult for me to demonstrate or promote models that are messy, or otherwise difficult to understand. As such, I struggle to find File Exchange entries that involve Simulink, especially because Simulink is designed as a graphical tool, and therefore the form of the simulation model is at the forefront.

If David and his colleagues somehow believe that the world of development and design will always be well thought-out and documented, they will soon be in for a shock.

Generally deadlines of overlapping projects loom, and we must produce functionality to meet specific needs now, and worry about documentation later (if at all).

As all of us know, being a user of a tool we want clear, well though-out documentation that not only has details of how to use a particular tool or feature, but examples of its use. Finally understanding the underlying nature of a particular feature helps us better use features appropriately. But as many of us are some sort of developer or designer, we also are aware of the effort required to produce such rigorous documentation.

The Quad-Sim not only has well written documentation on specifically how to use the included tools (e.g. Modeling GUI):

It also has detailed explanations of various experimental processes so they can be reproduced, and the tools can be understood relative to a specific design context (e.g. Motor Modeling Hardware).

The fact that David and his colleagues bothered to attempt and verify their system models really excited me. I know it sounds like maybe I should reconsider my priorities in life, but I couldn’t help showing some of my colleagues in my hall the Quad-Sim entry. Nor could I hide my excitement in seeing users spend the time to make sure their models are correct.

A model is an invaluable tool for design and analysis of complex systems. But what good is the use of that model for design if it is “wrong”?

Certainly “wrong” is a relative term. A model is always wrong, otherwise it’s not a model. I am specifically referring to the notion that a model must capture the important attributes of a system sufficiently with enough accuracy to be useful for the purposes of the design and analysis.

Maybe I should say: what good is using a model that is “too wrong”?

It wasn’t clear to me that David and his colleagues used automated C-code generation from Simulink to implement the controller on the Arduino processor. Control systems like this are ideal candidates for our C-code generation technology from Simulink.

In my opinion, the Quad-Sim entry covers the more important aspects of the design process. It’s more important to get it right, without damaging people or equipment, than to implement quickly.

Seeing as that part of the process has been taken care of, it seems natural to me to automate the process of generating the control algorithm from the Simulink model to be implemented on the Arduino processor.

This enables two key design process attributes:

- Prevent the introduction of errors when translating designs from Simulink to C-code
- Permit design iterations to be tested in simulation, and then quickly implemented

Both of these help save demand on software development resources, especially as designs become more complex, and design cycle iterations become shorter. We do offer an Arduino Hardware Support Package that enables you to include connections of the model to hardware peripherals such as analog-to-digital converters (ADC), and go all the way from the Simulink model to the executable running on the hardware with a single button click.

This is fine if you can represent the entire system software in Simulink. But in most cases, engineers want to integrate the control algorithm with some existing software framework. For this, you can leverage Embedded Coder to generate the C-code and then integrate the resulting code into your software project.

I have lots of other comments and notions about this entry, but I’ll save that for other discussions. Let us know what you think.

**Comments**

If you would like to leave any comments regarding this post, please click here.

Get
the MATLAB code

Published with MATLAB® R2014b

We are very excited to announce that MATLAB r2014b introduces a new graphics system. Graphics developer, Mike Garrity, is launching a graphics specific blog. Here is the landing page for this release of MATLAB and for the new graphics system in general.]]>

I wanted to show you a glimpse of some of the new math functionality available in R2014b.

Recently on the MATLAB newsgroup, Christoph asked this question:

*I have a vector A shown below, which has 6 elements. the elements are already sorted in descending order. now i want to create vector C by deleting elements from A, starting with element a1, until the sum of the vector equals or is smaller the value B*

A= 26 23 20 19 15 14

B=70

*So, the output should be*

C= 20 19 15 14

*Any idea how to do this?*

*I would do something like this. Get the numbers in increasing order and perform a cumulative sum.*

Aflip = flipud(A); Asums = cumsum(Aflip);

*Find the first sum that is GREATER than b. The next element is where you want to start.* startIndex = find(Asums > b, 1, 'first');

*But things are reversed now. So we need to find the starting index into the original array.* Aind = length(A) - startIndex + 1;

*Now delete the leading values you don’t want.*

Afiltered = A((Aind+1):end);

In R2014b, the function `cumsum` has a new option, `direction`. (So do some other functions in the "|cum|" category.) You can specify to do the calculation in the default or `'forward'` direction, or in `'reverse'`. Here's what the new solution looks like.

A = [ 26 23 20 19 15 14]; B = 70; tmpA = cumsum(A,'reverse') ind = find(tmpA > B, 1,'last') C = A; C(1:ind) = []

tmpA = 117 91 68 48 29 14 ind = 2 C = 20 19 15 14

Have you had a chance to explore some of the math functionality in R2014b? Which features did you find helpful? Let me know here.

Get
the MATLAB code

Published with MATLAB® R2014b

I'm terrible at rote memorization.

I thought about that recently when I happened to see a plot (I think it was in a documentation example) that compared $sin(\theta)$, $cos(\theta)$, and $sin(\theta) cos(\theta)$.

theta = linspace(0,6*pi,500); plot(theta,sin(theta)) hold on plot(theta,cos(theta)) plot(theta,sin(theta) .* cos(theta)) hold off legend({'sin(\theta)','cos(\theta)','sin(\theta) cos(\theta)'})

Hmm, I thought. I guess $sin(\theta) cos(\theta)$ is just a sinusoid with twice the frequency and half the amplitude.

How is that related to memorization, you might ask? Well, I always had a very hard time remembering trig identities in high school and college. But, as an electrical engineering student, I had the following formula engraved on my brain cells:

$$e^{j\theta} = cos(\theta) + j sin(\theta)$$

(Note that I am following the proud electrical engineering tradition of using *j* instead of *i* as the imaginary unit. Because *i* obviously stands for current.)

This is called Euler's formula, and it's used all over the place in electrical engineering. I did not know at the time that Richard Feynman called it "the most remarkable formula in mathematics."

Here are a couple of closely related formulas:

$$cos(\theta) = \frac{1}{2} (e^{j\theta} + e^{-j\theta})$$

$$sin(\theta) = \frac{1}{2j} (e^{j\theta} - e^{-j\theta})$$

At some point during my undergraduate days, it dawned on me that I could derive many of the common trig identities directly from these formulas. Here's how it goes for $sin(\theta) cos(\theta)$:

$$ sin(\theta) cos(\theta) = \frac{1}{2j} (e^{j\theta} - e^{-j\theta}) \frac{1}{2} (e^{j\theta} + e^{-j\theta}) $$

$$ sin(\theta) cos(\theta) = \frac{1}{4j} (e^{j2\theta} - e^{-j2\theta} + 1 - 1) $$

$$ sin(\theta) cos(\theta) = \frac{1}{2} \frac{1}{2j} (e^{j2\theta} - e^{-j2\theta}) $$

$$ sin(\theta) cos(\theta) = \frac{1}{2} sin(2\theta) $$

Presto! Like magic.

I used to think that I must be kind of weird to prefer this manipulation of Euler's formula to simply memorizing the trig identities. But while writing this blog today, I noticed that the Wikipedia article on Euler's formula describes this concept exactly. "Complex exponentials can simplify trigonometry, because they are easier to manipulate than their sinusoidal components."

OK!

PS. In my plot creation code above, I exploited a subtle improvement in MATLAB R2014b graphics that customers have been requesting for many years. Can you spot it?

Get
the MATLAB code

Published with MATLAB® R2014b

**Fast Restart**

If you are dealing with large models, you are definitely aware that every time you click play to simulate a model, a good amount of time can be spent in initialization. In R2014b, if you know you will simulate your model multiple times without making structural changes, you can click the Fast Restart button:

At the end of the simulation, the model will remain in an Initialized state, and the toolbar will indicate this, telling you: "I am ready for fast restart".

At this point you can change the value of any tunable parameter and click play. The simulation will start instantly.

**Quick Insert**

Those who prefer typing over using the mouse will love the Quick Insert feature in R2014b. Type whatever you need in the canvas and the blocks magically appear! Once the block is there, a pop up offers to set the value of the most common parameter for this block:

**Model Templates**

You can now create models from templates:

This will launch the Simulink Template Gallery:

Of course, you can add your own templates to the gallery using the **File** menu of any model:

**Now it's your turn**

As you can see, R2014b contains many new features designed especially to speed up your workflow.

What is your favorite feature in R2014b? Let us know by leaving a comment here

]]>- Connect MATLAB Mobile to your computer with the MATLAB Connector. For more information on how to do this, refer to the getting started instructions here. This feature is only supported on MATLAB R2014a and later, so make sure you are on a compatible version.
- Install the MATLAB Support Package for Android Sensors. Choose
*Add-ons*from the MATLAB Toolstrip, and then choose*Get Hardware Support Packages*. This will open the support package installer. Choose*Android Sensors*from the list and follow the instructions. - To establish communication between the sensors on your device and MATLAB, create a mobiledev object, as follows:

m = mobiledev;

Once you have completed the 3 steps from the above section, go to MATLAB Mobile and turn on the accelerometer. You should see something akin to this:

You can also enable the sensor directly from MATLAB, by executing the following command:

m.AccelerationSensorEnabled = 1;

Did you notice the enabled

m.Logging = 1;You can verify this in MATLAB, note the Current Sensor Values in the result:

m

m = mobiledev with properties: Connected: 1 Logging: 1 InitialTimestamp: '02-Oct-2014 21:53:26.707' AccelerationSensorEnabled: 1 (20 Logged values) AngularVelocitySensorEnabled: 0 MagneticSensorEnabled: 0 OrientationSensorEnabled: 0 PositionSensorEnabled: 0 Current Sensor Values: Acceleration: [0.2631 5.9226 8.1850] (m/s^2)

Walk around your campus/home/floor with your device. Once you are satisfied, stop sending this data to MATLAB. You can either tap on the

m.Logging = 0;To retrieve the data, use the accellog variable:

[a, t] = accellog(m);

Once you have retrieved the logged acceleration data, you can plot it in MATLAB:

plot(t, a); legend('X', 'Y', 'Z'); xlabel('Relative time (s)'); ylabel('Acceleration (m/s^2)');

Calculate the magnitude to convert your X, Y and Z vectors to scalar values. Then, plot it.

x = a(:,1); y = a(:,2); z = a(:,3); % Calculate and plot magnitude. mag = sqrt(sum(x.^2 + y.^2 + z.^2, 2)); plot(t, mag); xlabel('Time (s)'); ylabel('Acceleration (m/s^2)');

To remove constant effects such as gravity, you can subtract the mean from this data.

% Accounting for gravity. magNoG = mag - mean(mag); % Plot magnitude. plot(t, magNoG); xlabel('Time (s)'); ylabel('Acceleration (m/s^2)');

The plotted data is now centered on zero, and shows peaks which correspond to a step taken while walking.

To determine the number of steps taken, you can use to FINDPEAKS function from Signal Processing Toolbox. In this example, we are treating only peaks with a minimum height above one standard deviation as a step. This threshold should be tuned experimentally to match a person’s level of movement while walking, hardness of floor surfaces etc.

% Use FINDPEAKS to determine the local maxima. minPeakHeight = std(magNoG); [pks, locs] = findpeaks(magNoG, 'MINPEAKHEIGHT', minPeakHeight);The number of steps taken is the number of peaks:

numSteps = numel(pks)

numSteps = 15Finally, you can also identify these locations on your plot of acceleration magnitude data:

hold on; % Place a red marker on the locations that correspond to peaks. plot(t(locs), pks, 'r', 'Marker', 'v', 'LineStyle', 'none'); title('Counting Steps'); xlabel('Time (s)'); ylabel('Acceleration Magnitude, Gravity Removed (m/s^2)'); hold off;

Once you are done, make sure you turn off the acceleration sensor and clear the mobiledev object.

```
m.AccelerationSensorEnabled = 0;
clear m;
```

- Acquiring Data from Sensors - MATLAB Mobile
- Android Sensor Support from MATLAB
- Documentation for mobiledev Object

Let us know by leaving a comment below.

Published with MATLAB® R2014b

]]>Mike is a computer graphics guru, well informed on the state of the art and knee-deep in the operational knowledge of How It All Works. He helped develop many of the features he’ll be describing, and he is a thoroughly entertaining writer. Beyond all this goodness, MATLAB graphics have taken a tremendous leap forward with the latest release, and Mike’s blog will give you an ideal perch to learn all about it. I recommend heading over there and seeing what he has to say.

]]>In the meantime, I'd like to invite you to head over to see Mike Garrity's new blog about MATLAB graphics. The MATLAB graphics system received a substantial update in R2014b, and Mike has a lot of great stuff to show you.

Welcome to the MATLAB Central blogs, Mike!

]]>Sean's pick this week is makeinstall and Make Install Technology by Norbert Marwan.

My pick this week is a tool that makes it easy to package your own toolboxes into a single installer file to distribute. For 12.5 years this file has served a great purpose.

Starting with MATLAB R2012b, we introduced Apps. Apps packaging allow you to take your user interface and package it into a single file. A dependency analysis is done and you can include any auxiliary files that might be necessary such as data. The generated app has version control and can be installed into another user's Apps tab. Apps, however, only have one entry point and are thus not good for toolbox packaging.

Now, beginning with R2014b, you can package your own toolbox from the MATLAB toolstrip and it will build a single install file like Norbert's makeinstall.

Let's look at the two approaches: here is a directory of a few of my utilities I'd like to package up.

currentdir = cd('myutils'); % navigate makeinstall % Build install cd(currentdir) % restore

Three files are generate in my directory:

I copied these to a new directory and ran `install`. It took me through some steps, including how to add it to path at startup, a few tips on what it does and how to use it.
At the end it recommends running the following to see the contents:

`helpwin myutils`

Which reveals the H1 lines for my functions!

I'm impressed. In addition to this, there's also the Make Install Technology white paper that explains how this all works in great detail.

Now let's use the R2014b functionality.

I'll select the folder and it will grab all of the files and their dependencies. Then I can put in any other meta information, documentation example etc. that I might have:

This builds a single install file which I can give to anyone:

This can be installed into the custom toolboxes.

It's automatically added to the path so those utilities can be used anywhere. I do hazard a word of caution that this may
make it easier to have shadowed files. For example: I now have three versions of those functions, which one is actually going
to be called? If I make an update to one, how do I ensure it's the right one? `which -all` will be your friend here.

which -all FEX

C\Documents\MATLAB\FEX.m C\Documents\MATLAB\myutils\FEX.m % Shadowed C\Documents\MATLAB\Toolboxes\myutils\FEX.m % Shadowed

So thank you Norbert for being ahead of the curve and providing us with the make install technology as well as for maintaining it over the last decade.

Give both approaches a try and let us know what you think here or leave a comment for Norbert.

Get
the MATLAB code

Published with MATLAB® R2014b

Today I’d like to introduce a guest blogger, David Garrison, who is a MATLAB Product Manager here at MathWorks. This is the first in a series of blogs over the next few weeks describing the new graphics system in R2014b and how some of the changes will affect you.

**Part 1: Features of the New Graphics System**- Part 2: Using Graphics Objects
- Part 3: Compatibility Considerations in the New Graphics System

Here is Part 1 of the series.

- Big Changes in R2014b
- The New MATLAB Graphics System
- The New Look of MATLAB Graphics
- Rotatable Tick Labels
- Automated Updating of
`datetime`Tick Labels - Animated Plots
- Multilingual Text and Symbols in Plots
- User Interfaces with Tab Panels
- Improved Histograms
- Have you tried the new graphics system in R2014b?
- Next up -- Part 2: Using Graphics Objects

There are a number of big changes in R2014b. Some of the new features include:

- New MATLAB graphics system
- Date and time data types with time zone and display options
- Git and Subversion source control integration and access to projects on GitHub from File Exchange
- MATLAB toolbox packaging as single, installable files for easy sharing and downloading of user-developed tools
- MATLAB MapReduce™ data analysis that scales Hadoop for big data
- Arduino and Android hardware support for interacting with motors and actuators, and for accessing sensor data

You can learn more about all these features in the MATLAB R2014b Release Notes. This post is all about the new graphics system.

R2014b includes a new MATLAB graphics system. The new graphics system includes many new features which we will describe in this blog post. In Part 2 of this series, we will describe how graphics handles have changed in R2014b and how to use graphics objects. Finally, there are some changes in the new system which may require changes in some existing graphics related code. We'll discuss compatibility considerations in Part 3 of this series.

When you create a plot in R2014b, you'll see that MATLAB graphics look different than previous versions.

The first thing that you will notice is that lines are plotted using a different set of colors. New line colors were selected to make it easier to distinguish lines from one another and to help people with certain types of color blindness. There is also a new default colormap in R2014b called *parula*. The colors in the parula colormap are ordered from dark to light and are perceptually uniform. Smooth changes in the data appear as smooth changes in color, while sharp changes in the data appear as sharp changes in color. Grid lines are now gray to make data stand out visually. Axis labels and titles are larger and more prominent. Lines and text are now anti-aliased (smoothed) to remove jagged edges.

If you're like me, you sometimes need to use long labels for the ticks in your plot. In previous versions of MATLAB, tick labels were always displayed horizontally as shown on the left in the example below. In the new graphics system, ticks can be rotated so you can create a plot like that shown on the right. Tick labels on both the x and y axes can be rotated.

For example, to rotate the X tick labels, use the `XTickLabelRotation` property of the `Axes` object:

ax = gca; ax.XTickLabelRotation = -40;

R2014b includes a new date and time data type called a `datetime` array. When you create a plot with a `datetime` array, the tick labels will automatically update as you pan or zoom in the plot. Tick labels change from days to hours to minutes to seconds as shown below.

MATLAB R2014b graphics introduces a new function called `animatedline`. The code below shows how to create an `animatedline` and add points to it.

x = 1:100; y = rand(1,100); myLine = animatedline; myLine.Color = blue; xlim([1 100]) ylim([0 1]) for i = 1:100 addpoints(myLine,x(i),y(i)) pause(0.1) end

The result is a plot that changes over time. The picture below shows how the plot looks at iterations 20, 50, and 80 as points are added.

R2014b graphics also supports the use of Unicode characters to show multilingual text in axis labels and titles and in user interface controls. Here is an example. See the `native2unicode` function for more information about Unicode strings.

Creating a user interface with tab panels has been a common request among MATLAB UI builders. In the past, people had done it using undocumented functions. In R2014b, those functions have been updated and are now fully documented and supported.

You can easily create a user interface like the one shown above using the new `uitabgroup` and `uitab` functions.

myFig = figure('Toolbar', 'none', ... % Create the figure 'Menubar', 'none', 'Name', 'Using Tab Panels'); tgroup = uitabgroup('Parent', myFig); % Create the tabgroup tab1 = uitab('Parent', tgroup, 'Title', 'Loan Data'); % Create the tabs tab2 = uitab('Parent', tgroup, 'Title', 'Amortization Table'); tab3 = uitab('Parent', tgroup, 'Title', 'Principal/Interest Plot');

The new `histogram` function plots histograms with data-dependent bin picking. It provides capabilities that the traditional `hist` function does not including options for bin control, normalization, and visualization.

The code below shows how to create two overlapping histograms and modify their properties after they are created.

data1 = randn(5000,1); data2 = randn(5000,1)+ 2;

h1 = histogram(data1); hold on h2 = histogram(data2); hold off legend show

h1.BinMethod = 'sturges'; h2.Normalization = 'countdensity';

Have you installed MATLAB R2014b? Have you tried the new graphics system? We'd love to hear your thoughts here.

Well, that's all for now. Be sure to check out Mike Garrity's new blog, Mike on MATLAB Graphics, for some cool ideas about what you can do with the new graphics system.

In my next post I'll talk about how graphics handles have changed in R2014b and how to use graphics objects.

Get
the MATLAB code

Published with MATLAB® R2014b

A fascinating question came to the Image Processing Toolbox development team recently from Brett Shoelson, a MathWorks application engineer, prolific File Exchange contributor, and former MATLAB Central blogger. Brett was trying to solve a problem for a customer. Part of the problem came down to this question:

Suppose you have a label matrix that represents a number of different regions. Are any of the regions disconnected?

Let me elaborate. A label matrix looks something like this:

```
L = [ ...
1 1 1 2 5 6
1 1 2 2 5 5
3 4 4 4 4 4
3 3 7 7 7 8
3 3 7 7 8 8
9 9 9 9 9 4]
```

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

For a positive integer k, the set of elements of `L` equal to k is the kth region. So, for example, the 5th region is:

L == 5

ans = 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Or, displayed as an image:

imshow(L == 5,'InitialMagnification','fit')

So my sample label matrix here represents 9 different regions. One of them, region 4 is disconnected, meaning it consists of more than one connected set (or connected component) or pixels.

L == 4

ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

imshow(L == 4,'InitialMagnification','fit') title('Region 4')

How can we programmatically determine whether any of the regions is disconnected? A brute force method is to compute, one at a time, the binary image associated with each region and then call `bwconncomp` to determine the number of connected components associated with it.

for k = 1:9 mask_k = (L == k); cc = bwconncomp(mask_k); if cc.NumObjects > 1 fprintf('Region %d is disconnected.\n',k) end end

Region 4 is disconnected.

If a brute force method gives you the desired answer and does it fast enough, then sometimes you just declare success and move on to solve the next problem of the day. Was this method fast enough? "No," said Brett. "My customer's images are about 25,000-by-25,000, and this method is way too slow."

"Hmm," I replied, trying to sound thoughtful while really just playing for time.

"OK, ask your customer this question: Is it known whether any of the regions in question contain any holes?"

Soon, the answer came back: "No, the labeled regions do not contain holes."

That was the answer I had hoped for. I started thinking about how to compute something called the Euler number for every labeled region, more or less simultaneously.

No, not that Euler number. I was thinking of the Euler number for a binary image, which is the number of objects in the image (or the number of connected components) minus the number of holes. (I don't know why this property of a binary image is called Euler number. If you know, please let us know in the comments.) The Image Processing Toolbox function `bweuler` computes this number for a binary image.

For a topological property of an image, the Euler number is unusual because it can be computed by one sweep of a purely local, 2-by-2 neighborhood operation. I first saw the method in the book *Robot Vision*, Berthold K. P. Horn, MIT Press, 1989. In this method, you count the number of occurrences of this 2-by-2 neighborhood pattern (called an *upstream convexity*):

0 0 0 1

and then subtract the number of occurrences of this 2-by-2 neighborhood pattern (called an *upstream concavity*):

0 1 1 1

(Note: this variation of the method works for 4-connected objects, not for 8-connected objects.)

But Brett has a label matrix representing many different regions, not a single binary image mask representing one region. How can we apply this technique?

To answer the question, let's look at the 2-by-2 neighborhood to the northwest of element (3,2):

L32 = L(2:3,1:2)

L32 = 1 1 3 4

I'll take the lower-right pixel as my reference point for the neighborhood. I can look at just the binary mask formed by comparing `L(3,2)` to its three neighbors to the north and west.

L(3,2) == L32

ans = 0 0 0 1

So, for region 4, we have an upstream convexity at this location. We could repeat this calculation in one loop over all the pixel locations, looking for (and counting) upstream convexities and upstream concavities associated with the label at the reference point for each 2-by-2 neighborhood.

I want to show you how to do this just with matrix operations and indexing. No loop is required. The technique is fairly general and can be applied to other local neighborhood operations.

Consider these two submatrices:

L1 = L(1:end-1,1:end-1)

L1 = 1 1 1 2 5 1 1 2 2 5 3 4 4 4 4 3 3 7 7 7 3 3 7 7 8

and

L2 = L(2:end,2:end)

L2 = 1 2 2 5 5 4 4 4 4 4 3 7 7 7 8 3 7 7 8 8 9 9 9 9 4

Note that the elements in `L1` are the northwest neighbors of the corresponding elements in `L2`. This basic idea gives us a general approach for doing certain kinds of neigbhorhood operations without explicit per-pixel loops.

Here's a short function I wrote that uses this technique to compute 4-connected Euler numbers for a label matrix:

```
type eulerNumbers4
```

function e = eulerNumbers4(L) % e = eulerNumbers4(L) % % For a label matrix L containing nonnegative integers, returns a vector e % such that e(k) is the 4-connected Euler number of the binary image (L == % k), for k = 1:max(L(:)). Lp = padarray(L,[1 1],0,'pre'); num_regions = max(Lp(:)); I_NW = Lp(1:end-1,1:end-1); I_N = Lp(1:end-1,2:end); I_W = Lp(2:end,1:end-1); is_upstream_convexity = (L > 0) & (L ~= I_N); is_upstream_convexity = is_upstream_convexity & (L ~= I_NW); is_upstream_convexity = is_upstream_convexity & (L ~= I_W); is_upstream_concavity = (L > 0) & (L ~= I_NW); is_upstream_concavity = is_upstream_concavity & (L == I_N); is_upstream_concavity = is_upstream_concavity & (L == I_W); upstream_convexity_labels = L(is_upstream_convexity); upstream_concavity_labels = L(is_upstream_concavity); total_upstream_convexities = accumarray(upstream_convexity_labels, ... 1,[num_regions 1]); total_upstream_concavities = accumarray(upstream_concavity_labels, ... 1,[num_regions 1]); e = total_upstream_convexities - total_upstream_concavities;

Under the assumption that the regions do not contain holes, the Euler number for each region is exactly the number of 4-connected components belong to that region.

eulerNumbers4(L)

ans = 1 1 1 2 1 1 1 1 1

So, we can finally answer the question: Are any of the labeled regions disconnected?

any(eulerNumbers4(L) > 1)

ans = 1

Yes.

True confessions, though. After I sent this version to Brett, and it didn't work for his problem, we realized that he needed the answer for 8-connected regions, not 4-connected.

It turns out that there is a local neighborhood-based Euler number algorithm for 8-connected regions, but it is significantly more complex to implement in this form. I did write it up, but I thought it would be too much uninteresting code to show here. If you are really interested, the method is described in Pratt, *Digital Image Processing*, 2nd edition, John Wiley & Sons, 1991, pp. 632-634. Pratt, in turn, credits the paper S. B. Gray, "Local Properties of Binary Images in Two Dimensions," *IEEE Transactions on Computers*, C-20, 5, May 1971, pp. 551-561.

The function `bweuler` is based on this algorithm, plus the binary image neighborhood lookup table approach that I have described here previously.

Have you used the Euler number measurement in your own work? Let us know in the comments.

Get
the MATLAB code

Published with MATLAB® R2014a

This is the third in a series of posts on the finite Fourier transform. The Fourier matrix produces an interesting graphic and has a surprising eigenvalue distribution.

The Fourier matrix of order $n$ is the $n$ -by- $n$ complex Vandermonde matrix $F$ whose elements $f_{k,j}$ are powers of the $n$ th root of unity

$$ \omega = e^{-2 \pi i/n} $$

$$ f_{k,j} = \omega^{k j} $$

The matrix can be generated with the MATLAB statements

k = (0:n-1)'; j = (0:n-1); F = exp(-2*pi*i*k*j/n);

Or, by taking the FFT of the identity matrix

F = fft(eye(n))

The statement

plot(F)

connects points in the complex plane whose coordinates are the real and imaginary parts of the elements of the columns of `F`, thereby generating a subgraph of the graph on `n` points. If `n` is prime, connecting the elements of all columns generates the complete graph on `n` points. If `n` is not prime, the sparsity of the graph of all columns is related to the computational complexity, and hence the speed, of the fast FFT algorithm. The graphs for `n` = 8, 9, 10, and 11 are generating and plotted by the following code.

for n = 8:11 subplot(2,2,n-7) F = fft(eye(n)); plot(F) axis square axis off title(['n = ' int2str(n)]) end

Because `n = 11` is prime, the corresponding graph shows all possible connections. But the other three values of $n$ are not prime. Some of the links in their graphs are missing, indicating that the FFT of a vector with that many points can be computed more quickly.

The remainder of this post examines `F12` in more detail. Here is its graph.

clf n = 12; F = fft(eye(n)); plot(F) axis square axis off title(['n = ' int2str(n)])

The program `fftmatrix`, available here, or included with the NCM App, allows you to investigate these graphs. `fftmatrix(n)` plots all the columns of the Fourier matrix of order `n`. `fftmatrix(n,j)` plots just one column.

Let's plot the individual columns of `F12`. The first column of `F12` is all ones, so its plot is just a single point.

for j = 1:12 p = mod(j-1,4)+1; subplot(2,2,p); fftmatrix_mod(12,j-1) title(['j = ' int2str(j)]) if p == 4, snapnow, end end

To see typical behavior, look at the third subplot, the red graph labeled `j = 3`, generated by the third column

F(:,3)

ans = 1.0000 + 0.0000i 0.5000 - 0.8660i -0.5000 - 0.8660i -1.0000 + 0.0000i -0.5000 + 0.8660i 0.5000 + 0.8660i 1.0000 + 0.0000i 0.5000 - 0.8660i -0.5000 - 0.8660i -1.0000 + 0.0000i -0.5000 + 0.8660i 0.5000 + 0.8660i

These are the powers of $\omega^2$. Because 2 divides 12 evenly these powers hit all the even-numbered points around the circle twice and miss all the odd-numbered points. Now look at the cyan graph labeled `j = 11`. These are the powers of $\omega^{10}$, which are the complex conjugates of the powers of $\omega^2$. So the two graphs lie on top of each other.

Six of the twelve columns in the graph of `F12` connect only a subset of the nodes and ten of the columns lie on top of their complex conjugate columns. As a result, when all of the columns are combined to form the full graph, it is sparse. This sparsity, in turn, makes it possible to construct a fast finite Fourier transform algorithm for `n = 12`.

I've always been curious about the eigenvalues and vectors of the Fourier matrices. In 1979, three friends of mine at the University of New Mexico, Gus Efroymson, Art Steger, and Stan Steinberg, posed the question in the SIAM Review problem section. They didn't know that Jim McClellan and Tom Parks had actually solved their problem seven years earlier, when Jim was a grad student working under Tom at Rice. I didn't learn about the McClellan-Parks paper until fairly recently.

The Wikipedia page on the discrete Fourier transform says the facts about the eigenvalues are well known, but that the eigenvectors are the subject of ongoing research.

Let's rescale $F$ so that its columns have unit length.

F = F/sqrt(n);

This makes $F' \ F = I$

round(F'*F)

ans = 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1

Now $F$ is *unitary*, the complex generalization of *orthogonal*. This implies that all of its eigenvalues lie on the unit circle in the complex plane.

Furthermore, it turns out that $F^4 = I$.

round(F^4)

ans = 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1

This implies that any eigenvalue $\lambda$ must satisfy $\lambda^4$ = 1. Consequently the only possible eigenvalues are 1, -1, i, and -i. You might guess that guess that if $n$ is divisible by 4, the eigenvalues would be equally distributed among these four values. But, surprisingly, to me at least, that doesn't happen.

lambda = eig(F)

lambda = 1.0000 + 0.0000i -1.0000 + 0.0000i -0.0000 + 1.0000i -1.0000 + 0.0000i -1.0000 - 0.0000i -0.0000 + 1.0000i 0.0000 - 1.0000i 1.0000 + 0.0000i -0.0000 - 1.0000i 1.0000 - 0.0000i 1.0000 - 0.0000i 0.0000 - 1.0000i

It's hard to pick through this unsorted output, but there are four +1's, three -1's, three -i's, and only two +i's.

Here is a tricky piece of code that uses `angle` and the counting feature of sparse indexing to count the number of each of the four possible eigenvalues.

```
type eigfftmat
```

function c = eigfftmat(n) % EIGFFTMAT Count eigenvalues of the Fourier matrix. % c = eigfftmat(n) is a 4-vector with counts for +1, -1, -i, +i. % Compute the eigenvalues. e = eig(fft(eye(n))); % Sparse repeated indexing acts as a counter. c = full(sparse(mod(round(angle(e')/(pi/2)),4)+1,1,1))'; c([3 2]) = c([2 3]); end

When we run this code for a sequence of values of `n`, we see the pattern predicted by the McClellan-Parks analysis. The number of each of the four possible eigenvalues depends upon `floor(n/4)` and `mod(n,4)`.

disp(' n +1 -1 -i +i') for n = 4:20 disp([n eigfftmat(n)]) end

n +1 -1 -i +i 4 2 1 1 5 2 1 1 1 6 2 2 1 1 7 2 2 2 1 8 3 2 2 1 9 3 2 2 2 10 3 3 2 2 11 3 3 3 2 12 4 3 3 2 13 4 3 3 3 14 4 4 3 3 15 4 4 4 3 16 5 4 4 3 17 5 4 4 4 18 5 5 4 4 19 5 5 5 4 20 6 5 5 4

The proofs in the McClellan and Parks paper involve the eigenvectors and are quite complicated. It turns out that this MATLAB expression

floor((n+[4 2 1 -1])/4)

generates a 4-vector of the multiplicities of the +1, -1, -i, and +i eigenvalues for any given value of `n`.

J. H. McClellan and T. W. Parks, "Eigenvalues and eigenvectors of the discrete Fourier transformation", IEEE Trans. Audio Electroacoust 20, 66-74, <http://dx.doi.org/10.1109/TAU.1972.1162342>, 1972.

G. Efroymson, A. Steger, and S.Steinberg, "A Matrix Eigenvalue Problem", SIAM Review, 21, 139-139, <http://dx.doi.org/10.1137/1021013>, 1979.

Wikipedia, "Discrete Fourier Transform", <http://en.wikipedia.org/wiki/Discrete_Fourier_transform>, retrieved 09/03/2014.

Get
the MATLAB code

Published with MATLAB® R2014a

Let’s see how this works using a simple example!

**The Problem**

Let's say we have a setup where a DC motor is used to move a load between two positions. The position is measured using an encoder, and a PID is used to control the angle of the DC motor. The model looks like:

And when it simulates, the load tracks the desired position in that manner:

How could we detect if something breaks in the motor, or if the dynamics of the load change?

**Online Identification**

When starting the system we will spend a few seconds to identify its dynamics online (while the system is running). Once we have a model estimated, we will compare the encoder reading with the output of the model. If the error between the real system and the model becomes too large, it's a sign something undesired is happening... and we want to detect that!

I wrapped the control loop shown above in a subsystem and connected to it the Recursive Polynomial Model Estimator block from the System Identification Toolbox library.

Using the enable input port, we enable estimation for the first few seconds (only a few periods of the motor motion are necessary to get a good model). Once the estimation is completed, we monitor the estimation error, using very simple logic blocks to make a signal pass from zero to one as soon as the error level becomes larger than an acceptable level.

**The Results**

To test the logic in simulation, I faked the failure by doubling the friction in my DC motor model. As you can see, when the system dynamics change around 10 seconds, the error flag instantly signals the problem:

**What's next?**

The online estimation block also supports code generation. To see that in action, I recommend watching the video Online Fault Detection for a DC Motor by my colleague Karthik Srinivasan.

In this video, Karthik uses a model similar to the one I describe above, except that instead of a simulated DC Motor, he used driver blocks from the Arduino Simulink Suppport Package, to detect a failure using a real DC motor.

**Now it's your turn**

As you can imagine, fault detection is only one application of online parameters estimation. For example the identified model could be used in any sort of adaptive control.

Let us know what you think of the new online estimation features of the System Identification Toolbox by leaving a comment here.

]]>Jiro's pick this week is Colormouse by Patrick.

In addition to examining the actual values or creating line plots, you can get a lot of insight by visualizing data with color.
MATLAB has a number of functions that let you represent your data with color, including `image`, `surf`, `mesh`, and `contourf`. Many of these functions make use of a `colormap` that represents a range of colors for the values.

contourf(peaks) colorbar

`colorbar` above shows the range of colors represented in the figure. You can programmatically change this color range by using the
`caxis` function.

caxis([-4 4])

Patrick's `colormouse` allows you to change the colormap range interactively using the mouse. Calling `colormouse` adds a toolbar button to the current figure.

Once you're in colormouse mode, you can click-n-drag to move (by moving the mouse up and down) or widen/narrow (left and right) the color range.

You can also select a specific colormap from a context menu.

**Comments**

Let us know what you think here or leave a comment for Patrick.

*Note that the actual sourcecode resides on GitHub, so please check the GitHub page for the license information.*

Get
the MATLAB code

Published with MATLAB® R2014a

I'd like to introduce this week's guest blogger Alan Weiss. Alan writes documentation for mathematical toolboxes here at MathWorks.

An economist I know (well, he's my son) asked me a question the other day. I was able to answer him quickly by using Symbolic Math Toolbox™ calculations. You might find the steps to be an interesting case of a somewhat nonstandard use of symbolic mathematics

My son was working through a paper that described a computational algorithm. He got stuck on one step of the algorithm:

If $A$ is a positive definite (symmetric) matrix, find a **lower triangular** matrix $G$ such that

$$GAG^{-1} = D, \mbox{ where }D\mbox{ is diagonal.}$$

He called me asking if there was some linear algebra that he just never learned. I knew how to find a matrix $G$ that would diagonalize $A$ using the `eig` function. But a lower triangular version? I didn't know how to do that. I told him I would call him back after some investigation.

After spending a couple of minutes fruitlessly looking through a linear algebra book, I decided to see what Symbolic Math Toolbox had to say about the equations that any such matrix $G$ would have to satisfy. Here were my steps.

syms a b c g1 g2 g3 real % Create symbolic variables G = [g1,0;g2,g3] % Lower triangular matrix H = inv(G) % Its inverse A = [a,b;b,c] % Symmetric matrix A D = G*A*H % D is supposed to be diagonal

G = [ g1, 0] [ g2, g3] H = [ 1/g1, 0] [ -g2/(g1*g3), 1/g3] A = [ a, b] [ b, c] D = [ a - (b*g2)/g3, (b*g1)/g3] [ (a*g2 + b*g3)/g1 - (g2*(b*g2 + c*g3))/(g1*g3), (b*g2 + c*g3)/g3]

Look at the upper right entry in D. If D is diagonal, then b*g1 = 0. This means b = 0 or g1 = 0. But if b = 0 then A is already diagonal. And if g1 = 0 then the first row of G is zero, so G is not invertible.

In other words, I had just proved that there is no such thing as a lower triangular matrix that diagonalizes a symmetric matrix in this way.

I called my son back and explained the calculation. We decided that the paper he was reading had a misstatement, and he went ahead and used `eig` to diagonalize the matrix in his application. He will get in touch with the author of the paper to explain the misstatement in the algorithm.

I was quite happy with this result. Symbolic math calculations had saved us both a lot of time searching for an algorithm that doesn't exist.

Do you use symbolic math calculations in cute or innovative ways? Tell us about it here.

Get
the MATLAB code

Published with MATLAB® R2014a

Idin's pick for this week is the SimRF Models of Analog Devices RF Agile Transceivers by MathWorks RF Team.

This is probably one of the best documented submissions on the File Exchange. The Simulink models in this submission can be used to accurately simulate a commercial radio frequency (RF) chip (Analog Devices AD9361). The work was the result of collaboration between MathWorks and Analog Devices engineers, and the final result is a software model that matches the behavior of the physical device on the lab bench.

You can get a lot of details on the hardware and these software models here.

I will not go into any further technical details on this page, but I want to point out why I believe this was a useful and interesting project.

First, this is just plain cool! Seeing any software model predict (or match) hardware behavior at the level of detail that we find in these models is always interesting (having run simulations for 15+ years, I’m still amused when what I see in the lab matches what I predicted in simulation; it feels magical every time!)

Second, it helps understand and program the chip. An RF chip like the AD9361 (or most any chip for that matter) is a complex piece of hardware, and you can't see inside of it. This is an "agile" RF transceiver, which means it can be programmed to perform many different tasks for many different radio applications. Understanding the functioning of such a chip, and programming it correctly are rather difficult tasks in hardware, especially when the results of one's work can't be seen immediately (at least not without the use of expensive lab equipment such as spectrum analyzers). Having a software model that accurately mimics the hardware allows me to understand the behavior of the hardware such that I can program it better.

Third, debugging. If you have programmed anything but the simplest piece of hardware, you know that things rarely do what you want on the first try. Understanding apparent anomalies and hardware "misbehavior" is made infinitely easier when we can go to a simulation model and trace different signals to our heart"s content! With this Simulink model I can look at intermediate signals that I could never access on the physical chip. These intermediate signals often hold the key to understanding why the system’s outputs look the way they do.

Fourth (this should really be 3.5), studying the effects of an RF chip on digital baseband signals. These effects are generally something that we observe once we put everything together in the lab. Being able to simulate the digital baseband with the RF front-end will be quite useful in designing a robust communication system before we put any hardware together. Remember that as cool as physical prototypes are, they're typically quite expensive and time-consuming to build.

Finally (or fifth!), it builds trust in Simulink and SimRF. This one explains why we invested the time to do this. The fact that we could construct this model in a relatively short amount of time, and that it matches lab measurements, makes us confident that the simulation software tools are both capable and accurate.

I encourage you to take a look at the documentation and videos on this page to get more information on this package.

As always, your thoughts and comments here are greatly appreciated.

Get
the MATLAB code

Published with MATLAB® R2014b

Since I am really enjoying the book, I decided that for this week's post I would use Simulink to validate some of the numbers provided by the author to one of the question in the book: If an asteroid was very small but supermassive, could you really live on it like the Little Prince?

**What would the gravity feel like on this little planet?**

For a tiny asteroid of 1.75m radius like the one in Le Petit Prince, by Antoine de Saint-Exupéry to have a gravitational field comparable to the one on Earth, it would need to weight about 500 million tons.

But how would this gravitational field change as you get farther form the surface? Since the radius is so small, the effect of the gravity would quickly go down, as illustrated by Randall:

*Source: http://what-if.xkcd.com/68/*

To validate that, I used a new R2014a feature in SimMechanics: The Gravitational Field block. To test the effect of the gravity, I simply linked together the Gravitational Field of 500 million tons and a Solid of 65kg with a Prismatic Joint, and slowly increased the distance between them. As you can see, the SimMechanics results match the book pretty well!

**What About the Escape Velocity?**

The next thing mentioned in the book is that if you were able to run fast enough (~5m/s) and jump, you could escape the asteroid gravity field:

*Source: http://what-if.xkcd.com/68/*

In SimMechanics, I replaced the Prismatic Joint by a Rectangular Joint:

I specified appropriate initial conditions for the joint:

I clicked play and I could see in the Mechanics Explorer that for the chosen speed, I would end up in a pretty steady almost circular orbit:

**How stable would that orbit be?**

The book mentions that the orbit would be *weird*...

*Source: http://what-if.xkcd.com/68/*

To test that, I replaced the Rectangular Joint by a Planar Joint, to allow rotation. Then I dug out the skater I had built in a previous post about the Olympics and without any modification made him move his arms while in orbit... definitely, the trajectory is interesting:

**Now it's your turn**

I guess this blog post highlights the fact that you should always trust what you read on XKCD!

Thanks to the Gravitational Field block of SimMechanics, I can now begin to plan space missions without digging out my notes from the orbital mechanics class I took during grad school!

Let us know what you think by leaving a comment here.

]]>This is the second in a series of three posts about the Finite Fourier Transform. This post is about the fast FFT algorithm itself. A recursive divide and conquer algorithm is implemented in an elegant MATLAB function named `ffttx`.

The acronym FFT is ambiguous. The first F stands for both "fast" and "finite". A more accurate abbreviation would be FFFT, but nobody wants to use that. In MATLAB, the expression `fft(x)` computes the finite Fourier transform of any vector `x`. The computation is fastest if the integer `n = length(x)` is the product of powers of small primes.

Direct application of the definition

$$ Y_k = \sum_{j=0}^{n-1} \omega^{jk} y_j, \ k = 0, \ldots, n-1 $$

requires $n$ multiplications and $n$ additions for each of the $n$ components of $Y$ for a total of $2 n^2$ floating-point operations. This does not include the generation of the powers of $\omega$. A computer capable of doing one multiplication and addition every microsecond would require a million seconds, or about 11.5 days, to do a million-point FFT.

Several people discovered fast FFT algorithms independently and many people have since contributed to their development, but it was a 1965 paper by John Tukey of Princeton University and John Cooley of IBM Research that is generally credited as the starting point for the modern usage of the FFT.

Modern fast FFT algorithms have computational complexity $O(n \mbox{ log}_2 n)$ instead of $O(n^2)$. If $n$ is a power of 2, a one-dimensional FFT of length $n$ requires fewer than $3 n \mbox{ log}_2 n$ floating-point operations. For $n = 2^{20}$, that's a factor of almost 35,000 faster than $2 n^2$. Even if $n = 2^{10}$, the factor is about 70.

With MATLAB 8.3 and a 2.7 GMHz Intel Core i7 laptop, the time required for `fft(x)` if `length(x)` is `2^23 = 8388608` is about 0.5 second. The built-in `fft` function is based on FFTW, "The Fastest Fourier Transform in the West," developed at M.I.T. by Matteo Frigo and Steven G. Johnson in 1998.

The key to the fast FFT algorithms is the fact that the square of the $2n$ th root of unity is the $n$ th root of unity. Using complex notation,

$$ \omega = \omega_n = e^{-2 \pi i/n} $$

we have

$$ \omega_{2n}^2 = \omega_n $$

The derivation of the fast algorithm starts with the definition of the finite Fourier transform:

$$ Y_k = \sum_{j=0}^{n-1} \omega^{jk} y_j, \ k = 0, \ldots, n-1 $$

Divide the sum into terms with even subscripts and terms with odd subscripts.

$$ Y_k = \sum_{even \ j} \omega^{jk} y_j \ + \ \sum_{ odd \ j} \omega^{jk} y_j $$

Assume that $n$ is even, that $k \leq n/2-1$ and repurpose the summation subscripts.

$$ Y_k = \sum_{j=0}^{n/2-1} \omega^{2jk} y_{2j} \ + \ \omega^k \sum_{j=0}^{n/2-1} \omega^{2jk} y_{2j+1} $$

The two sums are components of the FFTs of length $n/2$ of the portions of $y$ with even and odd subscripts. In order to get the entire FFT of length $n$, we have to do two FFTs of length $n/2$, multiply one of these by powers of $\omega$, and concatenate the results.

The relationship between an FFT of length $n$ and two FFTs of length $n/2$ can be expressed compactly in MATLAB.

omega = exp(-2*pi*i/n); k = (0:n/2-1)'; w = omega .^ k; u = fft(y(1:2:n-1)); v = w.*fft(y(2:2:n));

then

fft(y) = [u+v; u-v];

Now, if $n$ is not only even but actually a power of 2, the process can be repeated. The FFT of length $n$ is expressed in terms of two FFTs of length $n/2$, then four FFTs of length $n/4$, then eight FFTs of length $n/8$, and so on until we reach $n$ FFTs of length one. An FFT of length one is just the number itself. If $n = 2^p$, the number of steps in the recursion is $p$. There is $O(n)$ work at each step, so the total amount of work is

$$ O(np) = O(n \mbox{ log}_2 n) $$

If $n$ is not a power of two, it is still possible to express the FFT of length $n$ in terms of several shorter FFTs. An FFT of length 100 is two FFTs of length 50, or four FFTs of length 25. An FFT of length 25 can be expressed in terms of five FFTs of length 5. If $n$ is not a prime number, an FFT of length $n$ can be expressed in terms of FFTs whose lengths divide $n$. Even if $n$ is prime, it is possible to embed the FFT in another whose length can be factored. We do not go into the details of these algorithms here.

The textbook function `ffttx` combines two basic ideas. If $n$ is a power of 2, it uses the $O(n \mbox{ log}_2 n)$ fast algorithm. If $n$ has an odd factor, it uses the fast recursion until it reaches an odd length, then sets up the discrete Fourier matrix and uses matrix-vector multiplication.

`ffttx` is one of my favorite MATLAB functions. Here is the entire code.

function y = ffttx(x) %FFTTX Textbook Fast Finite Fourier Transform. % FFTTX(X) computes the same finite Fourier transform % as FFT(X). The code uses a recursive divide and % conquer algorithm for even order and matrix-vector % multiplication for odd order. If length(X) is m*p % where m is odd and p is a power of 2, the computational % complexity of this approach is O(m^2)*O(p*log2(p)).

x = x(:); n = length(x); omega = exp(-2*pi*i/n);

if rem(n,2) == 0 % Recursive divide and conquer k = (0:n/2-1)'; w = omega .^ k; u = ffttx(x(1:2:n-1)); v = w.*ffttx(x(2:2:n)); y = [u+v; u-v]; else % The Fourier matrix j = 0:n-1; k = j'; F = omega .^ (k*j); y = F*x; end

Here is a plot of the Fourier matrix of order 16. Plots like this are the subject of the next post in this series. For now, just notice that this is *not* a complete graph with 16 nodes.

plot(fft(eye(16))) axis equal axis off

Steve Eddins has written extensively about the FFT in his blog. Select "Categories more" and then click on "Fourier transforms". <http://blogs.mathworks.com/steve>

James W. Cooley and John W. Tukey, "An algorithm for the machine calculation of complex Fourier series", Math. Comput. 19: 297-301, 1965. <http://dx.doi.org/10.2307%2F2003354>

Matteo Frigo and Steven G. Johnson, "FFTW, Fastest Fourier Transform in the West", <http://www.fftw.org>.

Get
the MATLAB code

Published with MATLAB® R2014a

*Marta is a senior application engineer at MathWorks UK, helping MATLAB users to design, develop and deploy production-quality code. Toshi is a Senior Marketing Manager at MathWorks.*

*by Marta Wilczkowiak and Toshi Takeuchi*

In the previous post, you learned how we built an RFID-based people tracking system from scratch in three months, a “technology stack” that spanned multiple hardware and software disciplines. We ended the post with how we managed to survive the moment of the truth: the first integration test, live at a conference with hundreds of delegates.

Now we had raw data to analyze to improve our system.

**The Twist – Making Sense of RFID Data**

Once we examined the real word data collected during the event more closely, we faced a new dilemma – we had originally assumed we could use signal strength to estimate distances. However, in reality there was so much fluctuation in signal strength that you could not make a confident estimate this way. This, incidentally, is why RFID-based localization is a hot research topic today.

We were determined to work on an engineering solution, but it would take more time and effort. In the meantime, was there a way to get meaningful insights from the messy data at hand?

**Web Analytics Meets Engineering **

As we discussed the conference data internally at MathWorks, this issue came to the attention of Toshi Takeuchi, a US-based analyst. He immediately saw a similarity with the familiar problems he faced in analyzing massive volume of page views and clicks – always messy. What would happen if you applied web analytics techniques to this engineering problem? We embraced this suggestion and immediately shared our code and data with Toshi. Our collaboration was possible because we all spoke the same language – MATLAB.

**Use of Proxy Metric**

Non-physical measurements are often used in web analytics. For example, the recommendations you see on online shopping or movie streaming websites are computed using distance measured in terms of similarity or affinity. Was there anything in the dataset, besides the signal strength, that we could turn into a proxy metric? The new approach Toshi came up with was to recast connection frequency as this virtual distance. If two tracking devices are connecting to one another more frequently than others in a given time interval, then we could think of them being “close” in the virtual space defined by connection frequency. This metric seemed far more reliable than signal strength. Our original goal was to understand the social interactions of the delegates, and physical distance metric was just a means to that end. If a virtual metric did the same job, then it seemed a good substitute.

**The Result**

With the support from Toshi, we were able to reconstruct the estimated movements of the delegates using the connection frequencies with stationary anchors such as demo stations (orange circles) and catering (green circles). Note that no personally identifiable information was collected or used for this analysis.

In the plot above, you see such an estimated movement for Delegate #120. This person appeared to have approached a nearby catering station first, then spent quite some time around BLOODHOUND Super Sonic Car in the middle, but in the afternoon stayed “closer” to software related demos than hardware-based demos as measured in this “virtual” space.

We can use this information to map the frequencies of connections between delegates and demo stations, and visualize it as a heat map. The vertical bands with lighter red stripes indicates which demos had more connections than others, and horizontal stripes show to which demo stations each delegate connected frequently. You can see that people who connected with “RT” (a real-time testing exhibit) and “HW” (a hardware testing exhibit) had very little overlap with the people who connected with “Robo” (a robotics exhibit).

**The Insight – the Three Factors of Delegate Behavior **

Ultimately, the data is only as good as the insights you get from it. When we computed the estimated movements of the delegates with a technique called Principal Component Analysis, we could summarize the complex data by a mere three factors.

- The “Web”: people who showed strong “connections” with the “Web” demo
- The “Robo”: people who showed strong “connections” with “Robo” demo
- The “Conventional”: people who showed strong “connections” with “HW” demo

Furthermore, “Web” people and “Robo” people tend to share similar connection patterns as compared to “Conventional” people – what is common to both? They are both what we might call “convergence” applications where software and hardware are being integrated in an unconventional way.

While this might not be a novel insight, it was still very useful for us to know that the data backed up what we had been intuitively sensing all along. It was also interesting to see which delegates showed which particular inclinations on an individual basis (each dot in the plot identify a specific delegate).

**Back to Reality – the Next Step**

It is tempting to speculate if the virtual proximity on this plot correlates to actual shared interest among those delegates. If so, it is also interesting that two of the main factors turned out to be the web and robotics, which are examples of convergence-oriented applications that cross the traditional virtual/physical divide.

If there was an exit survey that asked the participants which demos they found interesting, we could use predictive analytics techniques to validate how well this virtual metric predicts the actual interest. Alternatively, we could use additional stationary anchors to determine which sessions delegates attended. If we placed more anchors on the floors in a more evenly spaced grid, including the areas populated by third party exhibitors, it would give more data points and more accurate measurements.

If we validated this correlation, it could open up opportunities for new ways for our delegates to find like-minded people.

**The Journey Continues**

Our journey is not over yet – we still have challenges with the system. But our project proved how such an inter-disciplinary system could be built in three months using MATLAB, Simulink, and lots of team work.

]]>**The Problem**

If you ask me what is the difference between Simulink Coder and Embedded Coder, I would tell you that Simulink Coder allows you to generate code from a Simulink model, and Embedded Coder allows you to configure how the code looks like.

To illustrate that, we will start with a simple hand-written program, and see how to configure a model so that the code generated from it integrates without modification.

In this simple `main`, at every second, the program reads data from a text file and stores it in a global variable `u`. We want our auto-generated code to access this variable, and use it to compute value of another global variable `y` to be used later in the code.

**Configuring the Model**

For this example, let's use the following simple model:.

By default, if we generate code for this model we get something that looks like:

It is obvious that modifications are necessary to map the values of `u` and `y` in the hand-written code to the input and output of the generated code.

To configure the look of the code, the first step is to name the input signal `u` and the output `y`. Then we need to tell Embedded Coder that `u` and `y` are external variables, already defined outside of the generated code, in the hand-written code.

For that, right-click on the signal line and select Properties:

Go to the Code Generation tab, and set the storage class to `ImportedExtern` (Note that the storage class and a lot more could also have been specified using a data object).

**Generated Code**

Generate code for the model, and you should get something which looks like this:

Where `u` and `y` are declared as `extern`:

With this modification, the generated code integrates in the hand-written application, I can build the `main` program without errors.

**Now it's your turn**

There are many of ways to customize the code generated by Embedded Coder. This example is probably the simplest one possible, but I hope this gives you a good idea of how to get started. Let us know how you got started with the coders by leaving a comment here.

]]>Many of us carry around smartphones that can track our GPS positions and that's an interesting source of data. How can we analyze GPS data in MATLAB?

Today's guest blogger, Toshi Takeuchi, would like to share an analysis of a public GPS dataset from a popular ride sharing service Uber.

Uber is a ride sharing service that connects passengers with private drivers through a mobile app and takes care of payment. They are in fact so popular that you hear about them in the news due to their conflicts with local traffic regulations and taxi business interests.

Uber’s ride sharing GPS data was available publicly on infochimps.com, so I used it for this analysis (Unfortunately it is not available anymore). What can we learn from this dataset?

Let's start by downloading the dataset from the link above (a zipped TSV file), which contains the GPS logs taken from the mobile apps in Uber cars that were actively transporting passengers in San Francisco. The data have been anonymized by removing names, trip start and end points. The dates were also substituted. Weekdays and time of day are still intact.

For the purpose of this analysis, let's focus on the data captured in the city proper and visualize it with Mapping Toolbox.

Run the script to load data. Check `loadData.m` to see the details.

loadData

Overlay the GPS points on the map.

states = geoshape(shaperead('usastatehi', 'UseGeoCoords', true)); latlim = [min(T.Lat) max(T.Lat)]; lonlim = [min(T.Lon) max(T.Lon)]; ocean = [0.7 0.8 1]; land = [0.9 0.9 0.8]; figure ax = usamap(latlim, lonlim); setm(ax, 'FFaceColor', ocean) geoshow(states,'FaceColor',land) geoshow(T.Lat,T.Lon,'DisplayType','Point','Marker','.',... 'MarkerSize',4,'MarkerEdgeColor',[0 0 1]) title('Uber GPS Log Data') xlabel('San Francisco') textm(37.802069,-122.446618,'Marina') textm(37.808376,-122.426105,'Fishermans Wharf') textm(37.797322,-122.482409,'Presidio') textm(37.774546,-122.412329,'SOMA') textm(37.770731,-122.440481,'Haight') textm(37.818276,-122.498546,'Golden Gate Bridge') textm(37.819632,-122.376065,'Bay Bridge')

Let's start with a basic question - how does the use of Uber service change over time. We can use `grpstats` to summarize data grouped by specific categorical values, such as `DayName` and `TimeOfDay`, which were added in the data loading process.

Get grouped summaries.

byDay = grpstats(T(:,{'Lat','Lon','DayName'}),'DayName'); byDayTime = grpstats(T(:,{'Lat','Lon','TimeOfDay','DayName'}),... {'DayName','TimeOfDay'});

Reshape the count of entries into a 24x7 matrix.

byDayTimeCount = reshape(byDayTime.GroupCount,24,7)';

Plot the data by day of week and by hours per day of week.

figure subplot(2,1,1) bar(byDay.GroupCount); set(gca,'XTick',1:7,'XTickLabel',cellstr(byDay.DayName)); subplot(2,1,2) plot(byDayTimeCount'); set(gca,'XTick',1:24); xlabel('Hours by Day of Week'); legend('Mon','Tue','Wed','Thu','Fri','Sat','Sun',... 'Orientation','Horizontal','Location','SouthOutside')

It looks like the usage goes up during the weekend (Friday through Sunday) and usage peaks in early hours of the day. San Francisco has a very active night life!

Is there a way to figure out where people go during the weekend? Even though the dataset doesn't contain the actual starting and ending points of individual trips, we may still get a sense of how the traffic flows by looking at the first and last points of each record.

We can extract the starting and ending location data for weekend rides. Click `getStartEndPoints.m` to see how it is done. If you would like to run this script, please download `districts.xlsx` as well.

% Here we load the preprocessed data |startEnd.mat| to save time and plot % their starting points. % getStartEndPoints % commented out to save time load startEnd.mat % load the preprocessed data instead figure ax = usamap(latlim, lonlim); setm(ax, 'FFaceColor', ocean) geoshow(states,'FaceColor',land) geoshow(startEnd.StartLat,startEnd.StartLon,'DisplayType','Point',... 'Marker','.','MarkerSize',5,'MarkerEdgeColor',[0 0 1]) title('Uber Weekend Rides - Starting Points') xlabel('San Francisco') textm(37.802069,-122.446618,'Marina') textm(37.808376,-122.426105,'Fishermans Wharf') textm(37.797322,-122.482409,'Presidio') textm(37.774546,-122.412329,'SOMA') textm(37.770731,-122.440481,'Haight') textm(37.818276,-122.498546,'Golden Gate Bridge') textm(37.819632,-122.376065,'Bay Bridge')

When you plot the longitude and latitude data, you just get messy point clouds and it is hard to see what's going on. Instead, I broke the map of San Francisco into rectangular blocks to approximate its districts. Here is the new plot of starting points by district.

dist = categories(startEnd.StartDist); cc = hsv(length(dist)); figure ax = usamap(latlim, lonlim); setm(ax, 'FFaceColor', ocean) geoshow(states,'FaceColor',land) for i = 1:length(dist) inDist = startEnd.StartDist == dist(i); geoshow(startEnd.StartLat(inDist),startEnd.StartLon(inDist),... 'DisplayType','Point','Marker','.','MarkerSize',5,'MarkerEdgeColor',cc(i,:)) end title('Uber Weekend Rides - Starting Points by District') xlabel('San Francisco') textm(37.802069,-122.446618,'Marina') textm(37.808376,-122.426105,'Fishermans Wharf') textm(37.797322,-122.482409,'Presidio') textm(37.774546,-122.412329,'SOMA') textm(37.770731,-122.440481,'Haight') textm(37.818276,-122.498546,'Golden Gate Bridge') textm(37.819632,-122.376065,'Bay Bridge')

This is a step in the right direction. Now that we have the starting and ending points grouped by districts, we can represent the rides as connections among different districts - this is essentially a graph with districts as nodes and rides as edges. To visualize this graph, we can use a popular social networking analysis tool Gephi, which was also used in another post, Analyzing Twitter with MATLAB.

You can export `StartDist` and `EndDist` as the edge list to Gephi in CSV format.

writetable(startEnd(:,{'StartDist','EndDist'}),'edgelist.csv',... 'WriteVariableNames',false)

Once you export the edge list, you can plot the connections (edges) between districts (nodes) in Gephi. Now it is much easier to see where people went during the weekend! To see a bigger image, check out the PDF version.

- The size of the district nodes represents their in-degrees, the number of incoming connections, and you can think of it as measure of popularity as destinations. SOMA, Haight, Mission District, Downtown, and The Castro are the popular locations based on this measure.
- The districts are colored based on their modularity, which basically means which cluster of nodes they belong to. It looks like people hang around set of districts that are nearby - SOMA, Downtown, Mission District are all located towards the south (green). The Castro, Haight, Western Addition in the center (purple) and it is strongly connected to Richmond and Sunset District in the west. Since those are residential areas, it seems people from those areas hang out in the other districts in the same cluster.
- The locals don't seem to go to Fisherman's Wharf or Chinatown in the north (red) very much - they are probably considered not cool because of tourists?

Now you know where to go in San Francisco during the weekend if you want to experience an active night life there. We just looked at the overall weekend data, but you can explore more by time slicing the data to see how the traffic pattern changes based on the time of day or day of the week. You may be able to find traffic congestions by calculating the speed using the timestamps. Try it yourself and share what you find here!.

Get
the MATLAB code

Published with MATLAB® R2014a

*The story is broken into two separate blog posts. This week’s installment is told by Marta Wilczkowiak. Marta is a senior application engineer at MathWorks UK, helping MATLAB users to design, develop and deploy production-quality code.*

*by Marta Wilczkowiak*

Technology can be the kiss of death to a conversation. But can we instead, through the thoughtful convergence of hardware, software, and real-time streaming data, actually encourage social interactions?

The barrier that existed between physical and virtual in engineering is disappearing in the world where internet companies buy robotics companies, fly drones, test self-driving cars and wearable devices, while the industrial players are moving into sensor networks and Internet of Things. While those new approaches are promising, such inter-disciplinary engineering collaboration is very difficult to achieve.

This is a story about our attempt to build such a convergent system in just three months, through cross-disciplinary collaboration. This blue-sky project is called Expo Conversations.

**EXPO Conversations: The Idea**

MATLAB EXPO is the pinnacle of MathWorks UK events. It is a unique opportunity for the MATLAB community to network with peers across a range of industries and academia. Planning starts early in the year and involves engineers from all technical groups at MathWorks. During one of those meetings, we set ourselves a challenge: to create a system bringing together all the new technologies we tell our customers about, including low cost hardware, data analytics and cloud computing. At the same time, we saw an opportunity to understand how ideas and influences flow throughout EXPO.

Our idea was to use low cost hardware to detect the interactions between the EXPO participants, and provide a web interface displaying the live analysis of the social activity. To encourage social interactions we would run a contest and give a prize to the most prolific networker at the event.

*Figure 1: Simulated social interactions in the EXPO Exhibition Hall*

This conversation took place in June. EXPO was scheduled for October 8th so we had little time to develop something from nothing!

To summarize, the objective was to:

- Demonstrate end-to-end technology stack integration, starting from low cost hardware to the cloud-based visualization
- Enable cross-disciplinary team collaboration
- Move from a blue sky idea to cloud deployment in 3 months

**July: First tag**

14 weeks to go

July started with discussion about the best approach to detect the attendee interactions: Cameras? Apps on mobile phones? RF tags? We wanted participants to remain anonymous so we opted for RF tags, with attendees opting-in simply by picking up a tag on their way into the exhibition hall. Conor, our 2nd year summer intern, with support from Mark (embedded systems) and Graham (communications) developed the first tags using Arduino Uno microcontrollers (ATmega328P) and RF12 433MHz transceivers. They showed that two early tag prototypes can communicate with each other and the received signal strength is somewhat related to the distance between tags. Characterizing this relationship remained a holy grail throughout the project! You can read more about this part of the project in Conor’s MakerZone blog post.

**August: Algorithms**

9.5 weeks to go

With promising early results, we kicked off the data analysis part of the project. We started by brainstorming the elements we would like to understand about the social interaction patterns. Some of our favorite suggestions: Can we quantify the relationship between caffeine intake and social activity? Can we detect wallflowers, conversation boosters… or conversation killers? Can we see if visitors to a demo station recommend it to others? We quickly realized that most of those questions needed a definition of a conversation. We decided that two tag holders are in a conversation if they stay near to each other for 30 seconds.

But how do we know whether two tags are close to each other? We decided to use Received Signal Strength Indication (RSSI) as a proxy for distance between two tags. On the simplicity-robustness scale, it is closest to simplicity and that’s what we needed given the project timeframe.

By detecting conversations between tags, we can learn the popularity of demo and catering stations (see Figure 2) by monitoring tags located at those places of interest. All the tag data would be transmitted to Raspberry Pi base stations, which would send it in batches to the analytics server every 10 seconds.

*Figure 2: EXPO Exhibition Hall floor plan with marked demos and catering stations (click to enlarge)*

We had 10 engineers helping on the data analysis challenge in August. However this posed its own challenges:

- Many short term contributors. Median time to spend on the project was a week including understanding the project, developing algorithms, integration with the system and testing.
- Simultaneous development. Regardless of which part of the system they would work on, people only had time in August, so we needed to find a way to work simultaneously rather than consecutively.
- Experience. In general, one’s time allocated to the project was inversely proportional to ones’ experience.
- Communication. The “low cost hardware” team did not understand the “analytics” team and vice versa.
- Data. We had no real data to test our algorithms on – the hardware was in development!

We needed to design a system that would allow numerous people with different backgrounds and levels of experience to:

- write algorithms simultaneously that would depend on each other
- integrate their algorithms with main system without disrupting other’s work
- test their algorithms without the data they relied upon

We put in place a 5-point plan to address these challenges:

- Strong core. Following the philosophy that a faulty design would hurt us more than a few bugs in the algorithmic code, we asked a couple of our Consultants to spend a few days designing the architecture and coding the core classes of our system.
- Small independent modules. Representing the system by small modules with clear responsibilities and interfaces was key to being able to absorb code contributions from a large team. See Figure 3 for the main system components. Creating sensible modules was facilitated by support in MATLAB for Object Oriented Programming. Keeping them small was possible thanks to the rich existing capability for reading, manipulating, analyzing and saving data. This includes built-in support for communications with databases, UDP, and different file formats as well as processing functions such as logical indexing, set operations, group statistics and graph analysis functions. All the modules were small: the largest file implementing a data processing algorithm contained 100 lines of code. Short code was obviously faster to write, but it was also faster to test.
- Self-service model repository. We designed the system as a series of “analysis agents”: classes tasked to perform a very specific analysis operation. Some agents would be responsible for cleaning entry data, some for tracking conversations; some would search for the most active networkers, etc. We combined Observer and Template Method patterns to create agents that had a relatively sophisticated way of interacting with the main infrastructure, while allowing team members with no Object Oriented experience to contribute code. In practice, this meant that the algorithm developer only had to write very simple procedural, code specific to her or his algorithm: how to get input data, process it and save output. Then with one line of code, she/he could set his agent to “observe” and act on changes on the database. As an example, on the appearance of new raw data an agent would clean it and save it back to the database.
- Data: fake it till you make it. What might be a disputable strategy in everyday social interactions was essential to enable the algorithm developers to write and test their code. Almost all parts of the analytics system were developed using simulated data. In the first few days of the project, in parallel with designing and implementing the analytics system, we created a high-level simulator that would create synthetic data representing a number of tags following journeys around the exhibition hall. Key to developing the first simulator in just a couple of days was that MATLAB offered all the required components in a single environment: random number generation, visualization and an expressive programming language. As we needed to switch to real data late in the project, we wrote our code accessing data (“Data Access Layer”) in such a way that it was possible to use data from the simulator, test logs and database without any change in the algorithmic code.
- Code generation. The next step was to add a layer to the simulator representing the logic of the system that would relay information received from the tags to the analytics server: how the tags synchronize to acquire and share the information, how they send this to base stations, and finally how base stations form and send packets to the analytics server. This “executable specification” enabled early detection of the differences between how the “hardware” and “analytics” team understood the requirements. But even more importantly once this logic has been expressed in MATLAB, we could use MATLAB Coder to automatically generate C code that would be used on the Raspberry Pi base stations. This saved us time and increased our confidence that the data we were using to test our algorithms was representative of what our data acquisition system would produce when ready. It also allowed us to test the embedded algorithms while the tag hardware was still being developed, since they were included in the simulator. This is a big benefit for embedded systems designers.

*Figure 3: High-level system architecture: uniform access to real and simulated data; analysis agents representing system logic and producing aggregate data about the event to be displayed in MATLAB or via web interface*

In the meantime, the hardware team was busy creating in reality what we had simulated: the networking system that would relay the tag data to a server.

By the end of August we had early versions of all the analysis agents, tested on data from the simulator. We also had a long list of known issues, including the fact that our prototype code could not keep up to the predicted data rate.

**September: “Production”**

5 weeks to go

September saw the transformation of our prototype system into an early version of what could be called a production system.

First, we harnessed our code in a series of unit tests (better late than never!) to allow faster iterations when fixing bugs and optimizing the code.

Second, we optimized the prototype implementations for speed. A couple of the biggest wins were:

Vectorization of code of calculation of the social network graph Implementation of some of the tag data filtering in SQL instead of MATLAB to reduce time overhead of data transfersBoth operations reduced the time needed to analyze a batch of data from 10s of seconds to fractions of seconds.

Third, we re-implemented our analysis engine to run each analysis agent on a separate thread. Out of many choices for doing so, we decided to recode the thin analysis engine skeleton in Java and deploy each analysis agent as a separate function on MATLAB Production Server. Thanks to using the loosely coupled agents in the first place, this required only changing one line in the agent’s code.

As a final step we used a machine on EC2 with an elastic IP to set up a Tomcat webserver with our webpage and a directory that would be regularly synchronized with our local directory containing all the results.

At the end of September we had a system analyzing the tags every 10 seconds and updating the webpage every minute but still working only with the simulator. Throughout September attempts to connect the two failed – due to an Endian disagreement, too slow a processing rate in the early version of the analysis engine, data overflow in the prototype data gathering stage, and more.

But we were not ready to cancel the project. In the end, we would not often have an opportunity to test it in a hall with 500 people, 8 catering stations, 9 demo stations and a BLOODHOUND SSC life-size model in the middle (yes, a big pile of metal in the middle of our RF experiment).

**October: EXPO or The First Integration Test**

1 week to go

First week of October has gone between integration attempts and bug fixing. You can see that the moods were mixed:

*Figure 4: Rainbow of emotions; Note as moods are improving when we move from bug chasing to handling some real hardware*

The day before EXPO (October 7th) three of us set off to Silverstone, the conference venue. Two hours after arriving everything was connected. The moment the webpage started to show the two test tags in “conversation” with each other was worth all those sleepless nights. We stared in awe. By 7pm all 8 base stations were connected and receiving. At that point we put the system to sleep and left the site.

The morning of EXPO passed in a flash. We started with a pile of 200 tags (see Figure 5). By 9.30 am all of them has been distributed. Attendees were curious about the system and wanted to participate. The webpage started to update with plausible results (see sample in Figure 7). All team members, most of whom so far have seen only a small part of the system, were hypnotized by the updating web page. Throughout the day the demo drew a big crowd who asked about everything from working with low cost hardware, through machine learning and “big data” to web deployment. We also highlighted parts of the system during all relevant talks: Machine Learning, Accelerating MATLAB algorithms, Deployment, Modelling, Simulation, and Real-Time Testing for Model-Based Design. Finally, in the wrap up talk we announced the tag number that won the networking competition prize. We knew the system was working when we discovered that the winner was our ex-colleague Michael, known for being a social butterfly.

So through this project we achieved:

- end-to-end technology stack integration, starting from low cost hardware to the cloud-based visualization
- cross-disciplinary team collaboration
- blue sky idea to cloud deployment in 3 months

But as EXPO was also the first time the system ran together, we still had a lot of work to understand what it was really doing.

Learn more as the story continues: Expo Conversations – Part 2

*Figure 5 Tags waiting to be handed to attendees*

*Figure 6: Queue to pick up the tags*

*Figure 7: Screenshot of the live social networking analysis*

We all use Fourier analysis every day without even knowing it. Cell phones, disc drives, DVDs, and JPEGs all involve fast finite Fourier transforms. This post, which describes touch-tone telephone dialing, is the first of three posts about the computation and interpretation of FFTs. The posts are adapted from chapter 8 of my book, *Numerical Computing with MATLAB* .

Touch-tone telephone dialing is an example of everyday use of Fourier analysis. The basis for touch-tone dialing is the Dual Tone Multi-Frequency (DTMF) system. *Synthesis* is the generation of analog tones to represent digits in phone numbers. *Analysis* is the decoding of these tones to retrieve the digits. The program `touchtone`, which is available here, or which is included with my NCM App, demonstrates both these aspects of DTMF.

The telephone dialing pad acts as a 4-by-3 matrix. Associated with each row and column is a frequency. These basic frequencies are

fr = [697 770 852 941]; fc = [1209 1336 1477];

These frequencies were chosen so that none is an integer multiple of any other. In fact the ratios are 21/19. This avoids interfering overtones.

If `s` is a character that labels one of the buttons on the keypad, the corresponding row index `k` and column index `j` can be found with

switch s case '*', k = 4; j = 1; case '0', k = 4; j = 2; case '#', k = 4; j = 3; otherwise, d = s-'0'; j = mod(d-1,3)+1; k = (d-j)/3+1; end

An important parameter in digital sound is the sampling rate.

Fs = 32768

A vector of points in one-fourth of a second at this sampling rate is

t = 0:1/Fs:0.25

The tone generated by the button in position `(k,j)` is obtained by superimposing the two fundamental tones with frequencies `fr(k)` and `fc(j)`.

y1 = sin(2*pi*fr(k)*t); y2 = sin(2*pi*fc(j)*t); y = (y1 + y2)/2;

If your computer is equipped with a sound card, the statement

sound(y,Fs)

plays the tone.

For example, this is the display produced by `touchtone` for the "1" button. The top subplot depicts the two underlying frequencies and the bottom subplot shows a portion of the signal obtained by averaging the sine waves with those frequencies.

How is it possible to determine a phone number by listening to the signal generated? This involves the FFT, the Finite Fourier Transform, or, more likely, some specialized version of the FFT adapted to this task. I will discuss the FFT in my next post.

The data file `touchtone.mat` contains a recording of a telephone being dialed. The file is available from

http://www.mathworks.com/moler/ncm/touchtone.mat

Go to the `Home` tab above the MATLAB command window, click on the `Open` tab and insert this URL into the `File name` slot. This will load a structure `y` into your MATLAB workspace. Then the statement

y

will produce

y = sig: [1x74800 int8] fs: 8192

This shows that `y` has two fields, an integer vector `y.sig`, of length 74800, containing the signal, and a scalar `y.fs`, with the value 8192, which is the sample rate.

max(abs(y.sig))

reveals that the elements of the signal are bounded in absolute value by 127. So the statements

Fs = y.fs; y = double(y.sig)/128;

save the sample rate, rescale the vector and convert it to double precision. The statements

n = length(y); t = (0:n-1)/Fs

reproduce the sample times of the recording. The last component of `t` is `9.1307`, indicating that the recording lasts a little over 9 seconds. The next figure is a plot of the entire signal.

This recording is noisy. You can even see small spikes on the graph at the times the buttons were clicked. It is easy to see that 11 digits were dialed, but, on this scale, it is impossible to determine the specific digits.

Here is the magnitude of the FFT of the signal, which is the key to determining the individual digits.

p = abs(fft(y)); f = (0:n-1)*(Fs/n); plot(f,p); axis([500 1700 0 600])

The *x* -axis corresponds to frequency. The `axis` settings limit the display to the range of the DTMF frequencies. There are seven peaks, corresponding to the seven basic frequencies. This overall FFT shows that all seven frequencies are present someplace in the signal, but it does not help determine the individual digits.

The `touchtone` program also lets you break the signal into 11 equal segments and analyze each segment separately. The next figure is the display from the first segment.

For this segment, there are only two peaks, indicating that only two of the basic frequencies are present in this portion of the signal. These two frequencies come from the "1" button. You can also see that the waveform of a short portion of the first segment is similar to the waveform that our synthesizer produces for the "1" button. So we can conclude that the number being dialed in `touchtone` starts with a 1. We leave it as an exercise to continue the analysis with `touchtone` and identify the complete phone number.

`touchtone` has a hidden "Easter egg" that, as far as I know, has never been uncovered. If anyone finds it, submit a comment and let us know.

Cleve Moler, *Numerical Computing with MATLAB*, Electronic Edition, MathWorks, <http://www.mathworks.com/moler/index_ncm.html>,

Print Edition, SIAM Revised Reprint, SIAM, 2008, 334 pp., <http://bookstore.siam.org/ot87> .

Get
the MATLAB code

Published with MATLAB® R2014a

But what about Simulink? A lot of work is done using the mouse, so the touch typing configuration is less relevant. Are there guidelines or tips to build, edit, and navigate simulink models as efficiently as possible?

This week I share mine!

**Navigating**

When navigating, I typically keep my thumb on the **Alt** key, my middle finger on the **Esc** key, and the index on the **1** key.

This allows me to easily do:

**Esc**to go up one level with the middle finger**Atl+1**to zoom 100%, using thumb on**Alt**and index on**1****Space bar**to zoom to fit with the thumb**Alt+Tab**to switch between windows, typically the Library Browser and the model, using thumb on**Alt**and index on**Tab**.

Here is a typical example of me exploring a model. This typically ends up being a pattern of hitting the **Space bar** when I enter a new subsystem to get the big picture, then **Atl+1** to see the blocks zoomed at 100% and **Esc** to go up one level when I am done.

**Editing**

When editing, what I need the most are the **Ctrl** and **Ctrl+Shift** shortcuts. So I slightly move my left hand to keep my small finger on the **Ctrl** key, and my ring finger on the **Shift** key.

Here is the list of shortcuts I typically use in this configuration:

Bottom line:

**Ctrl+Z**to Undo**Ctrl+X**to Cut**Ctrl+C**to Copy**Ctrl+V**to Paste**Ctrl+B**to Build the model using Simulink Coder

Middle Line:

**Ctrl+A**to Select all blocks in the current system**Ctrl+S**to Save**Ctrl+D**to update the diagram (ALWAYS keep your model as close as possible to an updatable state, this will save you lots of debugging time!)**Ctrl+G**to Group selected blocks in one subsystem**Ctrl+H**to open the Model Explorer (am I the only one to pronounce it Model HHHexplorer because of that?)

Top line:

**Ctrl+E**to open the model configuration**Ctrl+R**to Rotate blocks**Ctrl+T**to start a simulation**Ctrl+Y**to redo something previously undone**Ctrl+U**to open a masked subsystem. With the introduction of the arrow badge, I use this shortcut less often since R2012b**Ctrl+I**to flip a block left-right. Ok, this one requires some intense stretching to be done with only one hand... be careful to not hurt yourself.

In addition go the **Ctrl** shortcuts, this hand position allows you to access a set of **Ctrl+Shift** shortcuts, by using the little finger for **Ctrl**, and the ring finger for **Shift**:

**Ctrl+Shift+X**to Comment out and uncomment blocks**Ctrl+Shift+Y**to Comment through blocks**Ctrl+Shift+H**to remove highlighting in the model**Ctrl+Shift+G**to Expand a subsystem**Ctrl+Shift+L**to open the Simulink Library Browser. (This one cannot be done only with the left hand, but I thought I should point out his existence. Most of the time, when I edit model I alternate between the Library Browser and the model using**Atl+Tab**)

You can see the full list of Simulink shortcuts here.

Since many **Ctrl+Shift** combinations do not have shortcuts assigned, I like to use those to assign custom shortcuts using an sl_customization file, as I previously highlighted in this post. One example for me is **Ctrl+Shift+Z** to display or not the name of blocks.

Here is an animation showing usage of the most common shortcuts:

**Now it's your turn**

Share your tricks to use Simulink as efficiently as possible by leaving a comment here.

]]>

I’d like to introduce you to this week’s guest blogger, Graham Dudgeon. Graham is with our Industry Marketing Team at MathWorks, and focuses on the Utilities & Energy Industry. In this blog, Graham shares a story about age being no barrier to exploring ideas and concepts in MATLAB and Simulink.

Hi Everyone, and thank you Loren for inviting me to be a guest blogger. This story begins with my dog, waiting patiently at the patio door expecting either of my two sons to let her in. She was out of luck. I am no dog psychologist, but I figure that because she could see the boys very clearly, that they in turn could see her. In her mind, she didn’t need to bark to raise their attention…. and so on she waited while my sons continued watching Phineas and Ferb. Once I was alerted to her predicament a mere fraction of a second after entering the family room (by seeing her hopeful eyes and her tail wagging with a vigor in direct proportion to how long she had been waiting), I let her in and then reinforced to my sons that they need to be more aware of their surroundings. This is where MATLAB and Simulink come into the story… To help pooch out, I made a slight modification to demoimaqsl_rgbhistogram_win.slx , aimed the web cam at the patio door and when pooch came to the door, we were alerted to here presence by a ‘gong’ that I took from ‘gong.mat’, and had a MATLAB callback to execute the sound when a color threshold was reached – it turns out red works best for a fox hound. Was pooch consequently let into the house by my boys? Unfortunately not. The boys decided that a much better use of my endeavor was to see if they could creep past the field of view without the gong going off. This kept them entertained for longer than I expected.

My eldest son, who is 9 years old, started asking questions about how the model worked, how the camera worked, what was doing what, how do you detect pooch, etc. etc. I must have answered with sufficient clarity and enthusiasm as he decided to take what I had done and modify it for his Science Expo.

My son watched the RGB plot as we passed different colors in front, and gained a feel for how to develop a detection algorithm. Our process was to let my son watch only the RGB plot, not the object in front of the camera, and tell me what color he thought was there. ‘Red!’, ’Blue!’, ’Kinda red and blue, and a bit of green!’. We would then confirm his observations by looking at the actual object and moving it again across the field of view to see the RGB response change with the movement. My son associated y-axis magnitude with detection and figured out that to detect, he need simply tell the computer to ‘gong’ when a color got to a certain magnitude. My son also figured that you could remain ‘invisible’ if the detector was tuned only for red and you wore another color, and that a black background would cause less noise on the color signals, as black caused only a small RGB signal. My youngest son also helped out, by wearing different colored shirts and walking (sometimes running) past the web cam. The list of insights my son gained go on and on. He presented his work at the Science Expo, and was met with great enthusiasm and interest by the High School Scientists who assessed his work.

So how does this story end? A few days ago my eldest son told me that when he grows up he wants to be an engineer… and that he wants to work for MathWorks. On hearing this, my 5 year old son said he wanted to be a scientist (as he ‘wants to see into stuff’)… and that he too wants to work for MathWorks. So the story continues, and our future is in good hands

And what of pooch? She’s standing outside the patio door (with hopeful eyes and tail wagging vigorously) waiting for me to finish this post.

Before I sign off, let me give you my main takeaway from this story… You can use MathWorks tools to help kids of any age explore and evolve their curiosity, by connecting them directly to the world of Science and Engineering. Give it a try!

Get
the MATLAB code

Published with MATLAB® R2014a

Since the video at the center of this post is about a rap battle between MATLAB and Simulink, I thought it would be interesting to bring that to a higher level and see how the video processing could have been implemented in Simulink instead of MATLAB.

Seriously... who wants to write code when you can simply connect blocks?

**Overview**

Since I have absolutely zero knowledge in image processing, I simply re-implemented the algorithm suggested by Matt in Simulink. The resulting model looks like:

Let's look at the details.

**Importing the video**

To begin, I used the From Multimedia File from the DSP System Toolbox to directly read the AVI-file with our rappers in front of a green background. Then I use the resize block from the Computer Vision System Toolbox. This allows me to work on images of a relatively small size. After, I normalize the data, as suggested by Matt.

**The Greenness factor**

I really like the simple equation the Matt figured out to calculate the "greenness" of each pixel in the image. In MATLAB, it looks like:

```
% Greenness = G*(G-R)*(G-B)
greenness = yd(:,:,2).*(yd(:,:,2)-yd(:,:,1)).*(yd(:,:,2)-yd(:,:,3));
```

**Thresholding**

The next piece of code is something that many users struggle to figure out:

`thresh = 0.3*mean(greenness(greenness>0));`

To implement logical indexing (e.g. `A(A>0)`) in Simulink, we need to combine the Find block with a Selector block. This results in a variable-size signal

**Combining the images**

The last line we need to implement is to replace the pixels identified as green in the foreground image but their values in the background image. In MATLAB, this is accomplished using this line:

`rgb2(isgreen) = rgb1(isgreen);`

In Simulink, we need to extend the logic seen at the previous step. We first select the desired pixels from the background image, and then assign them to the foreground image:

Note that the Computer Vision System Toolbox has a Compositing block implementing this functionality.

Here is what the final result looks like:

**Now it's your turn**

What do you think? Is it more convenient to implement this algorithm in MATLAB or Simulink? How do you typically decide between a *code-based* or a *block-based* implementation? Let us know by leaving a comment here.

For several years we thought Hadamard matrices showed maximum element growth for Gaussian elimination with complete pivoting. We were wrong.

I want to continue the discussion from my previous blog post of element growth during Gaussian elimination. Suppose we are solving a system of linear equations involving a matrix $A$ of order $n$. Let $a_{i,j}^{(k)}$ denote the elements in the matrix after the $k$ th step of the elimination. Recall that the *growth factor* for the pivoting process we are using is the quantity

$$ \rho_n = { \max_{i,j,k} |a_{i,j}^{(k)}| \over \max_{i,j} |a_{i,j}| } $$

Complete pivoting is intended to control the growth factor by using both row and column interchanges at each step in the reduction to bring the largest element in the entire unreduced matrix into the pivot position. In his 1962 analysis of roundoff error in Gaussian elimination, J. H. Wilkinson proved that the growth factor for complete pivoting satisfies

$$ \rho_n \le g(n) $$

where the *growth function* is

$$ g(n) = (n \ 2 \ 3^{1/2} 4^{1/3} \cdots n^{1/(n-1)})^{1/2} \approx 1.8 \ n^{1/2} n^{1/4 \log{n}} $$

Wilkinson's analysis makes it clear that the growth can never actually be anywhere near this large.

We saw that the growth factor for partial pivoting is $2^{n-1}$. For complete pivoting, it's much, much less. To get a feeling for now much, express $g(n)$ as a MATLAB one-liner

g = @(n) sqrt(n*prod((2:n).^(1./(1:(n-1)))))

g = @(n)sqrt(n*prod((2:n).^(1./(1:(n-1)))))

Check out $n = 100$

g(100)

ans = 3.5703e+03

So here $g(n) \approx 35n$.

Jacques Hadamard was a French mathematician who lived from 1865 until 1963. He made contributions in many fields, ranging from number theory to partial differential equations. The matrices named after him have entries +1 and -1 and mutually orthogonal rows and columns. The $n$ -dimensional parallelotope spanned by the rows of a Hadamard matrix has the largest possible volume among all parallelotopes spanned by matrices with entries bounded by one. We are interested in Hadamard matrices in the MATLAB world because they are the basis for Hadamard transforms, which are closely related to Fourier transforms. They are also applicable to error correcting codes and in statistics to estimation of variance.

So MATLAB already has a `hadamard` function. Here is the 8-by-8 Hadamard matrix.

H = hadamard(8)

H = 1 1 1 1 1 1 1 1 1 -1 1 -1 1 -1 1 -1 1 1 -1 -1 1 1 -1 -1 1 -1 -1 1 1 -1 -1 1 1 1 1 1 -1 -1 -1 -1 1 -1 1 -1 -1 1 -1 1 1 1 -1 -1 -1 -1 1 1 1 -1 -1 1 -1 1 1 -1

Check that the rows are perpendicular to each other.

H'*H

ans = 8 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 8

The orthogonality of the rows leads to element growth during Gaussian elimination. If we call for the LU factorization of `H`, no pivoting actually takes places, but the same result would be produced by complete pivoting that settles ties in favor of the incumbent.

[L,U] = lu(H)

L = 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 1 1 1 0 0 0 0 1 0 0 0 1 0 0 0 1 1 0 0 1 1 0 0 1 0 1 0 1 0 1 0 1 1 1 1 1 1 1 1 U = 1 1 1 1 1 1 1 1 0 -2 0 -2 0 -2 0 -2 0 0 -2 -2 0 0 -2 -2 0 0 0 4 0 0 0 4 0 0 0 0 -2 -2 -2 -2 0 0 0 0 0 4 0 4 0 0 0 0 0 0 4 4 0 0 0 0 0 0 0 -8

This, then, is one matrix at $n = 8$ where

$$ \rho_n = n $$

For almost 30 years everybody believed that Hadamard matrices were the extreme cases and that the growth factor $\rho_n$ for complete pivoting with real matrices was equal to $n$.

At various times in the 1960s and '70s I personally did numerical experiments looking for matrices with large pivot growth. The computers I had at the time limited the order of my matrices to a few dozen. I never found any cases of $\rho_n > n$.

One of the people who worked on the problem was Leonard Tornheim at Chevron Research in California. He also did a lot of numerical experiments. He found a complex 3-by-3 matrix with $\rho_3 > 3$. He also proved that for real matrices, $\rho_3 = 2 \ 1/4$.

In 1968, Colin Cryer made an official conjecture that $\rho_n$ was equal to $n$ for real matrices.

A 1988 paper by Jane Day and Brian Peterson provides a good survey of what was known up to that time about element growth with complete pivoting.

Then, in 1991, Nick Gould, from Oxford, surprised us all by announcing the discovery of real matrices for which the complete pivoting $\rho_n$ is greater than $n$. For an $n$ -by- $n$ matrix, with $n$ between 13 and 16, Nick set up a large, sparse, nonlinear optimization problem involving roughly $n^3/3$ variables, and tackled it with the LANCELOT package that he and his colleagues Andy Conn and Philippe Toint had developed. Most of his paper is about a 13-by-13 matrix with

$$ \rho_{13} = 13.0205 $$

This just barely exceeds the $\rho_n = n$ conjecture. But he also found some slightly larger matrice that do a better job.

$$ \rho_{14} = 14.5949 $$

$$ \rho_{15} = 16.1078 $$

$$ \rho_{16} = 18.0596 $$

The 16-by-16 is particularly dramatic because it is not a Hadamard matrix.

So, as far as I know, this is where the matter still stands today. Nobody is interested in running larger experiments. I have no idea what $\rho_{100}$ or $\rho_{1000}$ might be. Of course, we don't really need to know because we don't use complete pivoting in practice.

This reminds me of a story that I like to tell. I'm not sure it's actually true.

Peter Businger was a grad student at Stanford in Math and Computer Science at the same time I was in the early '60s. He worked with Gene Golub to write the first published Algol procedure for the QR matrix factorization using Householder reflections. Peter left Stanford to join Joe Traub's group at Bell Labs that was producing what they called the Port library of mathematical software. They were importing programs from everybody else, running the routines on Bell Labs' machines, thoroughly testing the software, deciding which codes were best, and producing a consolidated library.

Peter was in charge of matrix computation. He checked out all the imported linear equation solvers, including mine. He decided none of them was satisfactory, because of this pivoting growth question. So he wrote a program of his own that did partial pivoting, and monitor the growth. If the growth was too large, the program switched to complete pivoting.

At the time, Bell Labs was trying to patent software, so Peter named his new subroutine "P^8". This stood for "Peter's Pseudo Perfect Perfect Partial Pivot Picker, Patent Pending". In my opinion, that is one of the great all time subroutine names.

The experience inspired Peter to go to law school at night and change jobs. He switched to the Bell Labs patent department.

I've lost contact with Peter, so I can't ask him if he really did name the routine P^8. If he didn't, he should have.

Les Foster, who I mentioned in the last post for finding examples of exponential growth with partial pivoting, promotes what he calls "Rook Pivoting". Search down the column for the largest element, as in partial pivoting, but then also search along that row for the last element. This is intermediate between partial and complete pivoting. Les is able to prove some results about pivot growth. But the algorithm does not generalize well to a block form.

I was present at a milestone in the history of Hadamard matrices. The orthogonality conditions readily imply that the order $n$ of a Hadamard matrix must be 1, 2, or a multiple of 4. But the question remains, does a Hadamard matrix of order $n = 4k$ exist for every $k$? This is an open question in number theory.

It is fairly easy to create Hadamard matrices of some sizes. A nifty way to generate a Hadamard matrix of order a power of two is the recursion.

H = 1; for n = 2, 4, 8, ... H = [H H H -H]; end

If `A` and `B` are Hadamard matrices, then so is

kron(A,B)

By 1961, with these and other similar constructions, it turned out that it was known how to construct Hadamard matrices for all orders $n \le 100$ that are multiples of four except $n = 92$. In 1961 I had a summer job at JPL, Caltech's Jet Propulsion Lab. My office was in a temporary trailer and my trailer mate was Len Baumert. Len proudly showed me a color graphic that he had just made and that he was proposing for the cover of *Scientific American*. It was a 92-by-92 matrix made up from 26-by-26 blocks of alternating light and dark cells representing +1 and -1. The graphic didn't make the cover of *Scientific American*, but I kept my copy for a long time.

Len had done a computation on JPL's machines that filled in the last value of $n$ less than 100. His paper with his colleagues Sol Golomb, a professor at USC, and Marshall Hall, Jr., a professor at Caltech, was published in the prestigious Bulletin of the AMS. It turns out that I had taken a course on difference sets, the mathematics involved in generating the matrix, from Hall at Caltech.

Here is a MATLAB picture of `Baumert92`. You can get the function from this link.

H = Baumert92; pcolor92(H);

Let's check that we get the expected pivot growth.

[L,U] = lu(H); unn = U(end,end)

unn = -92

This reference to the book by Trefethen and Bau should have been included in my previous post about partial pivoting. Lecture 22 has an explanation the stability of partial pivoting in practice.

Lloyd N. Trefethen and David Bau, III, *Numerical Linear Algebra*, <http://bookstore.siam.org/ot50>, SIAM, 1997, 362pp.

Peter Businger, "Monitoring the numerical stability of Gaussian elimination", <http://link.springer.com/article/10.1007/BF02165006>, Numerische Mathematik 16, 4, 1971, 360-361.

Jane Day and Brian Peterson, "Growth in Gaussian elimination", <http://www.jstor.org/stable/2322755>, American Mathematical Monthly, 95, 6, 1988, 489-513.

Nick Gould, "On growth in Gaussian elimination with complete pivoting", <http://www.numerical.rl.ac.uk/people/nimg/pubs/Goul91_simax.pdf> SIAM Journal on Matrix Analysis and Applications 12, 1991, 351-364.

Leslie Foster, "The growth factor and efficiency of Gaussian elimination with rook pivoting", <http://www.math.sjsu.edu/~foster/gerpwithcor.pdf>, Journal of Computational and Applied Mathematics, 86, 1997, 177-194.

Leslie Foster, "LURP: Gaussian elimination with rook pivoting", matlabcentral/fileexchange/37721-rook-pivoting, MATLAB Central File Exchange, 2012.

Leonard Baumert, S. W. Golomb, and Marshall Hall, Jr, "Discovery of an Hadamard Matrix of Order 92", <http://dx.doi.org/10.1090%2FS0002-9904-1962-10761-7>, Bulletin of the American Mathematical Society 68 (1962), 237-238.

Get
the MATLAB code

Published with MATLAB® R2014a

Here is my answer.

]]>

Have you ever needed to move a figure? Perhaps because it is offscreen>? Or at least the menus or close box are not visible? This happens to me from time to time when I get code from someone else who is programming specifically to the size of their display. Working on a laptop, I often have fewer available pixels. What to do?

I used the Function Browser to search for "move figure" and found `movegui`. It might be just the ticket for you.

As you can see from the documentation, `movegui` gives you a fair amount of latitude for moving your figure around. The two choices I find most useful and easiest to remember for my purposes: `onscreen` and `center`. I think you can guess what they each do :-)

Do you use the Function Browser to help you? Do you ever need to reposition "lost" figure windows, like I have? Let me know here.

Get
the MATLAB code

Published with MATLAB® R2014a

Our latest update has two features that should be of interest to blog readers:

- Searching with Categories: Blog authors can tag their posts with Categories; keywords, or phrases that organize or describe a particular post. It’s a helpful tool for finding related content. Like a particular post? Just click the Category link at the bottom to find similar posts.
- Print a Post: Several blog readers requested this feature. Now you can easily print out posts by your favorite bloggers and read them when you’re off the grid. After spending all day staring into a monitor, there’s nothing better than being outside, reading a nicely formatted post on paper.Bliss.

We’re always looking to improve the blogs. Tell us what you think about these new features or let us know what we can do to make the blogs section better for you. While you’re thinking about all this, drop us a line in the Comments section below.

Thanks!

]]>In rare cases, Gaussian elimination with partial pivoting is unstable. But the situations are so unlikely that we continue to use the algorithm as the foundation for our matrix computations.

I almost hesitate to bring this up. Gaussian elimination with partial pivoting is potentially unstable. We know of a particular test matrix, and have known about it for years, where the solution to simultaneous linear equations computed by our iconic backslash operator is less accurate than we typically expect. The solution is contaminated by unacceptably large roundoff errors associated with large elements encountered during the elimination process. The offending matrix is generated by the following function, which is a special case of the `gfpp` function in Nick Higham's collection of MATLAB test matrices.

```
type gfpp
```

function A = gfpp(n) % A = gfpp(n) Pivot growth of partial pivoting. A = eye(n,n) - tril(ones(n,n),-1); A(:,n) = 1;

Here is the 8-by-8.

n = 8 A = gfpp(n)

n = 8 A = 1 0 0 0 0 0 0 1 -1 1 0 0 0 0 0 1 -1 -1 1 0 0 0 0 1 -1 -1 -1 1 0 0 0 1 -1 -1 -1 -1 1 0 0 1 -1 -1 -1 -1 -1 1 0 1 -1 -1 -1 -1 -1 -1 1 1 -1 -1 -1 -1 -1 -1 -1 1

Gaussian elimination with partial pivoting does not actually do any pivoting with this particular matrix. The first row is added to each of the other rows to introduce zeroes in the first column. This produces twos in the last column. As similar steps are repeated to create an upper triangular `U`, elements in the last column double with each step. The element growth in the final triangular factor `U` is `2^(n-1)`.

[L,U] = lu(A); U

U = 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 2 0 0 1 0 0 0 0 4 0 0 0 1 0 0 0 8 0 0 0 0 1 0 0 16 0 0 0 0 0 1 0 32 0 0 0 0 0 0 1 64 0 0 0 0 0 0 0 128

So if the matrix order `n` is set to the number of bits in the floating point word (plus one), the last diagonal element of `U` will grow to be `1/eps`.

n = 53 A = gfpp(n); [L,U] = lu(A); unn = U(n,n)

n = 53 unn = 4.5036e+15

If we try to use backslash to solve a linear system involving this matrix and a random right hand side, the back substitution will begin with `U(n,n)`. The roundoff errors can be expected to be as large as `eps(U(n,n))` and swamp the solution itself.

b = randn(n,1); x = A\b;

Ordinarily, we expect the solution computed by Gaussian elimination to have a residual on the order of roundoff error. But in this case, the residual is many orders of magnitude larger. Backslash has failed.

r = b - A*x; relative_residual = norm(r)/(norm(A)*norm(x))

relative_residual = 3.7579e-05

Anything we do to this matrix to provoke pivoting during the elimination will prevent the growth of the elements in `U`. One possibility is to swap the first two rows.

A = gfpp(n); A([1 2],:) = A([2 1],:); [L,U] = lu(A); unn = U(n,n)

unn = 2

That's much better. Now use backslash and look at the residual.

x = A\b; r = b - A*x; relative_residual = norm(r)/(norm(A)*norm(x))

relative_residual = 1.7761e-17

So the residual is now on the order of roundoff relative to `A` and `x`. Backslash works as expected.

Another way to provoke pivoting is to perturb the matrix with a little white noise.

A = gfpp(n); A = A + randn(n,n); [L,U] = lu(A); unn = U(n,n)

unn = -3.8699

So `U` is well behaved. Let's see about backslash.

x = A\b; r = b - A*x; relative_residual = norm(r)/(norm(A)*norm(x))

relative_residual = 1.2870e-16

Again the residual is on the order of roundoff and backslash works as expected.

Suppose we are using Gaussian elimination with any pivoting process to solve a system of linear equations involving a matrix $A$. Let $a_{i,j}^{(k)}$ denote the elements in the matrix after the $k$ th step of the elimination. The *growth factor* for that pivoting process is the quantity

$$ \rho = { \max_{i,j,k} |a_{i,j}^{(k)}| \over \max_{i,j} |a_{i,j}| } $$

This growth factor is a crucial quantity in Wilkinson's roundoff error analysis. If $\rho$ is not too large, then the elimination algorithm is stable. Unfortunately, the matrix we have just seen, `gfpp`, shows that the growth factor for partial pivoting is at least $2^{n-1}$.

At the time he was presenting his analysis in early 1970s, Wilkinson knew about this matrix. But he said it was a very special case. He reported:

It is our experience that any substantial increase in size of elements of successive An is extremely uncommon even with partial pivoting. ... No example which has arisen naturally has in my experience given an increase as large as 16.

During the development of LINPACK in the late 1970s, a colleague at Argonne, J. T. Goodman, and I did some experiments with random matrices to check for element growth with partial pivoting. We used several different element distributions and what were then matrices of modest order, namely values of $n$ between 10 and 50. The largest growth factor we found was for a matrix of 0's and 1's of order 40. The value was $\rho$ = 23. We never saw growth as large as linear, $\rho = n$, let alone exponential, $\rho = 2^{n-1}$.

Nick Trefethen and Rob Schreiber published an important paper, "Average-case Stability of Gaussian Elimination", in 1990. They present both theoretical and probabilistic models, and the results of extensive experiments. They show that for many different types of matrices, the typical growth factor for partial pivoting does not increase exponentially with matrix order $n$. It does not even grow linearly. It actually grows like $n^{2/3}$. Their final section, "Conclusions", begins:

Is Gaussian elimination with partial pivoting stable on average? Everything we know on the subject indicates that the answer is emphatically yes, and that one needs no hypotheses beyond statistical properties to account for the success of this algorithm during nearly half a century of digital computation.

The Higham Brothers, Nick and Des Higham, published a paper in 1989 about growth factors in Gaussian elimination. Most of their paper applies to any kind of pivoting and involves growth that is $O(n)$. But they do have one theorem that applies only to partial pivoting and that characterizes all matrices that have $\rho = 2^{n-1}$. Essentially these matrices have to have the behavior we have seen in `gfpp`, namely no actual pivoting and each step adds a row to all the later rows, doubling the values in the last column in the process.

In two different papers, Stephen Wright in 1993 and Les Foster in 1994, the authors present classes of matrices encountered in actual computational practice for which Gaussian elimination with partial pivoting blows up with exponential growth in the elements in the last few columns of `U`. Wright's application is a two-point boundary value problem for an ordinary differential equation. The discrete approximation produces a band matrix. Foster's application is a Volterra integral equation model for population dynamics. The discrete approximation produces a lower trapezoidal matrix. In both cases the matrices are augmented by several final columns. Certain choices of parameters causes the elimination to proceed with no pivoting and results in exponential in growth in the final columns. The behavior is similar to `gfpp` and the worst case growth in the Highams theorem.

Complete pivoting requires $O(n^3)$ comparisons instead of the $O(n^2)$ required by partial pivoting. More importantly, the repeated access to the entire matrix complete changes the memory access patterns. With modern computer storage hierarchy, including cache and virtual memory, this is an all-important consideration.

The growth factor $\rho$ with complete pivoting is interesting and I plan to investigate it in my next blog.

`lugui` is one of the demonstration programs included with Numerical Computing with MATLAB. I will also describe it in my next blog. For now, I will just use it to generate a graphic for this blog. Here is the result of a factorization, the lower triangular factor `L` in green and the upper triangular factor `U` in blue.

```
A = gfpp(6);
lugui(A,'partial')
```

Nicholas J. Higham and Desmond J. Higham, "Large Growth Factors in Gaussian Elimination with Pivoting", <http://www.maths.manchester.ac.uk/~higham/papers/hihi89.pdf>, SIAM J. Matrix Anal. Appl. 10, 155-164, 1989.

Lloyd N. Trefethen and Robert S. Schreiber, "Average-case Stability of Gaussian Elimination", https://courses.engr.illinois.edu/cs591mh/trefethen/trefethen_schreiber.pdf, SIAM J. Matrix Anal. Appl. 11, 335-360, 1990.

Stephen J. Wright, A Collection of Problems for Which Gaussian "Elimination with Partial Pivoting is Unstable", <ftp://info.mcs.anl.gov/pub/tech_reports/reports/P266.ps.Z>, SIAM J. Sci. Statist. Comput. 14, 231-238, 1993.

Leslie Foster, "Gaussian Elimination Can Fail In Practice", <http://www.math.sjsu.edu/~foster/geppfail.pdf>, SIAM J. Matrix Anal. Appl. 15, 1354-1362, 1994.

Get
the MATLAB code

Published with MATLAB® R2014a

Photo by Sara Kerens

Nick is what you might call a “heroic user” of MATLAB. He was there at the first class in which Cleve introduced the original FORTRAN MATLAB as a teaching tool. Near the lobby of the main MathWorks building here in Natick, we have a framed copy of the very first purchase order for MATLAB. Who bought that first copy of MATLAB? None other than Nick Trefethen. And it’s not merely a question of precedence or longevity. Across his long career, Nick has been a contributor of the highest order to the ecosystem of mathematical code written in MATLAB. He is the author of many books, an influential educator, and most recently he has been leading an innovative open source software project, a project that grew out of work he initiated back in 2002: Chebfun.

I can personally attest that the examples associated with Chebfun are among the most beautiful and thoughtful you will find documenting any mathematical project on the web. We recently conducted a virtual interview with Nick to learn about the latest release of Chebfun. Here’s what we learned.

Q: When did you first start using MATLAB?

A: I took a course from Cleve Moler when I was a graduate student at Stanford, so I must have typed x = A\b for the first time around 1978. When MathWorks started selling MATLAB, I bought the very first copy! That was in February, 1985. It was about ten years later when MathWorks told me I had been the first.

Q: What is Chebfun?

A: Chebfun is a system written in MATLAB that looks and feels like MATLAB, except that vectors and matrices become functions and operators. Whenever MATLAB has a discrete command, Chebfun aims to have the same command overloaded to execute a continuous analogue. For example, SUM(f) gives you the integral of a chebfun f over its interval of definition, MAX(f) gives you its maximum, ROOTS computes the roots, and so on.

Q: Why is it called Chebfun?

A: Well, Chebfun stands for Chebyshev-functions. Our algorithms are based on representations of functions by piecewise Chebyshev polynomial interpolants. The details are given in my book Approximation Theory and Approximation Practice.

Q: This Chebyshev fellow must have been an interesting character. Can you tell us a little more about how his name came to be associated with your favorite functions?

A: Pafnuty Lvovich Chebyshev was the preeminent Russian mathematician of the 19th century, born near Moscow and a professor in St. Petersburg. He was remarkable for contributing in a big way to both pure and applied mathematics. The crucial thing for Chebfun was that he laid the foundations of how one can approximate nonperiodic functions by polynomials. What Fourier is to periodic functions, Chebyshev is to nonperiodic functions.

Q: A new version of Chebfun has recently been released. What is different in this release?

A: The big thing is that the code has been completely rewritten in a more modular structure and released on GitHub; we are a fully public open-source project. There’s also a beautiful new web site at www.chebfun.org, with a user’s guide and hundreds of examples showing how to use the code. Besides that there are many speedups and a number of new features, notably the capability of taking advantage of periodicity of functions (i.e., Fourier as opposed to Chebyshev methods). And the two-dimensional side of Chebfun, for representing functions f(x,y), keeps advancing.

Q: Why have you released Chebfun on GitHub?

A: We’d like to encourage new users and new developers. Broadening the Chebfun community seems the best way to make the software develop as strongly as possible and ensure its long-term future. (We have thousands of users already, but we don’t know as much about them as we would like.)

Q: What benefit do you see in the MATLAB File Exchange / GitHub integration?

A: So far as we are aware, Chebfun is the most advanced project written in MATLAB. It’s completely natural that it should be easily available from the File Exchange. We are very pleased at this new development.

Q: How might MATLAB users apply Chebfun?

A: We think every MATLAB user is potentially a Chebfun user, for every MATLAB user works with functions. It’s so much simpler and more natural to call SUM, MAX, ROOTS and so on than to worry about the calling sequences of QUAD, FMINBND, FZERO. In 2D, the same commands find integrals and extrema and roots of functions f(x,y). Chebfun is also the easiest way to solve ODE boundary value problems. We’ve overloaded x=A\b to u=N\f, where N is a differential operator with boundary conditions and f is a right-hand side.

The big advantage in all these computations is that we’re always dealing with functions (chebfuns), they all take the same form, and the user doesn’t have to know anything about the underlying grids. Input to SUM or MAX or ROOTS — just a chebfun. Output from \ for ODEs — just another chebfun. This is very different from the traditional style of computing in MATLAB and other languages with codes like BVP4C/BVP5C where for each code you need to learn a calling sequence that may involve specification of an initial grid and the code’s special structures like “DEVAL” for representing output functions.

Q: How can MATLAB users access Chebfun?

A: Just go to the File Exchange, or google Chebfun. Be sure to check out the People page with its photo-bios of Chebfun developers including Nick Hale, Toby Driscoll, Asgeir Birkisson, Alex Townsend, Anthony Austin, Grady Wright, and Hrothgar.

]]>Denormal floating point numbers and gradual underflow are an underappreciated feature of the IEEE floating point standard. Double precision denormals are so tiny that they are rarely numerically significant, but single precision denormals can be in the range where they affect some otherwise unremarkable computations. Historically, gradual underflow proved to be very controversial during the committee deliberations that developed the standard.

My previous post was mostly about normalized floating point numbers. Recall that normalized numbers can be expressed as

$$ x = \pm (1 + f) \cdot 2^e $$

The *fraction* or *mantissa* $f$ satisfies

$$ 0 \leq f < 1 $$

$f$ must be representable in binary using at most 52 bits for double precision and 23 bits for single precision.

The exponent $e$ is an integer in the interval

$$ -e_{max} < e \leq e_{max} $$

where $e_{max} = 1023$ for double precision and $e_{max} = 127$ for single precision.

The finiteness of $f$ is a limitation on *precision*. The finiteness of $e$ is a limitation on *range*. Any numbers that don't meet these limitations must be approximated by ones that do.

Double precision floating point numbers are stored in a 64-bit word, with 52 bits for $f$, 11 bits for $e$, and 1 bit for the sign of the number. The sign of $e$ is accommodated by storing $e+e_{max}$, which is between $1$ and $2^{11}-2$.

Single precision floating point numbers are stored in a 32-bit word, with 23 bits for $f$, 8 bits for $e$, and 1 bit for the sign of the number. The sign of $e$ is accommodated by storing $e+e_{max}$, which is between $1$ and $2^{8}-2$.

The two extreme values of the exponent field, all zeroes and all ones, are special cases. All zeroes signifies a denormal floating point number, the subject of today's post. All ones, together with a zero fraction, denotes infinity, or `Inf`. And all ones, together with a a nonzero fraction, denotes Not-A-Number, or `NaN`.

My program `floatgui`, available here, shows the distribution of the positive numbers in a model floating point system with variable parameters. The parameter $t$ specifies the number of bits used to store $f$, so $2^t f$ is an integer. The parameters $e_{min}$ and $e_{max}$ specify the range of the exponent.

If you look carefully at the output from `floatgui` shown in the previous post you will see a gap around zero. This is especially apparent in the logarithmic plot, because the logarithmic distribution can never reach zero.

Here is output for slightly different parameters, $t = 3$, $e_{min} = -5$, and $e_{max} = 2$. Howrever, the gap around zero has been filled in with a spot of green. Those are the denormals.

Zoom in on the portion of these toy floating point numbers less than one-half. Now you can see the individual green denormals -- there are eight of them in this case.

Denormal floating point numbers are essentially roundoff errors in normalized numbers near the underflow limit, `realmin`, which is $2^{-e_{max}+1}$. They are equally spaced, with a spacing of `eps*realmin`. Zero is naturally included as the smallest denormal.

Suppose that `x` and `y` are two distinct floating point numbers near to, but larger than, `realmin`. It may well be that their difference, `x - y`, is smaller than `realmin`. For example, in the small `floatgui` system pictured above, `eps = 1/8` and `realmin = 1/32`. The quantities `x = 6/128` and `y = 5/128` are between `1/32` and `1/16`, so they are both above underflow. But `x - y = 1/128` underflows to produce one of the green denormals.

Before the IEEE standard, and even on today's systems that do not comply with the standard, underflows would simply be set to zero. So it would be possible to have the MATLAB expression

x == y

be false, while the expression

x - y == 0

be true.

On machines where underflow flushes to zero and division by zero is fatal, this code fragment can produce a division by zero and crash.

if x ~= y z = 1/(x-y); end

Of course, denormals can also be produced by multiplications and divisions that produce a result between `eps*realmin` and `realmin`. In decimal these ranges are

format compact format short e [eps*realmin realmin] [eps('single')*realmin('single') realmin('single')]

ans = 4.9407e-324 2.2251e-308 ans = 1.4013e-45 1.1755e-38

Denormal floating point numbers are stored without the implicit leading one bit,

$$ x = \pm f \cdot 2^{-emax+1}$$

The *fraction* $f$ satisfies

$$ 0 \leq f < 1 $$

And $f$ is represented in binary using 52 bits for double precision and 23 bits for single recision. Note that zero naturally occurs as a denormal.

When you look at a double precision denormal with `format hex` the situation is fairly clear. The rightmost 13 hex characters are the 52 bits of the fraction. The leading bit is the sign. The other 12 bits of the first three hex characters are all zero because they represent the biased exponent, which is zero because $emax$ and the exponent bias were chosen to complement each other. Here are the two largest and two smallest nonzero double precision denormals.

format hex [(1-eps); (1-2*eps); 2*eps; eps; 0]*realmin format long e ans

ans = 000fffffffffffff 000ffffffffffffe 0000000000000002 0000000000000001 0000000000000000 ans = 2.225073858507201e-308 2.225073858507200e-308 9.881312916824931e-324 4.940656458412465e-324 0

The situation is slightly more complicated with single precision because 23 is not a multiple of four. The fraction and exponent fields of a single precision floating point number -- normal or denormal -- share the bits in the third character of the hex display; the biased exponent gets one bit and the first three bits of the 23 bit fraction get the other three. Here are the two largest and two smallest nonzero single precision denormals.

format hex e = eps('single'); r = realmin('single'); [(1-e); (1-2*e); 2*e; e; 0]*r format short e ans

ans = 007fffff 007ffffe 00000002 00000001 00000000 ans = 1.1755e-38 1.1755e-38 2.8026e-45 1.4013e-45 0

Think of the situation this way. The normalized floating point numbers immediately to the right of `realmin`, between `realmin` and `2*realmin`, are equally spaced with a spacing of `eps*realmin`. If just as many numbers, with just the same spacing, are placed to the left of `realmin`, they fill in the gap between there and zero. These are the denormals. They require a slightly different format to represent and slightly different hardware to process.

The IEEE Floating Point Committee was formed in 1977 in Silicon Valley. The participants included representatives of semiconductor manufacturers who were developing the chips that were to become the basis for the personal computers that are so familiar today. As I said in my previous post, the committee was a remarkable case of cooperation among competitors.

Velvel Kahan was the most prominent figure in the committee meetings. He was not only a professor of math and computer science from UC Berkeley, he was also a consultant to Intel and involved in the design of their math coprocessor, the 8087. Some of Velvel's students, not only from the campus, but also who had graduated and were now working for some of the participating companies, were involved.

A proposed standard, written by Kahan, one of his students at Berkeley, Jerome Coonen, and a visiting professor at Berkeley, Harold Stone, known as the KCS draft, reflected the Intel design and was the basis for much of the committee's work.

The committee met frequently, usually in the evening in conference rooms of companies on the San Francisco peninsula. There were also meetings in Austin, Texas, and on the East Coast. The meetings often lasted until well after midnight.

Membership was based on regular attendance. I was personally involved only when I was visiting Stanford, so I was not an official member. But I do remember telling a colleague from Sweden who was coming to the western United States for the first time that there were three sites that he had to be sure to see: the Grand Canyon, Las Vegas, and the IEEE Floating Point Committee.

The denormals in the KCS draft were something new to most of the committee. Velvel said he had experimented with them at the University of Toronto, but that was all. Standards efforts are intended to regularize existing practice, not introduce new designs. Besides, implementing denormals would reguire additional hardware, and additional transistors were a valuable resource in emerging designs. Some experts claimed that including denormals would slow down all floating point arithmetic.

A mathematician from DEC, Mary Payne, led the opposition to KCS. DEC wanted a looser standard that would embrace the floating point that was already available on their VAX. The VAX format was similar to, but not the same as, the KCS proposal. And it did not include denormals.

Discussion went on for a couple of years. Letters of support for KCS from Don Knuth and Jim Wilkinson did not settle the matter. Finally, DEC engaged G. W. (Pete) Stewart, from the University of Maryland. In what must have been a surprise to DEC, Pete also said he thought that the KCS proposal was a good idea. Eventually the entire committee voted to accept a revised version.

Denormal floating point numbers are still the unwanted step children in today's floating point family. I think it is fair to say that the numerical analysis community has failed to make a strong argument for their importance. It is true that they do make some floating point error analyses more elegant. But with magnitudes around $10^{-308}$ the double precision denormals are hardly ever numerically significant in actual computation. Only the single precision denormals around $10^{-38}$ are potentially important.

Outside of MATLAB itself, we encounter processors that have IEEE format for the floating point numbers, but do not conform to the 754 standard when it comes to processing. These processors usually flush underflow to zero and so we can expect different numerical results for any calculations that might ordinarily produce denormals.

We still see processors today that handle denormals with microcode or software. Execution time of MATLAB programs that encounter denormals can degrade significantly on such processors.

The Wikipedia page on denormals has the macros for setting trap handlers to flush underflows to zero in C or Java programs. I hate to think what might happen to MATLAB Mex files with such macros. Kids, don't try this at home.

William Kahan, Home Page.

Charles Severance, An Interview with the Old Man of Floating-Point, Reminiscences for *IEEE Computer*, February 20, 1998.

Wikipedia page, Denormal number.

Get
the MATLAB code

Published with MATLAB® R2014a

This summer my mother-in-law is renting a house on a lake in New Hampshire. Looking at the calendar, my wife said: "The ten-day forecast makes it look like it's going to be pretty hot up at the lake next week." This led to a more general discussion of the merits of ten-day forecasts.

It's funny how we can make decisions based on long-term predictions of weather even though we rarely go back and verify that the forecast was any good. Somehow the fact that the forecast exists at all gives it value. I'm left pondering this question: how much should we trust a ten-day prediction? As it happens, I have some data that can be useful here. For some time, I have been collecting some relevant data on Trendy: the ten day forecast for Natick, Massachusetts (hometown for MathWorks). So let's run some numbers.

Here's the trend: Ten Day Forecast Highs for Natick, MA.

Once a day this trend collects ten data points: today's high temperature and the predicted high temperature for the next nine days. In MATLAB, we'll be working with a matrix with one row for each day and ten columns.

Let's get the data into MATLAB so we can play around with it. I can retrieve (and so can you) the data from Trendy as a JSON object using the following call:

http://www.mathworks.com/matlabcentral/trendy/trends/1655/trend_data.json

In order to read this into MATLAB, I'm going to use Joe Hicklin's JSON parser.

```
url = 'http://www.mathworks.com/matlabcentral/trendy/trends/1655/trend_data.json';
json = urlread(url);
raw = JSON.parse(json);
```

t = zeros(length(raw),1); d = zeros(length(raw),10); for i = 1:length(raw) t(i) = raw{i}{1}; predictions = raw{i}{2}; for j = 1:10 d(i,j) = str2num(predictions{j}); end end firstTenRows = d(1:10,:)

firstTenRows = 44 51 59 66 46 48 53 50 51 54 50 58 63 49 48 52 48 46 57 54 59 61 49 47 53 49 48 43 41 48 62 49 48 54 49 39 39 44 47 46 49 48 54 50 40 38 39 47 51 54 47 55 50 39 40 48 52 53 53 53 54 50 40 39 48 53 54 52 52 50 49 40 38 50 55 54 56 56 53 49 40 39 50 56 54 52 56 54 47 43 39 50 55 55 55 59 58 40 41 46

Now I have a temperature prediction matrix that's structured like this.

I want to re-order this matrix so that each line shows the prediction trajectory for a single day in time. That means picking off the diagonals highlighted in the diagram above. So let's write some code that does this shift. I'm going to end up with two new matrices, d1 and d2

rowIndex = 1:10; colIndex = 10:-1:1; sz = size(d); len = (size(d,1)-10); d1 = zeros(len,10); d2 = zeros(len,10); t1 = zeros(len,1); for i = 1:len ind = sub2ind(sz,rowIndex+i-1,colIndex); trend = d(ind); d1(i,:) = trend; d2(i,:) = trend-trend(end); t1(i) = t(i+9); end firstTenRows = d1(1:10,:)

firstTenRows = 54 57 43 39 38 40 39 38 39 39 54 41 44 39 48 48 50 50 50 49 48 47 47 52 53 55 56 55 54 57 46 51 53 54 54 54 55 55 57 58 54 53 52 56 52 55 57 57 60 58 53 52 56 56 59 59 61 63 63 64 50 53 54 58 45 42 45 44 44 43 49 47 40 39 37 42 43 43 43 43 43 41 41 42 44 48 47 49 48 47 46 44 48 49 46 52 51 50 49 48

In d1, each row is the evolving temperature prediction for each day. So when we plot the first row of d1, we're getting the predictive arc for November 13th of last year.

i = 1; plot(-9:0,d1(i,:)) title(sprintf('Predicted Temperature for %s',datestr(t1(i),1))) xlabel('Time of Prediction (Offset in Days)') ylabel('Predicted Temperature (Deg. F)')

In d2, we just subtract from each row the last value in each row. Since this last value is the final (and presumably correct) temperature, this difference gives us the predictive error across the ten days. Here's the error for the November 13th prediction.

i = 1; plot(-9:0,d2(i,:)) title(sprintf('Error in Predicted Temperature for %s',datestr(t1(i),1))) xlabel('Time of Prediction (Offset in Days)') ylabel('Prediction Error (Deg. F)')

Notice that it shrinks to zero over time. That's good! Our predictions get more accurate as we approach the actual day in question. But the early predictions were off by as much as 18 degrees. Is that good or bad? You tell me.

Now let's look at all the days.

plot(-9:0,d2','Color',[0.5 0.5 1]) title('Error in Predicted Temperature') xlabel('Time of Prediction (Offset in Days)') ylabel('Prediction Error (Deg. F)')

It's hard to get a sense of the error distribution. So let's finish with a histogram of the absolute value of the error. Out of 240 measurements in this data set, the median error for a ten-day prediction is six degrees.

hist(abs(d2(:,1)),1:25) title('Histogram of Error in the Ten-Day Forecast') xlabel('Error (deg. F)') ylabel('Number of Samples')

That seems pretty good. Most of the time that error is going to be less than seven or so degrees Fahrenheit (or four degrees Celsius). I probably don't need to pack a sweater for the weekend at the lake.

Get
the MATLAB code

Published with MATLAB® R2014a

I teach the MATLAB 101 class to every new hire at MathWorks. Inevitably, someone will ask me how to make a global variable. I then Google up: “How to make MATLAB Doug Cry” and find the above post. I think this is a post I should bring up every now and then because these things should be avoided today just like they should have been when I first wrote it.

I will be back from vacation soon with fresh content.

Thanks,
Doug]]>
*Today we have a guest post by Gareth Thomas. Gareth is Academic Marketing Manager for the Education Business in the MathWorks. He aims to inspire people of all ages about technology and share the magic that MathWorks brings to the table when learning, teaching, researching and solving problems that accelerate the pace of engineering and science. You can find Gareth on LinkedIn and Google+.
*

*by Gareth Thomas*

Today I want to talk about sharing tables in reports (a table is a relatively new data type in MATLAB suitable for column-oriented or tabular data). In the current release of MATLAB (R2014a) it can be a little tricky to represent a table variable in an HTML document. So I’m going to show you how to take the reporting capabilities that MATLAB currently offers, combine them with some HTML 5 and JavaScript, and end up with a beautiful report that functions nicely as a standalone document.

Today there are basically two different ways to create HTML reports from MATLAB.

- MATLAB’s publish command
- The MATLAB Report Generator

The publish command ships with MATLAB, while the Report Generator is a product dedicated to the creation of sophisticated reports from MATLAB. Beyond this, the key difference between the two methods really comes down to when you want to deploy your algorithm. That is, do you want the report to update itself in real time? If you want to create a standalone executable, you can really only use the Report Generator.

I’m going to focus on MATLAB’s publish command. Suppose you want to pass a lot of MATLAB-generated HTML into your published document. It turns out that you can use the disp command to add that HTML code into the generated document. However it is not really easy to pass arguments, as this HTML is hard coded at publish time. To get around this I created a function (makeHTMLTableFilter) that creates HTML code that, when published, is also dynamic and responsive. It achieves this by adding Javascript code. This allows you to do things like filter the table interactively.

See the example below:

This is a very nice way of sharing your results in a form of a table while allowing others to filter the fields that they wish.

One question that I get is, “Okay, but how would I deploy this?” As it happens, this idea of using

`disp(makeHTMLTableFilter(mytable))`

is also supported in the MATLAB Report Generator.

So download makeHTMLTableFilter and give it a try. You’ll soon be making interactive reports of your own.

]]>