Matlab and Ptolemy/Matlab Interface

November 9, 1995, DSP Design Group Meeting

Introduction - Demonstrations - Testing Schedulers - Ptolemy/Matlab Interface - Conclusion

Nov. 16th Talk: Filter Design with Matlab

Nov. 30th Talk: Tool Collaboration with Matlab in Design

Introduction

vectorized interpreted code

>> [1 2 3 4 5]  +  [1 1 1 1 1]

ans =

     2     3     4     5     6

cross between a scripting language and C, plus matrix operations

>> fid = fopen('x.dat','wb')
>> fwrite(fid, x, 'double')
>> fclose(fid)

user sees one data type: complex double-precision floating-point matrix
```
>> dot(x, j*y)

ans =

	0 +15.0000i
```

no pointers, no data structures, automatic memory allocation/deallocation

>> A = [1 2 1; 2 1 3; 4 0 5];
>> y = [3; 2; 1];
>> A \ y

ans =

         0
    1.4000
    0.2000

single and multiple return values for the same function

>> eigvalues = eig( [ 1 2; 2 1 ] )

eigvalues =

   -1.0000
   3.0000

>> [eigvectors, eigvalues] = eig( [ 1 2; 2 1 ] )

eigvectors =

    0.7071    0.7071
   -0.7071    0.7071

eigvalues =

   -1.0000         0
         0    3.0000

Matlab will try to load undefined functions, e.g.
```
>> m=mandelbrot(100,.75,.25,.01,.01);
```
Matlab will go looking for the file "mandelbrot.m" on the Matlab path. You can append to the Matlab path by either setting the MATLABPATH Unix environment variable or by using Matlab's path command.

Toolboxes
- signal: 1-D/2-D FFT, 1-D analog/digital filter design, filtering, 1-D/2-D convolution, LPC, resampling, spectrogram
  - FFT
```
>> fft( [1 1 1 1] )

ans =

     4     0     0     0
```
  - Design a lowpass fourth-order elliptic analog filter with a cutoff of 10 Hz with 1 dB passband ripple and 20 db of stopband attenuation:
```
>> [poles,zeroes,gain] = ellip(4, 1, 20, 10, 's')

poles =

        0 +11.2433i
        0 -11.2433i
        0 +20.3918i
        0 -20.3918i


zeroes =

  -4.0026 + 6.5087i
  -4.0026 - 6.5087i
  -0.5162 +10.0365i
  -0.5162 -10.0365i


gain =

    0.1000
```
- image: DCT, Radon, 2-D FIR filter design, morphology, median filtering, enhancement, statistics
```
>> dct2( [1 0; 0 1] )
ans =
    8.0000         0
         0    4.0000
```
- optimization: linear, quadratic, and sequential quadratic programming
- symbolic: interface to Maple (e.g., computing symbolic derivatives)
Mex Files
- C functions dynamically linked into Matlab (or C++ functions with extern "C" specifier)
Simulink
- run "matlab" and then type "simulink" at the Matlab prompt
- slow GUI block diagram interface to Matlab
- cannot use Matlab routines directly: one must append four arguments to every Matlab call to update state in the block diagram
- DSP block set (C code generation)
- Fixed-point block set (8-, 16-, and 32-bit word lengths)
Future
- Matlab 5: object-oriented
- Compiled scripts
Third-party Uses
- 20 controls, signal processing, and communications textbooks
- Cosimulating MATLAB and COSSAP: converts Matlab script to COSSAP netlist
- MAT2DSP: RASSP Project at U.C. Davis to profile primitive operations for approximating cost of implementation on a DSP architecture
- Ptolemy/Matlab interface for simulation

Demonstrations

For Praveen's demonstrations, you will need to set the Matlab path. In Unix, evaluate the following before running "matlab":

setenv MATLABPATH /users/murthy/matlab:/users/murthy/matlab/graphAlgs
setenv MATLABPATH ${MATLABPATH}:/users/murthy/matlab/numberAlgs
setenv MATLABPATH ${MATLABPATH}:/users/murthy/matlab/maxAlgebra
setenv MATLABPATH ${MATLABPATH}:/users/murthy/matlab/dagMinbufS

filtdemo (from the Signal Processing Toolbox)
Wondrous numbers (Praveen)
This is a very simple example that illustrates how a mathematical recreation can be quickly programmed and visualized.

A wondrous number (from "Godel, Escher, and Bach" by Hofstadter) is the number that causes the following algorithm to never terminate:
```
>> help wondrous

   function x=wondrous(m)
   Computes the sequence 
   While (m > 1)
   	(if m = ODD), x(i++)=3m+1, m=3m+1, else x(i++)=m/2, m=m/2,  end
   end
   and plots x.  This function is a classic example of one where it is unknown
   if it terminates for all m or whether there are m for which is provably
   does not terminate

>> !more wondrous.m

x=[];
while ( m > 1)
	if (m/2 ~= ceil(m/2))
		x = [x,3*m+1];
		m = 3*m+1;
	else
		x = [x,m/2];
		m = m/2;
	end
end
plot(x)
```
Notice that the above example uses the list-as-dynamically-increasing-vector data structure. This is one reason why prototyping is so easy with Matlab. There is no need for any fancy data-structures for handling these routine occurrences of array sizes not being known in advance.
The plot enables one to quickly see (by invoking the algorithm for different m) that a proof of whether there exists a wondrous m that will cause this algorithm to never terminate is likely to be very difficult, if not impossible (since the values that m takes on during the execution of the algorithm traces very different curves for different starting values).

Mandlebrot images (Praveen)

This is another example of a simple mathematical algorithm that can be easily coded and visualized in Matlab. The algorithm uses complex numbers but because the only data-structure in Matlab is a complex-valued matrix, this is no problem.

% more mandelbrot.m

function m=mandelbrot(n,cx,cy,stepx,stepy)
 
 s2=sqrt(2);
 k=0;
 c=cx+i*cy;       % i = sqrt(-1)
 m=zeros(n,n);
 for x=cx:stepx:cx+(n-1)*stepx
	k = k+1;
	j=0;
	for y=cy:stepy:cy+(n-1)*stepy
		j=j+1;
		z=x+y*i;
		numIts=0;
		while (abs(z) < s2 & numIts < 50)
			z = z*z + c;
			numIts = numIts + 1;
		end
		m(k,j)=numIts;
	end
end

Then, in Matlab, we can invoke this function and visualize it.

>> m=mandelbrot(100,.75,.25,.01,.01);
>> surf(m);
>> view(2);

The colormap can be changed easily to see different color assignments:

>> colormap(hot); %"Hot" colors
>> colormap(rand(64,3));  % Set the colormap to random rgb values

Testing Schedulers (Praveen Murthy)

I have used Matlab for a lot of scheduling related prototyping. It is the language I use most often because it is a) interpreted, not compiled, and b) has extensive plotting and visualization tools making it unnecessary to hack anything on your own. While Matlab's use in signal processing and controls is well known (because of the extensive toolboxes that are available), it is perhaps surprising that it should be useful for graph-theoretic algorithms and such but it is really a simple programming language that is useful for prototyping for me. In fact, Friedlander at UC-Davis has been prototyping a compiler using Matlab!

Most schedulers are rather boring since they take graph incidence matrices as input and return...a list containing when nodes should be fired.

However, I had an idea one morning, prototyped it one afternoon, and forgot about pursuing it that night after it became apparent that the problem was too hard. Anyway, this scheduler just does a dynamic multiprocessor schedule for an HSDF graph greedily (assign tasks to processors when they become available, break ties using some reasonable criterion which escapes me now). The point though was to see how the dynamic schedule evolves, and what its steady state period is, and whether it is possible to predict theoretically, what the optimum "blocking factor" is when the number of processors is restricted. The idea then would be to take as the schedule, the periodic pattern that emerges. An attraction of this approach is that IPC cost can also be taken into account, and as long as the procedure is deterministic (i.e, it doesn't start flipping coins anywhere), the periodic steady state behavior holds. It would be interesting to see how this periodic pattern compares to, say, Sih's heuristic which only takes one iteration into account. The Matlab prototype showed that some very interesting patterns do emerge, but they do not seem to correspond in any way to the graphs (or max-algebra matrices) that arise in my blocking factor paper where the question was explored for the case of an unrestricted number of processors.

In Matlab, I invoke:

>> help asap

function schedule=asap(A0,A1,execTimes,m,numTicks)
Generate an ASAP schedule for a graph G given m processors.
G is a unirate, SDF graph.
Inputs: A0 = G with 0 delay arcs
	A1 = G with 1 delay arcs
	execTimes = execution times of nodes
	m = # procs
	numTicks = number of clock ticks to simulate
Outputs:schedule

>> a0=randdag(6);  % Create a random 6-node acyclic graph
>> a1=zeros(6,6);  % For simplicity, this graph does not have any delay edges
>> et=[1,3,2,1,2,3]; % Arbitrarily chosen execTimes
>> s=asap(a0,a1,et,2,70) % Simulate for 70 ticks, on 3 processors

s =
 
 1b  ac e f  b  ac e f  b  dac e f  b  dac e f  b  dac e f  b  dac e f  b  
 2db  dac e f  db  dac e f  b  dac e f  b  dac e f  b  dac e f  b  dac e

The crude gantt chart shows that a cyclostatic pattern emerges after a short transient.

Ptolemy Interface to Matlab

SDF domain
- multi-input/multi-output SDFMatlab_M and SDFMatlabCx_M stars
- Matlab star as source, input/output, and sink block
- Parameters
  - ScriptDirectory
  - DeleteOldFigures
  - MatlabFunction
- The star will automatically build up a Matlab command: a star with two outputs and one input and a MatlabFunction parameter set to "eig" for eigenanalysis would be built as
```
[output#1, output#2] = eig( input#1 )
```
- The MatlabFunction parameter can contain arbitrary Maltab scripts as well, with commands separated by semicolons.
Implemented as a base class, so usable by all of Ptolemy
- DE Domain
- Tcl/Tk
- Parameter computations

Conclusion

Last modified 04/10/97.