March 29, 2012

simple algorithm for 2D-Peakfinding

When it comes to analyzing data, a common task is to identify peaks in an x-y dataset. There are several ways to do that. Maybe the first idea that would come into one's mind is to look for zeros in the derivatives. This works fine for smooth theoretical curves but fails when noise is present (as it is certainly all kinds of experimental data).
To get rid of noise you can apply a low-pass filter on your data, but in some cases this is not desired. A simple approach is to use an algorithm that can does

  1. find peak candidates that are above some user-defined threshold level
  2. find out for each candidate whether it is the local maximum of a user defined environment
  3. return the valid peaks
the algorithm is able to find peaks even in noisy data


I've implemented such piece of code for octave that does the job. The code below created the picture shown. The function expects three input parameters: the environment for maxima, the data list and a threshold value.



function peaklist = SimplePeakFind ( environment, data,  thresh)
%
% ----------------------------------------------------------
%
% function peaklist = SimplePeakFind (  environment, data, thresh)
%
% finds the positions of peaks in a
% given data list. valid peaks
% are searched to be greater than 
% the threshold value.
% peaks are searched to be the maximum in a certain environment
% of values in the list.
%
%  ----------------------------------------------------------

  listlength = length( data );  
  peaklist = zeros(listlength,1); % create blank output
  SearchEnvHalf = max(1,floor( environment/2));
  %
  % we only have to consider date above the threshold 
  %
  dataAboveThreshIndx = find (data >= thresh);
  for CandidateIndx = 1 : length(dataAboveThreshIndx)
    Indx = dataAboveThreshIndx(CandidateIndx);  
    %
    % consider list boundaries
    %
    minindx = Indx - SearchEnvHalf;
    maxindx = Indx + SearchEnvHalf;
    if (minindx < 1)
      minindx = 1;
    end
    if ( maxindx>listlength)
      maxindx = listlength;
    end
    %
    % data( CandidateIndx ) == maximum in environment?
    %
    if (data(Indx) == max( data( minindx : maxindx)))
      peaklist(Indx) = Indx;      
    end
  end 
  % finally shrink list to non-zero-values
  peaklist = peaklist( find( peaklist));
end


In many cases, this works fine for me. Let's look at a simple test case: we numerically create a number of gaussian pulses with random widths, positions and heights. Then we add noise. 


clear all;
more off;
N = 1500;
T = linspace(-20,20,N);  %x-vector
% create some random gaussian peaks
Npeaks = 9
ppos = T(round(rand(1,Npeaks)*N))
pheight = rand(1,15)*1;
pdur    = rand(1,15)*1;
data = zeros(1,N);
for i = 1:Npeaks
  data = data + pheight(i) * exp(- ((T-ppos(i))/pdur(i)).^2);  
endfor
% apply noisenoise = rand(1,N)*0.1;
datanoise = data+noise;
% Find peaks that are maxima in an
% environment of 25 data points and that
% bigger than 0.15 (noise level)
peaklist = SimplePeakFind(25, datanoise, 0.2);
%Plot result
figure(1)
plot(T,datanoise,"linewidth",2, 
     T(peaklist), datanoise(peaklist),
     "o","markersize",5,"markerfacecolor",[1,0,0])
legend("noisy data","found peaks")
xlabel "x"
ylabel "y"
set(gca(),"tickdir","out","linewidth",2,"ticklength",[0.005,0.005])

The result is shown in the picture above: all 9 peaks are identified correctly. Depending on your data, you'll have to play around with the parameters environment and thresh. Having the the peak positions, further data analysis (like curve fitting etc) is much easier.


Links: of course people already thought of this problem:



March 26, 2012

simulating nonlinear oscillators using the julia language

In my post 'simulating nonlinear oscillators using octave' i showed you how to ... -- i guess you can read. Today i discovered, that a package similar to the odepkg also is included in julia. You just have to find it.


Basically it works quite similar, with some little tweaks.
There is a file called 'ode.jl' in the 'extras' folder located in the julia source folder. We just have to load this thing and have some ode solvers. The drawback is that up to now there  the only way (i found) to adapt the accuracy is to change the lines 


    rtol = 1.e-5
    atol = 1.e-8

in the function ode23
and
    tol = 1.0e-7
in the definition of oderkf{T}, which is a quite dirty method.

Another issue is that the only arguments of the ode solvers are up to now

(F::Function, tspan::AbstractVector, y_0::AbstractVector)


which means we have to hide additional parameters in our solution vector (and make its derivatives zero).

A code example, that works for me (asymmetric nonlinear potential)



# add the path to your juliasource/extras folder
load("/Users/user/juliadir/extras/ode.jl")


function  rhs_asym(t, x)
  #
  # position and velocity are stored in x[1] and x[2]
  #
  xp = x[1]
  vx = x[2]  
  #
  # the parameters are stored in x[3] and x[4]
  #
  k=x[3]
  alpha =x[4]
  dxp = vx
  dv = - (k * xp - alpha * xp^2)
  dx = [dxp; dv;0;0]   # the two zeros ensure that k and alpha  
                       # will not change
end




# problems can occur when you forget to give these values 
# explicitly as Floats ...
k=1.
alpha =0.2
tend=42.


initialDisplacement = 0.0
initialVelocity     = 1.
xini = [initialDisplacement; initialVelocity; k; alpha]


#
# to get higher accuracy, change the values of rtol atol or tol 
# in the ode.jl-file (or better in a copy of it).
#
# now let's integrate
(T, X) = ode45( rhs_asym, [0; tend], xini)



# we store the whole data in one matrix and save it, 
# so we can load and plot with e.g. octave
#
M=[T1 Xlin]
csvwrite("oscil.csv",M)

March 22, 2012

Simulating nonlinear oscillators using octave

Hi there,

numerics are usefull especially when nonlinear problems are considered. Here, I want to show how to calculate the trajectory for some simple one-dimensional oscillators.
You can imagine such an oscillator as a little mass attached to a spring. There is a point of rest, where all forces vanish. When the oscillator is displaced, a force appears originating from the spring.

When energy is put in the system, the mass oscillates. Depending on the form of the potential (which basically is the integral of the force with respect to the displacement)  we can distinguish harmonic and anharmonic oscillators.

From textbooks (or wikipedia) we know the following items, which we want to check using numerics:
  • for the harmonic (linear) oscillator, the trajectory is a pure sine wave of some frequency
  • for the anharmonic (nonlinear) oscillator, more frequency components appear:
    • there are only odd components for the symmetric nonlinear potential
    • there are odd and even components for the asymmetric nonlinear potential


1. The linear (harmonic) oscillator
This is a well known textbook example. The potential energy of the harmonic oscillator is simply a parabula:

Epot = 1/2 k x^2
harmonic potential. The potential energy grows
quadratically with the displacement.

The resulting force is the negative derivative of the potential:


F = - dE/ds = - k x

This is Hooke's law which is a simple model for springs.

Force acting in the harmonic potential
It is also the reason why the harmonic oscillator is called the linear oscillator: the force is directly proportional to the displacement (linear dependency).
The equation of motion is can be derived from the force

d^2 x /dt^2 = F/m = -k/m x

where x is the displacement and k the spring constant.  We can easily give the solution for the equation:

x(t) = A1 cos( k^2/m^2   t   + A2) 

where A1 and A2 are two constants that have to be adapted to the initial values of the problem. The numerical solution using ODEs can be done following the steps described in one earlier post. As the linear oscillator is a special case of the more general problem i'll give the code later in this text.

After integration, we get a trajectory looking as the following  (k=1, initial position =0, initalspeed = 2)
trajectory of the harmonic oscillator
This looks like a nice, pure sine wave (as expected). To get more insight, we use the fourier fransform
of the trajectory x(t).  You have to be careful using results from octaves ode solvers, because in most cases the output points are not equally spaced in time. Using the interpolation function interp1 can help to linearize the data.
The fourier spectrum looks like the following:


fourier components for the harmonic oscillator:
 there is only one frequency present
As expected, we only see one single frequency component, which means we indeed have a harmonic oscillation here.

2. The anharmonic / nonlinear oscillator
Considering the example with the spring, Hooke's law is valid for small displacements. When the displacement gets bigger, additional terms can appear. We want to consider two cases: a symmetric and an antisymmetric nonlinear potential. We give the potential energy:


Epot,sym = 1/2 k x^2 + 1/3 alpha |x|^3

Epot,asym = 1/2 k x^2 + 1/3 alpha x^3


here alpha is a nonlinearity constant. The curves of the potential are shown below:

Two aharmonic potentials, compared with the
harmonic one. They differ in their symmetry
Again, the forces occuring are the derivatives of the potential with respect to the displacement:


Fsym = - dEpot,sym/ds = - (k x + alpha |x| x)
Fasym= - dEpot,asym/ds = - (k x + alpha x^2)


We see additional dependencies between the force and the displacement  which are of the order of x^2. This is why these oscillators are called nonlinear oscillators. The difference to the linear case becomes clear when we look at the force plot:

Forces derived from the potentials above. Note that the
asymmetric force is equal to the symmetric force
for positive displacements, but differs for negative
displacements.


For small displacements, the forces are almost identical. They begin to differ  considerably, when the displacement grows. For the symmetric potential, the absolute value of the force is equal for -x and x. For the asymmetric potential, the force for negative displacements is smaller than for positive values.

The equations of movements can be written as

(symmetric)
d^2 x /dt^2 = F/m = - (k/m x + alpha/m |x| x)

(asymmetric)
d^2 x /dt^2 = F/m = - (k/m x + alpha/m x^2)

Again, we can easily solve this with the odepkg (see code sample below)I used the same input energy as for the harmonic example resulting in  the following trajectories:

Trajectories for the nonlinear oscilators (two kinds), compared with the one
 of the harmonic oscillator


The two new trajectories basically also show  sine osciallations. As the force for the symmetric potential is stronger for big displacements, the maximum amplitude is smaller. This also results in a higher frequency for the oscillation. For the asymmetric potential, the oscillation is shifted  slightly towards negative values in the mean.
We can get the fourier components of this oscillations when analyzing the time traces x(t) of the signals.
We have a closer  look on the symmetric harmonic potential:
spectrum for the symmetric nonlinear oscillator.
Only odd multibles of the base frequency are generated
Besides the fact, that the base frequency is higher than for the harmonic oscillator, we see an additional frequency component! Its value is three times the one of the base frequency (third harmonic). Also, energy is converted into the fifth harmonic.

A similar behavior can be seen for the asymmetric potential:




spectrum for the asymmetric nonlinear oscillator.
Here, even and odd multibles of the base frequency are generated
Here we have generated additional frequencies at two, three and four times the base frequency. The conversion efficiency is higher than for the symmetric potential. 

As a result, the simulations show the expected behavior. We can see that nonlinear oscillators are capable of generating multiples of the base frequency.  This is extensively used in nonlinear optics, where light is sent through nonlinear crystals. The crystal structure is responsible for the potential, in most cases they have asymmetric potentials.
When e.g. light from a laser (at a wavelength of 1064 nm) of high intensity passes such a crystal, green light (532 nm) can be generated.

see 


  • https://en.wikipedia.org/wiki/Second-harmonic_generation



  • https://de.wikipedia.org/wiki/Frequenzverdopplung (german)


  • for further details.


    3. Code to solve the linear and nonlinear oscillators using octave
    clear all;
    more off;


    %
    % right hand side of the ODE for the symmetric harmonic potential
    %
    function dx = rhs_sym(t, x, k, alpha,tend)
      t/tend     %just prints the progress 
      xp = x(1);  % position
      vx = x(2); % velocity
      dxp = vx;   % position change
      dv = -(k * xp + alpha * xp * abs(xp) );  %velocity change
      dx = [dxp, dv]; %give the derivative
    endfunction




    %
    % right hand side of the ODE for the asymmetric harmonic potential
    %
    function dx = rhs_asym(t, x, k, alpha,tend)
      t/tend
      xp = x(1);
      vx = x(2);  
      dxp = vx;
      dv = -(k * xp + alpha * xp^2);
      dx = [dxp, dv];
    endfunction




    % program parameters
    k =1;
    alpha =0.2;
    initialDisplacement = 0;
    initialVelocity     = 2;
    xini = [initialDisplacement, initialVelocity];




    tend = 15;  %integrate up to that time
    % set integration options
    options = odeset( 'InitialStep',tend/1e3,
    'MaxStep',tend/1e3,
    'RelTol',1e-6,
    'AbsTol',1e-6,
    'NormControl','on')




    % linear case, using alpha = 0
    [T1, Xlin] = ode45( @rhs_sym, [0, tend], xini, options, k, 0, tend);


    %NL case, symetric potential:
    [T2, XnlS] = ode45( @rhs_sym, [0, tend], xini, options, k, alpha,\
    tend);
         
    %NL case, symetric potential:  
    [T3, XnlA] = ode45( @rhs_asym, [0, tend], xini, options, k, alpha,\    tend);


    % plot the result
    plot (T1, Xlin(:,1),'color',[0.5,0.5,0.5],
          T2, XnlS(:,1),'.','color',[0,0,1],
          T3, XnlA(:,1),'.','color',[1,0,0]);
    xlabel 'time'
    ylabel 'displacement'




    Hint: for the analysis of the frequency components shown above, the integration time has to be set to a higher value and the odeoptions have to be changed to higher accuracy...

    March 20, 2012

    How to solve ODEs using octave.

    The odepkg package for octave is a great tool to solve ordinary differential equations (ODEs). To solve some problem we only have to bring it into a form octave understands. We can do that following four steps:

    1. determine the variables of the problem
    2. figure out which variables change and how they do it
    3. if higher order derivatives do appear, reduce the order of the system
    4. write a function that can gives the first derivatives for all values considered and solve it with octave's ode solvers
    First example: Shooting a paperball using a slingshot (no friction).

    problems like the trajectory of a sphere fired by a slingshot can be solved using ODEs

    This is a very basic example. In fact, we can give the analytical solution (parabula). The paperball is accelerated by the catapult and flies some distance until it hits the groud. 
    1. The variables of the problem are
    • the position of the paperball in space r(t) = ( x(t), y(t) )
    • the velocity    v(t) = ( vx(t), vy(t) )
    • the acceleration a(t) = (ax(t), ay(t) )
    Here we use the fact that we can split each of the variables of motion into two independent parts, showing in the x or and in the y-direction.

    2. Now we have to consider how our variables change:
    • the change in position is proportional to the velocity at a certain moment:   dr/dt = (vx, vy)
    • the velocity is changed by the acceleration: dv/dt = (ax, ay)
    • for our case, the acceleration is constant. We consider the gravitational force: F=m a where a = g = 9.81 m/s^2 is the gravitational acceleration. It only acts on the y-component of the velocity (perpendicular to the ground)
    As we figured out which variables change in what way, we consider how our result should look like. Octave expects the problem to form a vector (1xN components). What we want to know is the position of the paperball and its velocities for all times t.

    The result should be of the form   result(t) = [rx, ry, vx, vy]

    What we now have to do is to give octave a vector containing the derivative of our result vector. We do this by defining a function.  [drx, dry, dvx, dvy] = deriv_result( time,position)

    3. if we did not already do so, we have to formulate the problem as a system of first order ordinary differential equations.

    What does this mean? Instead of writing the second derivative of some value, we define some additional variable that is  the first derivative of the initial value and consider its changes.

    example: 
    • the acceleration is the second derivative of the position  d^2r / dt^2 = a
    • we introduce the velocity. it is the first derivative of the position dr/dt = v
    • the change of the velocity is given by the acceleration dv/dt = a
    • now we have to solve the problem both for the position and the velocity and such the second order ODE has been reduced to two first order ODE
    4. solving the problem in octave

    We define derivative function

    function dr = dr_gravi(t,r)
      g  = 9.81; %m/s^2
      vx = r(3);
      vy = r(4);
      ax = 0;
      ay = -g;  
      dr = [vx,vy,ax,ay];
    endfunction

    In the program we have to set some initial value. For our problem this is the initial position r0=(x0,y0) and the initial speed v0 = (vx0,vy0)

    X0 = 0
    Y0 = 0
    VX0 = 2
    VY0 = 3

    initialVector = [ X0, Y0, VX0, VY0]

    Also,  some integration window is necessary: for our problem, this is the starting time and the ending time for time integration:

    StartT= 0 %s
    StopT = 0.7 %s

    then we let octave do the job:

    [T,Result] = ode45( @dr_gravi,[StartT, StopT], initialVector)

    Here we used the ode45  solver, but there are many more. See the odepkg documentation for details.
    Our result looks like the following

    T =

       0.00000
       0.07000
       0.14000
       0.21000
       0.28000
       0.35000
       0.42000
       0.49000
       0.56000
       0.63000
       0.70000
       0.70000

    Result =

       0.00000   0.00000   2.00000   3.00000
       0.14000   0.18597   2.00000   2.31330
       0.28000   0.32386   2.00000   1.62660
       0.42000   0.41369   2.00000   0.93990
       0.56000   0.45545   2.00000   0.25320
       0.70000   0.44914   2.00000  -0.43350
       0.84000   0.39476   2.00000  -1.12020
       0.98000   0.29231   2.00000  -1.80690
       1.12000   0.14179   2.00000  -2.49360
       1.26000  -0.05679   2.00000  -3.18030
       1.40000  -0.30345   2.00000  -3.86700
       1.40000  -0.30345   2.00000  -3.86700

    T is the vector that gives us the time in seconds for some row. Results contains the position and velocity components for the times stored in T.  Using this data, we can plot the trajectory and other dynamic information.

    The following code plots the trajectory of the ball and the time dependence of the velocity components:

    figure(1)
      plot(Result(:,1),Result(:,2),'o');
      xlabel 'x position / m'
      ylabel 'y position / m'

    figure(2)
      plot(T, Result(:,3),'.',
         T, Result(:,4),'x');
      xlabel 'time / s'
      ylabel 'velocity / (m /s)'
      legend('vx','vy')



    Is there a way to influence the calculation? Of course there is. We can do that using odeoptions

    options = odeset( 'RelTol',1e-4,
    'AbsTol',1e-4,
    'InitialStep',StopT/1e3,
    'MaxStep',StopT/1e3)
    [T,Result] = ode45( @dr_gravi,[StartT, StopT], initialVector , options )

    There are many options to play around with, of which most are self-explanatory . You have to find a tradeoff between accuracy and calculation time. Again, see the odepkg documentation.

    A more complex example: slingshot with friction.
    Additional to the gravity, let's consider the drag force influencing the flight of our paperball. It is a force proportional to the velocity  v of an object.

    F = -1/2 rho Cd A v^2

    Additional parameters appear, which we set here: 
    • rho, the density of air ( rho = 1.2 kg / m^3)
    • Cd - the drag coefficient . We consider a sphere here, which has Cd = 0.45
    • A is the crosssection of the sphere
    Again, we can split the force into two components: Fx and Fy. To calculate the acceleration components, we also need the mass of our object, as F=ma which leads to a = F/m.

    We can give the derivative function as follows with additional parameters (area and mass of the sphere):


    function dr = dr_gravi_friction(t,r,area,mass)
      g = 9.81; %m / s^2
      cdSphere = 0.45;
      rhoAir = 1.20; %kg / m^3
      frictioncoefficient = 1/2 * rhoAir * cdSphere * area / mass;  
      
      vx = r(3);
      vy = r(4);
      ax = - frictioncoefficient * vx^2;          % only friction

      ay = -( g + frictioncoefficient * vy^2 ) ;  % friction and  
                                                  %  gravitation
      dr = [vx,vy,ax,ay];
    endfunction

    Now we can try out how friction affects the trajectory.  We calculate this for a relative heavy sphere, a light sphere and without friction. The radius of the spheres shall be r=2 cm


    mass1 = 0.1 %kg
    mass2 = 0.001 %kg
    area1 = pi * 0.02^2   %m^2
    area2 = area1;
    options = odeset( 'RelTol',1e-4,
      'AbsTol',1e-4,
      'InitialStep',StopT/1e3,
      'MaxStep',StopT/1e3)

    % trajectory of heavy sphere
    [T1,Result1] = ode45( @dr_gravi_friction,[StartT, StopT], \
        initialVector , options, area1, mass1  )

    % trajectory of light sphere
    [T2,Result2] = ode45( @dr_gravi_friction,[StartT, StopT], \
        initialVector , options,area2, mass2 )
    % trajectory without friction
    [T3,Result3] = ode45( @dr_gravi,[StartT, StopT], \
        initialVector , options)

    When we look at the results, we see the following: The relative impact of the friction force on the light sphere is quite strong. It considerably gets slowed down and hits the ground earlier. The trajectory deviates more and more from the parabula shape while propagating. Contrarily, the heavy sphere almost acts as if  friction was absent.  This basically is what we expected.

    trajectories of spheres

    March 7, 2012

    The 'jet' colormap must die! Or: how to improve your map plots and create your own nice colormaps.


    Introduction.

    We all like the colorful map plots of data, produced by the imagesc function in octave (and matlab). Unfortunately, the standard colormap 'jet' is far from beeing optimal. 
    Octave and Matlab give other colormaps for some purposes, but the great thing is that you are able to freely create your own maps!

    some remarks:

    • i'm sorry for my english
    • the functions shown here should work in GNU octave, a great free and open source matlab clone. With some little changes, this should also be applicable to Matlab.
    • I made quick'n'dirty plots here, which means I disregarded some rules of scaling and size of axis, labels and ticks. I mostly concentrated on colormaps here.

    General structure of colormaps.


    A Octave colormap basically is just a  N x 3  matrix, which you can freely manipulate. The columns of the matrix represent the three channels Red, Green and Blue. N is the number of intermediate values for each color.
    As an example, jet(10) gives us 

    octave-3.4.0:1> m = jet(10)
    m =
       0.00000   0.00000   0.50000
       0.00000   0.00000   0.94444
       0.00000   0.38889   1.00000
       0.00000   0.83333   1.00000
       0.27778   1.00000   0.72222
       0.72222   1.00000   0.27778
       1.00000   0.83333   0.00000
       1.00000   0.38889   0.00000
       0.94444   0.00000   0.00000
       0.50000   0.00000   0.00000


    Each color has values between 0 and 1, representing the intensity of that channel.
    I wrote a little function that visualizes the channels of a given map. 






    function showRGBchannels(Fignr,Map); 

      x = linspace(0,1, size(Map)(1));
      figure(Fignr)
      lw = 4;
      plot( x, Map(:,1),'color',[1,0,0],'linewidth',lw,
            x, Map(:,2),'color',[0,1,0],'linewidth',lw,
            x, Map(:,3),'color',[0,0,1],'linewidth',lw,
            x, mean(Map,2),'color',[0.7,0.7,0.7],'o')
      xlabel 'fraction'
      ylabel 'intensity'
    end


     Let's have a look at the colormap 'summer'.


       showRGBchannels(1,summer(500))

    gives us the following plot:




    color channels of the 'summer' map
    For every map,  'fraction' ranges from 0 to 1 and is mapped to the lowest or highest value of your 3d data. For 'summer', we see that a constant blue value is present. Red and Green increase while fraction increases. 
    The gray line shows the mean intensity value. It's an indication of the 'brightness' of the map, if you want so. Here one disadvantage of the summer map becomes visible: it has a very low contrast. The 'brightness' only changes between 0.3 and 0.8, half of the possible value.


    Comparing colormaps

    To compare the pros and cons of some colormaps, we consider some data to plot. The figure below shows intensity distribution for the Airy Disk , which maybe some of you know from their physics lectures.
    A test function: The Airy disk
    Basically, the Airy distribution has the typical shape of some peak structure, but there are also little features in the wings.


    Why the jet colormap is no good choice here.


    This is the standard plot you would get using the imagesc(data) function
    Airy Disk, plotted with the jet colormap
    When we look at the channels of the jet colormap, we see that the colors appear and also disappear again. The grey curve again gives the mean brightness value. The different regions of the picture are mainly distinguished by the colors.
    Color channels of the jet colormap
    Besides the fact that not the full contrast is used: the average brightness of the colormap is not a monotonously growing or falling curve. After a first rise it drops again! This is a bad thing. Let me explain you why:

    • There are about 5-10 percent of people out there, who have some form of colorblindness. They only will see some reduced projection of the color space. It is no good idea to pronounce different regions of your plot only by using a certain color.
    • Depending on the form of presentation, you will loose color contrast. Consider a projector with poor color quality or some black-and-white printer.
    The two figures below show what happens, when the jet colormap is converted to grayscale. This was done using the freely available program gimp ('colors / desaturate')  with two different algorithms.
    Jet map converted to grayscale I

    Jet map converted to grayscale, II

    The main issue here is that in grayscale representation, information about the 'height' of map regions gets lost. The center of the plot (where the highest intensity is found) has the same gray value as the border regions.


    Using other colormaps (that can be printed in grayscale)
    One of these maps is already included in octave, it is the 'hot' map. Below is the Channelplot of this map:

    RGB channels of the standard hot colormap


    As you can see on the gray curve, the brightness of the map increases monotonously. Now let's see how the Airy disk looks like when plotted using the hot map:
    Airy Disk, plotted using the hot map.

    This already is an improvement, the plot clearly and intuitively points out where the most power of the Airy function is. Unfortunately, the fact that there are more bright rings besides the maximum gets lost. Let's look how we can improve that (see below).

    A simple way to create your own, specialized colormaps
    As you saw above, there often are reasons to differ from the standard map. I already mentioned, that colormap creation basically is easy, as it is only a Nx3 matrix. You can think of all kinds of curves for your different color channels.
    Nevertheless, in most cases it is sufficient to use piecewise linear gradients for the channel function. To make things catchier, i wrote a little function that makes the creation of a colormap simple:


    function mymap = colormapRGBmatrices( N, rm, gm, bm)
      x = linspace(0,1, N);
      rv = interp1( rm(:,1), rm(:,2), x);
      gv = interp1( gm(:,1), gm(:,2), x);
      mv = interp1( bm(:,1), bm(:,2), x);
      mymap = [ rv', gv', mv'];
      %exclude invalid values that could appear
      mymap( isnan(mymap) ) = 0;
      mymap( (mymap>1) ) = 1;
      mymap( (mymap<0) ) = 0;
    end


    It expects 4 input parameters: N is the number of intermediate points that your colormap should have. The other three are matrices that contain the transitions for each channel. Such a matrix should have the following form:

    M = [ 0, b1;
          x1,b2;
          ...
          x_n, b_(n+1);
          1, b_(n+1);
        ];

    the first column give the fractions, where a brightness value should be defined. The second column should contain the brightness levels. Make sure to start the first column at 0 and end it with 1!

    A simple, linear grayscale map can be created with


    M = [0,0;1,1;];
    simplegray = colormapRGBmatrices( 256, M, M, M);

    where we use the same matrix three times (which always results in gray). Now the trick is to figure out, which region of your data you want to emphasize. Knowing the fraction in which the magic happens, you can set accents using sharp color or brightness transitions. This mostly can be done by setting two or three extra points per matrix, which is no big effort.

    Simple example: improve the hot colormap.
    Let's go back to the airy disk. As we saw in the figure above, the hot colormap gives some natural feeling about intensities. Now we want to pronounce the second intensity ring a little bit more. This can be done by adding one extra node point in the red channel, which creates a strong gradient in color and brightness. We set the matrices as follows


    MR=[0,0; 
        0.02,0.3; %this is the important extra point
        0.3,1;
        1,1];
    MG=[0,0;
        0.3,0; 
        0.7,1;
        1,1];
    MB=[0,0; 
        0.7,0;
        1,1];
    hot2 = colormapRGBmatrices(500,MR,MG,MB);

    showRGBchannels(2,hot2);

    The resulting channelplot looks like this:
    color channels of 'hot2'
    You can see, that on the first some percent of the scale, the red channel is dramatically increased, and so does the overall brightness. This results in a better visibility of the second intensity ring of the airy function.

    Airy disk, plotted using 'hot2'
    There are more easy modifications. Assume, you make a presentation and are bound to some given corporate design (e.g. blue-white or so). You can switch the order of the matrices to generate colormaps exhibiting the overall same brightness gradient but with different colors. Lets make a blueish 'hot2' map:

    bluehot = colormapRGBmatrices(mcolors,MB,MG,MR);


    some blueish 'hot2' colormap
    That was easy, wasn't it? 
    So I hope I could  make clear that thinking about colormaps indeed is worth the time.  Using user-defined colormaps gives you the ability to create plots that fit to your data and pronounce the remarkable features. Besides, it's nice to play around with such things.

    [edit: code for plotting the Airy disk field / intensity]


      clear all;
      more off;
      x = linspace(-20,20,500);


    %compose a grid
      [xx,yy] = meshgrid(x,x);


    % calculate the radius r for xx and yy pairs
    % as the airy disc has a  circular symmetry
      r = sqrt( xx.^2 + yy.^2 );


    %calculate the electrical field
      AiryField = besselj(1,r)./r;
    % intensity = abs( field ) ^2
      AiryIntensity = abs(AiryField).^2;



    % plot them
      figure(1)
      imagesc(x,x,AiryField)


      figure(2)
      imagesc(x,x,AiryIntensity);


    EDIT 2013/13/04

    I stumbled upon a paper from a research group that discusses the choice of colormaps with regard to people that exhibit red-green perception deficiencies:


    EDIT 2013/23/04: 

    There's another nice page which has a dedicated view on the perception of plots/presentations by colorblind people:
    EDIT 2015
     have a look at a talk on the coming default colormap in matplotlib from scipy2015: https://www.youtube.com/watch?v=xAoljeRJ3lU
    Looks really nice to me.

    EDIT 2017 

    - matplotlib has introduced official new maps since 2.0 http://bids.github.io/colormap/
    - also, mathworks updated their standard colormap to 'parula' since version R2014b, which looks quite nice and clean. I may stop my nagging.