A KnowledgeBased Environment for the Design and Analysis of Multidimensional Multirate Signal Processing Algorithms
A THESIS Presented to The Academic Faculty By Brian Lawrence Evans
In Partial Ful llment of the Requirements for the Degree of Doctor of Philosophy in Electrical Engineering Georgia Institute of Technology Copyright c 1993 by Brian Lawrence Evans
June, 1993
A KnowledgeBased Environment for the Design and Analysis of Multidimensional Multirate Signal Processing Algorithms
Approved:
James H. McClellan, Chairman
Ronald W. Schafer
Russell M. Mersereau Date approved by Chairman
To my parents
Acknowledgments
I would like to thank my parents for their encouragement, support, motivation, and drive to search for truth and pursue knowledge. My parents were great examples of the personal bene ts of higher education. I would also like to thank my sisters, Melanie and Wendy, for their concern and support over the years. I would like to thank my advisor and mentor, Dr. James McClellan. He has motivated me to reach my potential as researcher and allowed me to pursue the educational impact of my work. Dr. McClellan has my utmost respect. I would like to thank the members of the reading committee Dr. Russell Mersereau, Dr. Ronald Schafer, and Dr. James McClellan and Dr. Alfred Andrew from the Mathematics Department for thoroughly reading the thesis. The Mathematica initiatives in the signal processing courses here at Georgia Tech would not have been possible without the Mathematicabased Calculus initiatives spearheaded by Dr. Andrew, Dr. Thomas Morley, and Dr. George Cain. I would also like to thank the fth member of my defense committee, Dr. Sudha Yalamanchilli. I have enjoyed the basketball games with and against him. I would like to thank Dr. Barnwell and Dr. Schafer for establishing and maintaining an excellent signal processing program. I would like to thank Dr. David Schwartz for suggesting the use of Mathematica as a platform upon which to build a symbolic signal processing environment. Later, Mathematica turned out to be useful in helping students learn signal processing theory. I would like to thank Dr. James Wiltsie for sponsoring my research at the Georgia
iii
Tech Research Institute GTRI and Dr. Chrystanos Papanicolopulos for hiring me at GTRI for two quarters after I completed my Master's degree. I would like to thank the State of Georgia for my primary nancial support as a research assistant. I would like to thank Dr. Sayle for supporting me as a teaching assistant to develop Mathematica courseware. The remainder of my funding came from the Joint Services Electronics Program under contract DAAL0390C0004. I would like to thank Dr. Je rey Froyd and Dr. Mark Yoder of the RoseHulman Institute of Technology for their feedback over the years. I have appreciated all the comments and suggestions by users of the signal processing packages and Notebooks. I would like to thank two professors from the University of Queensland, Dr. Keith Matthews of the Mathematics Department and Dr. George Havas of the Computer Science Department, for their insights into Smith form decomposition algorithms. I would like to express my appreciation for my closest friends Sam JanToua" Liu, Fatma" Ayhan Sakarya, and Gregory Lee" Schug. I have also valued the friendship of Terry Hope Sack" Crone, Kevin Happy" Hargaden, Faouzi Kossentini, Nancy Nance" Lillo, Keith Mad Dog" Matlack, and Kambiz Nayebi, as well as those in the Emmanuel Catholic Charismatic Prayer Group and the Georgia Tech Turkish Community. I would like to thank Kay Gilstrap and Stacy Schultz for their help. I would like to thank Hal^k Ayd
nolu, Wilson Chung, Jos
Crespo, and Faouzi u g e Kossentini for helping me proofread my thesis. I would also like to thank Faouzi and his wife Faten for their extraordinary help in setting up for the oral defense. I am so thankful to the Creator God for continuing to sustain me physically, emotionally, and spiritually. I am indebted beyond repayment for the love shown by Jesus Christ and the guidance given by the Holy Spirit. God gave me the strength to get through the nal push to nish the thesis. I am also thankful to God for communicating with us in many varied ways. I have tried to keep in mind that what is seen is transitory, and what is unseen lasts forever" II Corinthians 4:20.
iv
Trademarks
Connection Machine is a trademark of Thinking Machines Inc. ILS is a trademark of Signal Technology Inc. Macintosh is a trademark of Apple, Inc. MACSYMA is a trademark of MIT. Maple is a trademark of Waterloo University, Canada. Mathematica is a trademark of Wolfram Research Inc. Matlab is a trademark of The Math Works Inc. MS Windows is a trademark of MicroSoft. Monarch is a trademark of The Athena Group. NeXT is a trademark of NeXT Computer, Inc. N!Power is a trademark of Signal Technology Inc. PostScript is a trademark of Adobe Systems Inc. Signal Processing WorkSystem is a trademark of Comdisco Systems Inc. TEX is a trademark of the American Mathematical Society. Unix is a trademark of AT&T. X Windows is a trademark of MIT.
v
Copyrights
Khoros is copyrighted by the University of New Mexico. Blosim, Gabriel, and Ptolemy are copyrighted by the University of California at Berkeley. The multidimensional signal processing packages and Notebooks SPP&N for Mathematica are copyrighted by the Georgia Tech Research Corporation.
vi
Contents
Acknowledgments Trademarks Copyrights Contents List Of Figures List Of Tables Summary 1 Introduction
1.1 1.2 1.3 1.4 1.5 Rapid Prototyping : : : : : : : : : : : : : : : : : : : : : : : : Algorithm Development : : : : : : : : : : : : : : : : : : : : : The Multidimensional Signal Processing Packages MDSPPs : Layout of the Thesis : : : : : : : : : : : : : : : : : : : : : : : Scope of the Thesis : : : : : : : : : : : : : : : : : : : : : : : :
iv v vi x xii xiii xiv
: : : : : : : : : : : : : : : : : : : :
1
3 3 5 7 8
2 Background
2.1 Representing Signals as Objects : : : : : : : : : : : : : : : : : : : : : 2.2 Algorithm Design Environments : : : : : : : : : : : : : : : : : : : : : 2.3 Mathematica Programming for KBSP : : : : : : : : : : : : : : : : : :
10
13 14 17
vii
CONTENTS
2.3.1 Mathematica : : : : : : : : : : : : : : : : : : : : : : : : : : : : 2.3.2 Mathematica As a Platform for Algorithm Design : : : : : : : 2.4 Summary : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 18 19 23
3 Multidimensional Multirate Signal Processing
3.1 Lattice and Matrix Theory : : : : : : : : : : : : : : : : : : : 3.1.1 De nitions : : : : : : : : : : : : : : : : : : : : : : : : 3.1.2 Uniform Resampling of Lattices : : : : : : : : : : : : 3.1.3 Regular Unimodular Matrices : : : : : : : : : : : : : 3.1.4 Decomposing A Resampling Matrix Into Smith Form 3.1.5 Finding Common Resampling Matrix Factors : : : : 3.2 Resampling Operations in Cascade : : : : : : : : : : : : : : 3.3 The Role of Smith Form Matrices : : : : : : : : : : : : : : : 3.3.1 Computation of Coset Vectors : : : : : : : : : : : : : 3.3.2 Computing the Greatest Common Sublattice : : : : : 3.3.3 Simpli cation of Resamplers in Cascade : : : : : : : 3.4 Commutativity of Resamplers in Cascade : : : : : : : : : : : 3.5 A Comprehensive Set of Rules : : : : : : : : : : : : : : : : : 3.6 Summary : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : :
: : : : : : : : : : : : : :
: : : : : : : : : : : : : :
: : : : : : : : : : : : : :
: : : : : : : : : : : : : :
25
26 27 29 31 31 35 37 39 39 41 42 43 47 59
4 The Signal Processing Packages: Implementing Linear Systems Theory in Mathematica 60
4.1 4.2 4.3 4.4 4.5 De ning Signals : : : : : : : : : : : : De ning Systems : : : : : : : : : : : Plotting Signals and Systems : : : : Continuous and Discrete Convolution Summary : : : : : : : : : : : : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
60 62 64 65 69 71
5 Analysis of Multidimensional Signals and Systems viii
5.1 Extending Signal Properties in ESPLICE and ADE : : : : : : : : : :
70
CONTENTS
5.2 Linear Transforms : : : : : : : : : : : : : : : : : : : : : : : : : : : 5.3 Analysis of OneDimensional Systems Using Transforms : : : : : : 5.3.1 Solving Di erence and Di erential Equations : : : : : : : : 5.3.2 Generalized Signal Analysis : : : : : : : : : : : : : : : : : 5.4 Analysis of Multirate Systems Using Transforms : : : : : : : : : : 5.5 Multidimensional Stability Analysis : : : : : : : : : : : : : : : : : 5.5.1 Symbolic Analysis of Stability : : : : : : : : : : : : : : : : 5.5.2 Graphical Analysis of Stability : : : : : : : : : : : : : : : : 5.6 Analysis of Multidimensional Multirate Systems Using Transforms 5.6.1 Automated TwoDimensional Signal Analysis : : : : : : : 5.6.2 Visualization of Downsampling in Two Dimensions : : : : 5.6.3 Automatic Derivation of Transform Properties : : : : : : : 5.7 Summary : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : : : : : : :
: : : : : : : : : : : : : : : : : : :
72 74 74 77 79 83 83 85 87 87 89 93 94 98 100 100 101 102 103
6 Rearranging Multidimensional Systems
6.1 Extending System Properties in ESPLICE and ADE : : : : : : 6.2 Algorithm Rearrangement : : : : : : : : : : : : : : : : : : : : : 6.2.1 General Procedure : : : : : : : : : : : : : : : : : : : : : 6.2.2 Rearranging A Multidimensional Rational Rate Changer 6.3 Finding Optimal Algorithms : : : : : : : : : : : : : : : : : : : : 6.4 Summary : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
: : : : : :
97
7 Generating Equivalent Code for Algorithms
7.1 Generating TEX Code : : : : : : : : : : : : : : : : : : : : : : : : : : 104 7.2 Generating Complete Ptolemy Simulations : : : : : : : : : : : : : : : 106 7.2.1 Program Synthesis : : : : : : : : : : : : : : : : : : : : : : : : 106 7.2.2 Converting Algebraic Formulas to Working Ptolemy Simulations107 7.2.3 Code Generation After Algorithm Rearrangement : : : : : : : 111 7.3 Summary : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 114
104
ix
8 Interactive Design of TwoDimensional Decimation Systems
8.1 Theory Underlying Decimator Design : : : : : : : : : : : : : : : : : : 116 8.2 Design Examples : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 118 8.3 Summary : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 121 9.1 9.2 9.3 9.4 9.5 Student Handout : : : : : : : : : : : : : Interactive Tutorial Notebooks : : : : : : Notebooks Serving as OnLine Reference Notebooks for SelfEvaluation : : : : : : Summary : : : : : : : : : : : : : : : : :
115
9 Impact on the Engineering Curriculum
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
: : : : :
125
128 128 131 133 134
10 Conclusion
10.1 Contributions : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 140 10.2 Future Research : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 143 A.1 Code Generation for the TwoChannel NonUniform Filter Bank : : : 145 A.2 Code Generation for a More E cient Form of the Filter Bank : : : : : 149
138
Appendix A Ptolemy Simulation Code
145
Appendix B Glossary Bibliography Vita
153 157 167
x
List of Figures
2.1 Categories of Existing Algorithm Support Tools : : : : : : : : : : : : 2.2 Environments In uenced by the Signal Representation Language : : : 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 Iterative Algorithm in the MDSPPs to Compute the Smith Form : : Examples of Di erent Smith Forms : : : : : : : : : : : : : : : : : : : Equivalent Structures for " L = KL0 Followed by M = KM 0 : : : : Five Equivalent Forms of an Up downsampler Cascade : : : : : : : : Four Equivalent Forms of a Down upsampler Cascade : : : : : : : : Identities for Downsamplers : : : : : : : : : : : : : : : : : : : : : : : Identities for Upsamplers : : : : : : : : : : : : : : : : : : : : : : : : : Interactions between Up downsamplers and LTI Filters : : : : : : : : Identities for Cascades of Upsamplers and Downsamplers : : : : : : : Commutativity of Cascades of Upsamplers and Downsamplers : : : : Fundamental Identities Based on the Smith Form Decomposition : : Removing Redundancy in Cascades of Upsamplers and Downsamplers Switching Operations in Noncommutable Cascades : : : : : : : : : : Interaction between Up downsampler Cascades and Shifters : : : : : Polyphase Implementations of Rational Decimation Systems : : : : : 10 12 33 34 38 44 45 48 49 50 51 52 53 54 56 57 58 66 75 76
4.1 Visualizing TwoDimensional Sequences as Density Plots : : : : : : : 5.1 Structure of the OneDimensional Transform Rule Bases : : : : : : : 5.2 Interaction with the Di erential Equation Solver : : : : : : : : : : : :
xi
5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 7.1 7.2 7.3 7.4 8.1 8.2 8.3 8.4
Solving the Fibonacci Di erence Equation : : : : : : : : : : : : : : OneDimensional Analog Signal Analysis : : : : : : : : : : : : : : : OneDimensional DiscreteTime Analysis of a Multirate Signal : : Deriving InputOutput Relationship for a NonUniform Filter Bank Stability Analysis of a NonSeparable TwoDimensional Signal : : : TwoDimensional PoleZero Diagrams for an Unstable Signal : : : Signal Analysis of Aliasing in Two Dimensions : : : : : : : : : : : : Quincunx Downsampling Without Aliasing : : : : : : : : : : : : : : Quincunx Downsampling With Aliasing : : : : : : : : : : : : : : :
: : : : : : : : :
78 80 81 82 84 86 88 90 92 109 110 112 113 115 120 122 123
Algorithm to Convert Algebraic Expressions to Ptolemy Simulations Ptolemy Simulation of Filter Bank Run From Within Mathematica : Deriving the Input Bandpass Signal for the Filter Bank Simulation : Block Diagram Form of the TwoChannel Filter Bank : : : : : : : : : Flow Graph of a TwoDimensional Decimator : : : : : : : : : : : : : Automatic Design of a Quincunx Decimator : : : : : : : : : : : : : : Automatic Design of a Decimator for an ArbitrarilyShaped Passband Automatic Design of a Decimator for Circularly Bandlimited Signals
9.1 Animation of Filter Response for Corrupted Poles : : : : : : : : : : 132 10.1 Descendants of Kopec's Signal Representation Language : : : : : : : 139
xii
List of Tables
1.1 The Prototyping Process and Examples of Supporting Tools : : : : : 2.1 2.2 2.3 2.4 Signal Properties in SPLICE, ESPLICE and ADE System Properties in ESPLICE and ADE : : : : : Overview of Algorithm Design Environments : : : : Mathematica Operators That Mimic C Syntax : : : 4 14 15 16 19 61 63 68 73 95 99
: : : :
: : : :
: : : :
: : : :
: : : :
: : : :
: : : :
: : : :
: : : :
: : : :
4.1 Signals Introduced by the Signal Processing Packages : : : : : : : : : 4.2 New Operators and Their Parameters in the MDSPPs : : : : : : : 4.3 Modi cation to the Square Matrix Rule for Convolution : : : : : : : 5.1 Signal Properties Supported by the Signal Processing Packages : : : : 5.2 HighLevel Abilities of the Signal Processing Packages : : : : : : : : 6.1 System Properties in the MDSPPs : : : : : : : : : : : : : : : : : : :
9.1 Coverage of Linear Systems Topics by Tutorial Notebooks : : : : : : 130 9.2 Past and Present Uses of Our Signal Processing Extensions : : : : : : 136 10.1 Contributions of the Thesis : : : : : : : : : : : : : : : : : : : : : : : 141 10.2 Future Research : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 144
xiii
Summary
This thesis discusses the design and analysis by computer of algorithms composed of linear periodically timevarying LPTV multidimensional systems. Analysis of linear systems is based on linear transforms e.g. z and Laplace transforms. Algorithm design rewrites component systems to reduce the implementation cost. To support algorithm design for multidimensional systems, the thesis derives the rules for rewriting the interconnections of discretetime LPTV multidimensional systems, a.k.a. multidimensional multirate systems, as well as develops metrics to measure implementation costs. We encode the rewrite rules, number theoretic algorithms underlying the rules, and cost metrics in a set of multidimensional signal processing packages MDSPPs for the computer algebra program Mathematica. For algorithm analysis, the MDSPPs implement the multidimensional multisided forms of the commonly used linear transforms, and the transforms can justify their answers with natural language. Using the transforms, the MDSPPs can deduce ranges on free parameters to guarantee stability and generate inputoutput relationships. Engineers can use the MDSPPs to visualize signals in transform domains. The MDSPPs represent signals as functions and systems as operators. In such an algebraic framework, the many interconnections in complex systems cannot be captured. The MDSPPs can, however, convert an algorithm for layout in the Ptolemy block diagram environment, which ts the MDSPPs into a rapid prototyping process. Please direct any electronic correspondence concerning this research to Brian Evans at evans@eedsp.gatech.edu. A freely distributable version of the MDSPPs and Notebooks is available via anonymous FTP to gauss.eedsp.gatech.edu IP 130.207.226.24. An estimated 10,000 people use these extensions.
xiv
CHAPTER 1 Introduction
Signal processing algorithms are often expressed as either ow graphs or algebraic formulas. As ow graphs, algorithms are decomposed into blocks with each block consisting of a set of input signals, an operator, and a set of output signals. In an algebraic framework, algorithms are decomposed into a sequence of formulas containing nested calls to operators systems which take functions signals as arguments. Operators map one or more input signals to one or more output signals. A primary drawback of an algebraic representation of algorithms is that it cannot capture the complex interconnections characteristic of large systems that are easily represented by block diagrams. This thesis describes a software environment that takes advantage of both algebraic and block diagram representations of algorithms by automating the design and analysis of algebraic forms of signal processing algorithms and
converting the algebraic forms of algorithms to block diagram representations for simulation and prototyping.
In the new environment, analysis of signal processing formulas is primarily based on linear transforms of multidimensional signals the z, discretetime Fourier, and discrete Fourier transforms of discretetime signals, and the Laplace and Fourier transforms of continuoustime signals. Design capabilities include the interactive graphical design of twodimensional decimation systems and the symbolic rearrangement of linear, multidimensional, periodically timevarying, discretetime systems, i.e. multidimensional multirate systems. The conversion of algorithms to a block diagram
1
Introduction
representation enables designers to insert developed algorithms into a layout of a complex system. Our new environment encodes a comprehensive collection of rules and routines underlying signal processing in a set of multidimensional signal processing packages MDSPPs 1, 2, 3 for the MathematicaTM 4 computer algebra program. For example, the MDSPPs can apply commutative and associative rules to switch the order of linear shiftinvariant operators 5 . The MDSPPs implement many key routines such as the family of Smith form integer matrix decompositions which unlocks e cient implementations of the generalized multidimensional discrete Fourier transform when the decompositions are applied to the underlying sampling matrix 6 . In combining the rules and routines underlying signal processing into a single framework, we have been in uenced by the computing environments for
algorithm design especially SPLICE 7, 8, 9 and ADE 10, 11, 12 , symbolic mathematics especially Mathematica, and rapid prototyping especially Ptolemy 13, 14, 15, 16 .
SPLICE and ADE are interpretive environments for the design of onedimensional multirate signal processing algorithms represented as formulas. These two environments can analyze properties of signals and nd better implementations of algorithms. SPLICE and ADE as well as Mathematica and other computer algebra programs represent algorithms as formulas so they cannot capture the interconnections inherent in large complex systems. Using hierarchical block diagram representations, Ptolemy can represent and simulate complex systems as well as generate the equivalent C, DSP assembly language, and VLSI Hardware Descriptive Language VHDL code for complex systems. Essentially, our MDSPPs not only extend the ideas present in SPLICE and ADE to multidimensional multirate systems but also t into a rapid prototyping process for signal processing algorithms.
2
Rapid Prototyping
1.1 Rapid Prototyping
In designing signal processing systems, signal processing theory is translated into hardware and software. That is, the mathematical principles underlying the signal processing theory are combined with knowledge of computer architectures and programming to produce a working system. The design process typically begins by translating speci cations into signal processing operations. Once the component operations are identi ed, they are individually designed by modifying old algorithms or by developing entirely new ones. The complete system is usually simulated before it is assembled into a prototype. Table 1.1 outlines the prototyping process. Design tools, as shown in Table 1.1, often address particular aspects of prototyping, but they rarely have the ability to transfer design information directly to other computeraided design CAD tools. Hence, the prototyping process becomes fragmented because algorithms cannot be exchanged between environments. Furthermore, side information about algorithms representing lessons learned during one prototyping stage cannot be electronically transmitted to other stages. Today, many isolated" tools are used in prototyping signal processing systems, thereby driving up prototyping costs. For example, in a recent moving target Radar system, the signal processor accounted for over 50 of the total material cost of production" page 10 of 17 . Simply put, the complexity of prototyping signal processing processors has outpaced the sophistication of and interaction between CAD tools 17 .
1.2 Algorithm Development
In this thesis, we discuss a software environment that is useful in the algorithm development phase of prototyping. By algorithm development, we mean more than numerical simulation of algorithms, i.e. tweaking numeric parameters until the algorithm performs well enough. We also mean the symbolic reasoning about algorithms that occurs before and after simulation. Symbolic reasoning involves analyzing and
3
Algorithm Development
Development Typical Tool Activity Functions Availability System Requirements Traceability System Requirements Flowdown Sparse Engineering Algorithm Development Module Physical Design & Layout Design Upgrade Simulate Fabricate Integrate Final Test VHDL Synthesis Interactive Software Development Environment VHDL Simulation Module Simulation Pick and Place Systems, with CAD design download capability Tester Interfacing High delity simulation of outside world" for functional performance test Generation of boundary scan testing at up to the system level Extensive
Typical Sources GEC Marconi iLogix MathWorks Mentor Synopsis Centerline Vantage LSI Logic N. American Phillips Test Software Services Inc. TSSI XpertTest Hewlett Packard
Extensive Extensive Less Extensive Extensive Evolving
Table 1.1: The Prototyping Process and Examples of Supporting Tools adapted from Table 1 in 17 ; see Appendix A of 17 for more information
4
The Multidimensional Signal Processing Packages MDSPPs
nding optimal implementations for algorithms. Analysis characterizes signals and systems in terms of their properties which tend to be symbolic relationships rather than numeric values. Optimal implementation seeks to rearrange the form of the algorithm until some cost function is minimized. Chapter 2 reviews three environments for analyzing and nding optimal implementation of algorithms the Signal Processing Language and Interactive Computing Environment SPLICE, the Algorithm Design Environment ADE, and Meta Morphology MetaMorph 18, 19, 20, 21 . As previously mentioned, SPLICE and ADE manipulate onedimensional signals and multirate systems and attach properties to signals and systems. SPLICE and ADE can deduce signal properties such as bandwidth and extent. On the other hand, system properties are used to rearrange systems, e.g. the order of two systems in cascade can be switched if both systems are linear and shiftinvariant. By applying rules that rearrange cascaded and parallel combinations of systems, Extended SPLICE ESPLICE and ADE nd optimal implementations of onedimensional multirate systems by minimizing a simple cost function. MetaMorph manipulates morphological signals sets and systems nonlinear set operators in a similar way. In order to t into a rapid prototyping process, however, an environment should be able to generate equivalent code for algorithms that would be suitable for incorporation by tools at the next level of prototyping. Unfortunately, ADE, SPLICE, and MetaMorph cannot generate code for algorithms.
1.3 The Multidimensional Signal Processing Packages MDSPPs
The goal in creating our new algorithm development environment was to complement the abilities of the SPLICE, ADE, and MetaMorph environments. As previously mentioned, our new environment is a set of multidimensional signal processing packages MDSPPs for Mathematica that analyze and nd optimal implementations of
5
The Multidimensional Signal Processing Packages MDSPPs
linear, multidimensional, multirate signal processing algorithms. Unlike SPLICE and ADE, which only run on Symbolics machines, the MDSPPs are available for many computing platforms because they run wherever Mathematica does. Unlike SPLICE, ADE, and MetaMorph, the MDSPPs t into a rapid prototyping framework because it generates block diagram descriptions for algorithms which can be inserted into the layout of a complex system in Ptolemy. We chose Ptolemy because it can represent and simulate complex systems, as well as generate the equivalent C code, DSP assembly programs, and VHDL descriptions for complex systems. Although SPLICE, ADE, and MetaMorph have been written in Lisp, we have chosen to start at a higher level by using the symbolic mathematics program Mathematica. Mathematica o ers many bene ts. Mathematica is pro cient at partial fraction decomposition, power series expansions, symbolic integration, and in nite summations, which are fundamental operations in convolution and linear transforms. Mathematica also provides a mechanism for expressing and applying rules. Because of the various plotting and graphics formats in Mathematica, users of the MDSPPs can visualize in di erent domains signals generated by onedimensional and twodimensional multirate systems. Mathematica requires any routine accessible by the user to have usage information attached to it, so all new functions come with help information. Mathematica also supports the playing of sound and animation. On several platforms, Mathematica o ers a multimedia Notebook user interface 22 that organizes formatted text, sound, graphics, animation, and Mathematica code in a tree structure. Using the Notebook interface, we have written a user's guide, several tutorials, and online help in the form of interactive electronic documents 3, 23, 24, 25, 26 .
6
Layout of the Thesis
1.4 Layout of the Thesis
As previously mentioned, Chapter 2 discusses the SPLICE, ADE, and MetaMorph algorithm design environments and evaluates Mathematica as a platform upon which to build an algorithm design environment. Chapter 3 uni es the theory underlying the operation of and interaction between multidimensional multirate structures. Based on this new theory, we develop the rules that generalize the interaction between up downsamplers and linear shiftinvariant operators to multiple dimensions in such a way as to be easily implementable by computer. Chapter 4 introduces the representation of signals and systems in the MDSPPs. The subsequent chapters discuss the algorithm development capabilities of the MDSPPs. The MDSPPs 1. analyze linear multidimensional multirate structures Chapter 5, 2. nd optimal implementations for algorithms involving these structures Chapter 6, 3. generate equivalent TEXTM and Ptolemy code for these structures Chapter 7, and 4. assist in the interactive design of twodimensional rational decimation systems Chapter 8. The ability to generate complete Ptolemy simulations for algorithms links the MDSPPs to the rapid prototyping process. The higherlevel functions of the MDSPPs feature the ability to justify their answers in the form of textual and or graphical dialogue. For example, the linear transform routines can show the intermediate hand calculations involved in taking z, Laplace, and other transforms, whereas the convolution routines can illustrate by animation the ipandslide approach to convolving two piecewise functions. These
7
Scope of the Thesis
dialoguing capabilities have helped students learn the theory underlying signal processing 3, 23, 24, 25, 26 . The Notebooks accompanying the MDSPPs guide a student's exploration of signal processing theory. Chapter 9 assesses the impact of the our new environment on engineering education. Chapter 10 highlights the contributions of the thesis and summarizes areas of future research. Appendix A gives two examples of Ptolemy code generation for a twochannel onedimensional nonuniform lter bank. Appendix B is a glossary of terms.
1.5 Scope of the Thesis
This thesis treats several diverse topics that may not appeal to all readers. Readers with a signal processing background would probably be the most interested in Chapters 3, 5, and 8 as they characterize, analyze, and design multidimensional multirate systems, respectively. They would also be interested in Chapter 4 as it introduces signal processing with formulas by computer, which is a review of the algebraic representation of signals and systems as functions and operators, respectively. Mathematicians would probably be interested in the same chapters namely Chapters 3, 4, and 8 but for di erent reasons. Chapter 3 applies lattice theory and integer matrix algebra to describe the regular insertion and deletion of points on a uniform grid known as upsampling and downsampling, respectively. Chapter 4 derives an algorithm that performs convolution with functions described in a piecewise manner each interval consists of a formula that de nes the function on that interval and a pair of endpoints that can be symbolic e.g. in nite. Chapter 8 presents three design examples that require both exact precision and oating point arithmetic in the calculation of free parameters. Although not discussed in this thesis, a mathematician may be interested in the use of a computer algebra environment to derive, encode, and test new number theoretic algorithms for integer matrix decompositions which
8
Scope of the Thesis
we present in 27 . On the other hand, computer scientists would most likely be interested in Chapters 6 and 7. Chapter 6 discusses algorithm rearrangement using forwardchained contextfree rules. Chapter 7 explores automatic code generation program synthesis for formulas using bottomup parsing of the algorithm represented as a tree of operations and a hash table to prevent the generation of redundant operations. Educators would most likely be interested in Chapter 9 because it explores the use of computer algebra systems to teach linear systems transform theory. We have written Mathematica routines so that they can show students how to perform the calculations by hand. Using Mathematica's Notebook interface, we have developed interactive tutorials to enable students to explore new topics. Other Notebooks help students test themselves. In essence, this thesis weaves together concepts from many di erent elds. These concepts were necessary to make an algorithm development environment that is simultaneously useful to designers implementing signal processing algorithms, students learning linear systems theory, and researchers exploring new signal processing theory. Because of the multidisciplinary nature of this thesis, we provide a glossary of terms in Appendix B.
9
CHAPTER 2 Background
Many environments have been developed for the automatic simulation, interpretation, and design of algorithms, as shown in Figure 2.1. In testing algorithms, engineers simulate them by passing a variety of numerical signals into the algorithm until the output signals exhibit the desired characteristics. For some domains, signal interpretation can be performed on numeric signals to produce symbols as the output e.g. words in speech recognition. Sometimes, engineers need assistance deriving new algorithms or rearranging old algorithms to achieve a better performance maybe fewer multiplications and additions per output sample. In this case, the engineer begins with a symbolic possibly mathematical description of an algorithm and then tries to manipulate it into another symbolic description having a more desirable form. Many of the computer tools available for simulation, interpretation, and design of signal processing algorithms originated in the work of Kopec 28, 29, 30 . Kopec developed a general simulation environment for causal numeric signals, but he set
INPUT
Numeric Signal
H H H
Simulation

OUTPUT
Numeric Signal
Hj
KnowledgeBased Signal Processing
H H H H H
Symbols
H

Algorithm Design
Symbols
Figure 2.1: Categories of Existing Algorithm Support Tools
10
Background
the stage for future environments by introducing an objectoriented representation of signals. Signals, as objects, contain slots data and methods procedures. Slots include the signal's history and previously computed values. Methods include those describing how to compute signal values. Kopec's interpretive Signal Representation Language SRL in uenced the development of environments for the simulation, knowledgebased interpretation, and design of signal processing algorithms, as shown in Figure 2.2. All three types of environments, especially the latter two, are discussed in 31 . The rst half of 31 discusses algorithm design environments in detail, whereas this chapter only provides an overview of them. The rst algorithm design environment to follow Kopec's work is the Signal Processing Language and Interpretive Computing Environment SPLICE by Myers 7, 8, 9 . Myers extended Kopec's representation in the direction of abstract" signals. Abstract signals are de ned in terms of their properties rather than in terms of their sample values. Myers also introduced abstract systems. Combining abstract signals and systems led to the ability to simplify and rearrange onedimensional multirate signal processing algorithms. By enumerating equivalent implementations, Extended SPLICE ESPLICE can choose implementations that are optimal" according to a simple cost function. Covell 10, 11, 12 improved upon ESPLICE to create the Algorithm Design Environment ADE. ADE applies heuristics to reduce the number of alternative implementations generated, and it uses a more complicated cost function. Richardson 18, 19, 20, 21 applied the same concepts to morphological signals sets and systems set operators in his MetaMorph environment. Section 2.1 describes Kopec's work. Section 2.2 discusses the SPLICE, ADE, and MetaMorph algorithm design environments. Although these environments have been written in Lisp, we have chosen to write ours in the Mathematica symbolic mathematics program. Section 2.3 evaluates Mathematica as a platform for algorithm design.
11
Background
Signal Representation Language 30
J
Simulation
KBSP
?
Algorithm Design
J J J J
Integrated Signal Proc. System 29 Interactive Laboratory System 32
? ?
The KBSP Package 7 Pitch Detector's Assistant 35
? ?
Signal Proc. Language and Interactive Comp. Environment 8
J^
Blosim 33 Gabriel 34 Ptolemy 13
An Algorithm Design Environment 10
?
Interactive Signal Proc. Under Duress 36 Integrated Proc. and Understanding of Signals 37
?
Meta System for Morphology 20
?
Figure 2.2: Environments In uenced by the Signal Representation Language
12
Representing Signals as Objects
2.1 Representing Signals as Objects
In his thesis, Kopec identi es di erent models for computing and obtaining signal values 28 . As functions, signals can be computed pointbypoint and for e ciency, computed values are stored in a cache for future use. Finitelength signals can be captured in an array. Data streams 38 e.g. modems and statespace models 39 e.g. Kalman lters form two important classes of in niteduration signals. These four classes provide a way of abstracting signal computations. That is, a programmer could ask for values of a signal without having to know how they were computed. Hiding the method of computation is a principle of objectoriented programming. As an object, a signal has data slots and procedures methods attached to it. Under this paradigm, the description of a signal becomes more than a domain and range of values it can now contain information about its parameters, history, computational model, and so forth. This extended representation of a signal is implemented in Kopec's Signal Representation Language SRL 30 . SRL hierarchically classi es signals and then instantiates signals as needed. Once signals are created, their properties remain xed. This immutability of signals is the result of the closure model for signals 31 and is necessary during the generation of alternate implementations. SRL abstracts the representation of signals so as to aid in their computation but restricts signals to be onedimensional niteextent causal sequences beginning at time index 0. Systems are implemented as Lisp routines that map one numeric array into another. SRL can simulate algorithms that use these representations of signals and systems, but it cannot symbolically reason about signal properties. Because interaction with SRL was limited to a command line interface using a Lisp syntax, Kopec developed the Integrated Signal Processing ISP System 29 which provided a graphical user interface to SRL. The ideas behind SRL and ISP were recently implemented in the program N!PowerTM 40 .
13
Algorithm Design Environments
Signal Property Meaning nonzero extent of frequency domain center of symmetry in the time domain end of nonzero extent in the time domain end of bandwidth period value is REAL or COMPLEX start of nonzero extent in the time domain start of bandwidth nonzero extent in the time domain multivalued parameter describing timedomain symmetry; possible values are: SYMMETRIC, ANTISYMMETRIC, CONJUGATESYMMETRIC, and
CONJUGATEANTISYMMETRIC
Bandwidth CenterofSymmetry End EndBW Period RealorComplex Start StartBW Support Symmetry
Table 2.1: Signal Properties in SPLICE, ESPLICE and ADE
2.2 Algorithm Design Environments
In his Signal Processing Language and Interactive Computing Environment SPLICE, Myers further generalizes the representation of onedimensional discretetime signals begun by Kopec. SPLICE allows signals to be in nite in extent and to have arbitrary starting and ending points. Besides keeping track of timedomain characteristics of a signal, SPLICE captures information about the signal's bandwidth and other signal properties Table 2.1. In an extended version of SPLICE called ESPLICE, Myers introduces abstract" signals. These are signals that do not have signal values associated with them. Instead, they are de ned completely by their properties. This notion allows ESPLICE to manipulate signals at a higher level of abstraction. ESPLICE also applies abstraction to systems by assigning properties to them see Table 2.2. Unlike signals, however, systems do not have to be instantiated because their properties do not change with di erent values of parameters. For example,
14
Algorithm Design Environments
System Property Meaning ASSOCIATIVE can change grouping of inputs ADDITIVE distributes over addition COMMUTATIVE can change order of inputs HOMOGENEOUS scaled input gives scaled output LINEAR additive and homogeneous MEMORYLESS output does not depend on previous inputs or outputs; if a singleinput system, then SHIFTINVARIANT SHIFTINVARIANT shifted input gives shifted output Table 2.2: System Properties in ESPLICE and ADE shifting by L with respect to n is linear and shiftinvariant regardless of the value of L. Based on the properties of systems, Myers developed rules to simplify and rearrange parts of signal processing expressions. Using these rewrite rules, ESPLICE can generate equivalent forms of algorithms and compare their relative costs by measuring the number of additions and multiplications required. Unfortunately, ESPLICE blindly applies its rewrite rules. As a consequence, a combinatoric explosion in computer memory and computation time can result for complicated expressions. Furthermore, the equivalent forms may loose the regularity parallelism inherent in the original form. ESPLICE set the stage for Covell's Algorithm Design Environment ADE 10 . ADE is a recasting of ESPLICE from Symbolics Zetalisp into Symbolics Common Lisp. Unlike ESPLICE, ADE allows new abstract systems to be de ned. New systems are de ned in terms of their properties. Since rules are already associated with the properties, only rules not based on properties have to be added to handle special cases. Describing systems in terms of their properties prevents a combinatoric explosion in the number of rules required as new systems are added. Compared to ESPLICE, ADE not only encodes new simpli cation and rewrite
15
Algorithm Design Environments
Year 1985 1986 1989 1991 Name SRL ESPLICE ADE Field of Signal System Algorithm Knowledge Simulation Properties Properties Design 1D DSP X X 1D DSP X X X X 1D DSP X X X X Morphology X X X X
MetaMorph
Table 2.3: Overview of Algorithm Design Environments rules, but the engine that applies the rewrite rules is more sophisticated. The engine prunes the number of rewrite rules to apply to an algorithm. For example, the engine would apply a rule to a set of parallel branches in an algorithm only if the parallelism regularity in the branches would be maintained. ADE also uses a more complicated cost function. Instead of being able to count only the number of multiplications and additions, ADE can also count the amount of memory needed and express inequality relationships between individual cost measures, e.g. 1 complex multiply is always more expensive than 3 real multiplies" 41 . Like SRL and SPLICE, ADE can also simulate algorithms. Kopec, Myers, and Covell have focused their work on linear onedimensional signals and systems. Richardson 18, 19, 20, 21 has applied their ideas to the design, analysis, and simulation of morphological algorithms. Morphological signals are sets of points, and morphological systems nonlinearly map one set of points to another. Richardson identi ed signal and system properties based on set theory and encoded a comprehensive collection of simpli cation and rearrangement rules in his MetaMorph environment. To aid in the nding of optimal implementations for morphological algorithms, Richardson developed a new cost metric that takes into account the number of maximum, minimum, union, and intersection operations as well as the number of additions and multiplications. MetaMorph intelligently applies its rewrite rules.
16
Mathematica Programming for KBSP
2.3 Mathematica Programming for KBSP
SPLICE, ADE, and MetaMorph are implemented in Lisp, whereas we have decided to forge an algorithm design environment out of Mathematica. Mathematica 4 , released in 1988, is a general mathematics program like MapleTM 42 and MACSYMATM 43 . All three programs are standardized, documented, portable, programmable environments. Because they embody extensive knowledge about mathematical structures and operations, their forte is the manipulation of formulas. For example, they compute 7x + 3x as 10x. They are pro cient at factoring polynomials and performing partial fraction expansions as well as di erentiating and integrating expressions. Di erentiation and integration are examples of symbolic operations that need to know which symbols to treat as variables. For example, the expression C tanx contains two symbols, C and x. Integrating with respect to x yields ,C logcosx, but integrating with respect to C yields C tanx = 2. Symbolic mathematics programs can also represent and compute sequences from recursive formulas. Since signals can be represented as functions and systems as mathematical operators that map signals to other signals, many signal processing algorithms can be expressed in algebraic form as a sequence of formulas. Formulas can compute signal values point by point or by array operations. In general, symbolic mathematics environments do not support streams or statespace computational models 28 although they can be programmed to do so. All three mathematics programs can implement onetoone transformations as recursive function calls so that wellde ned operations like linear transforms can be conveniently de ned. Alternately, Mathematica and MACSYMA can implement onetoone transformations by applying a list of conditional rules, whereas Maple lacks even the ability to express conditional rules because it does not support pattern matching. Unfortunately, Mathematica and MACSYMA do not have builtin mechanisms to apply rules for transformations that are not onetoone e.g. searching for optimal implementations of an algorithm.
2
17
Mathematica Programming for KBSP
2.3.1 Mathematica
We have chosen Mathematica primarily because of its Notebook front end. The Notebook front end provides pulldown menus, screen and mouseoriented editing of previous commands, and a PostScriptTM driver for all graphics. Graphics, animation, sound, formatted text, and Mathematica code can coexist in the same Notebook. Notebooks typically arrange information as a hierarchical grouping of cells, so a user can quickly and randomly access information. A user can also choose to interact with these active documents by modifying the Mathematica code to try new examples 22 . Although the Notebook front end is a multimedia platform that groups ideas in a tree structure, it is not a true hypertext platform as individual ideas are not interconnected in a web structure 44 . As a multimedia platform, however, the o line documentation can be identical to the online documentation. The Notebook interface exists for several windowing systems e.g. MacintoshTM , MicroSoftTM , NeXTTM , and XTM . Without the Notebook facility, the front end defaults to a simple teletype tty terminal interface. On UnixTM machines, however, the Emacs editor can be con gured to run as a front end to the Mathematica kernel. Although the user interface is machinedependent, the backend kernel is the machineindependent computational engine. The kernel provides a powerful pattern matcher, conditional rule formats, arbitrary precision arithmetic, and a programming language. This programming language organizes code into packages modules much as Lisp does. Embedded in the kernel are packages that perform integration, di erentiation, matrix manipulation, and linear programming. The kernel provides powerful symbolic operations such as di erentiation and integration as well as hundreds of computational functions such as the multidimensional numerical discrete Fourier transform sampled on a rectangular grid. If initialized properly, Mathematica will automatically load new functions as they are invoked by the user 4, 45 . Making the transition to Mathematica from a highlevel programming language is not di cult. Internally, Mathematica maintains expressions in pre x form as Lisp
18
Mathematica Programming for KBSP
mathematical operations relational operators bit operations logical operators
+ =  &&
+= = & 
* =
*= =
++ = ==
!=
Table 2.4: Mathematica Operators That Mimic C Syntax does. Mathematica functions have the form
function argument , argument , : : :
1 2
Mathematica data structures have the same format:
dataType field , field , : : :
1 2
Mathematica's preprocessor enables the user to express function calls in in x or post x format. The preprocessor mimics most of the standard C operators, which are either in in x or post x notation, as shown in Table 2.4. Many of Mathematica's listprocessing primitives such as First and Rest are borrowed from Lisp, as well as the Mathematica primitives Print, Apply, Map, and MapAll. A user can customize the preprocessor so that it can recognize a syntax that is more familiar. The novice can learn more about Mathematica from online documentation written as Mathematica Notebooks or from the more than 20 books written on Mathematica we highly recommend 46 . The help facility built into the kernel, abbreviated as ?, is more useful when users already know the names of the routines they need.
2.3.2 Mathematica As a Platform for Algorithm Design
Mathematica presents an attractive platform upon which to build an environment for algorithm design. As previously mentioned, the kernel provides ways to write programs and encode rules 45 . In fact, many di erent programming paradigms 38 are o ered: procedural, functional, objectoriented, constraintbased, and rulebased.
19
Mathematica Programming for KBSP
The Block and Module constructs are useful for writing procedures, and the Function primitive generates the equivalent of function pointers 47 . Mathematica does not implement objectoriented programming per se but instead allows the programmer to attach data and code to a data type. When using such an datadirected approach, a programmer can write code that is incrementally extensible so that the addition of a new data type tag does not require modi cation of existing code. Constraint propagation, implemented by the Solve command, only applies for equality constraints. Rules in Mathematica are speci ed using an if: : : then form. When the conditional part of the rule is satis ed, its consequence is evaluated. The consequence expression can reference patterns matched during the evaluation of the conditional clause. This rule framework only supports reasoning with certainty. These programming styles o er two direct advantages over highlevel languages for implementing signal processing operations. First, as a consequence of support of the functional programming paradigm, cascaded systems operators become nested calls to Mathematica objects. Second, Mathematica provides several ways to encode rules and apply them to expressions. One way is to attach rules to object heads function names. Mathematica will re such rules whenever an expression contains that function. For example, the relationship
efx + yg = efxg + efyg
can be rendered in Mathematica by either attaching the rule to the Re or Plus builtin primitives:
Re : Plus : Re x + y Re x + y := Re x := Re x + Re y + Re y
The := de nes a procedure in Mathematica, but unlike many highlevel languages, several de nitions can be attached to an object to cover di erent combinations of arguments. The underscore character indicates a pattern to be matched x_ is analogous to ?x in Lispbased pattern matchers 48 . If the rst de nition above is evaluated,
20
Mathematica Programming for KBSP Mathematica's readevaluate loop will check this rule every time the real part of any quantity is sought, even though the rule does not apply. Therefore, the rst but not the second de nition will slow down computations involving the Re function such as taking the absolute value of a complex number because Mathematica's Abs x function is implemented as Sqrt Re x ^2 + Im x ^2 . For any complicated operation, it is desirable for a knowledgebased environment to justify its answers. Students may want to learn how to perform the operation by hand, whereas designers and researchers may want to verify the approach taken by the knowledge base. When the knowledge base evaluates a recursive operation such as onetoone transformations as recursive function calls, the order of operations may follow a preorder traversal of the tree form of the current expression. For example, if the Fibonacci sequence were de ned as
fn = fn, + fn,
1 8
2
f =0
0 7 6 7 7
f =1
1
then tracing calls to f for f would rst show f + f , but then it would show the computations involved in computing f . After f had been computed, then the computations for f would be shown. In this case, displaying the intermediate calculations would only show the transformation of nested subexpressions as they are decomposed i.e. a depth rst traversal of the tree of operations to compute. Such a dialogue would not appear natural or be easy to follow because all of the terms involved in the computation would not be shown at each intermediate step. For the purposes of displaying natural dialogue, a more desirable way to encode recursive operations is to collect related rules in a list called a rule base and apply them recursively e.g. using ReplaceRepeated. Since the rules are not attached to any function name head, this approach avoids slowing down unrelated computations, as was the case when we attached rules to the Re primitive. Another bene t is that a programmer can control how the rules are applied. For onetoone transformations, the programmer can take the expression to be transformed, apply the rst valid rule, and display the intermediate result. By repeating this method until the
6
21
Mathematica Programming for KBSP
transform is computed, the resulting dialogue would mimic how a human would take the transform. This method also allows the user to decide when to re the rules. Our multidimensional signal processing packages MDSPPs implement transforms using this rulebased approach. Being a symbolic mathematics environment, Mathematica does more than apply rules. For example, two primitives partial fraction decomposition Apart and power series expansion Series are critical for taking inverse transforms. The ability to factor, expand, and rearrange polynomials plays another key role in computing symbolic transforms. A programmer can combine these mathematical operations with the rulebased approach to encode powerful routines to perform the ztransform and other onetoone transformations. The user can extend the routines by simply adding transform pairs as rules. Besides implementing a wide variety of programming styles, this environment provides other advantages for algorithm design: 1. dynamic data typing, 2. diverse graphics capabilities, and 3. code generation. Dynamic data typing allows valueless symbols and permits variables to assume any data type, although Mathematica 1.0 generally assumes that unassigned symbols represent a realvalued quantity whereas Mathematica 2.0 generally assumes that unassigned symbols represent a complexvalued quantity. The kernel can graph onedimensional functions and parametric relationships. It also supports scatter, density, and contour plotting as well as threedimensional graphics. On general purpose computers, Mathematica's ability to generate the equivalent Fortran, C, and TEX code for mathematical formulas by default and signal processing expressions by extension becomes useful, because the new code can be spliced into existing programs and documents.
22
Summary
The Mathematica kernel does come with some disadvantages for algorithm design. First, the level of abstraction means that numerical computations are slower, although Version 2.0 introduces the notion of compiled functions to speed this up. Second, Mathematica can only chain rules in one direction it is not an inferencing engine. Third, the rule rewriting mechanism applies a list of production rules in sequential order to transform one expression into another. That is, Mathematica does not generate a tree of all possible new expressions, and therefore, it does not provide heuristic techniques to search through a solution space. Finally, common signal processing operators are, for the most part, not available in Mathematica. The multidimensional signal processing packages overcome the last two de ciencies.
2.4 Summary
This chapter summarizes previous algorithm design environments which have their origins in Kopec's work 28, 29, 30 . Kopec's key contribution was the abstraction of signals as selfcontained objects. In subsequent work, Myers applied abstraction to systems and enhanced the abstraction for signals 7, 8, 9 . Combining abstract signals and systems lead to the ability to simplify and rearrange onedimensional multirate signal processing algorithms and, ultimately, to nd optimal implementations of algorithms based on a simple cost function. Inspired by Myers' work, Covell developed e cient inferencing strategies to apply rearrangement rules in the nding of optimal implementations with the option of maintaining parallelism in the algorithm 10, 11, 12 . Richardson has applied the same concepts to morphological signals and systems in his MetaMorph environment 18, 19, 20, 21 . In developing their algorithm design environments, Myers, Covell, and Richardson had to extend the Lisp language to implement pattern matching, inferencing strategies, mathematical abilities, and objectoriented programming. Because Myers and Covell encoded SPLICE and ADE respectively in Symbolics dialects of Lisp,
23
Summary
SPLICE and ADE are not available for general use. SPLICE, ADE, and MetaMorph do not support graphics and o er only a command line interface. Instead of beginning with Lisp, we have chosen to start at a higher level by using the symbolic algebra program Mathematica. The Mathematica programming language is essentially Lisp without garbage collection but with pattern matching, forward chaining of rules, and a simple implementation of objects built in. The Mathematica kernel adds hundreds of mathematical operations and many plotting routines. Although the Mathematica kernel does not possess much knowledge about signals and systems, we have programmed it to do so. The result is a set of multidimensional signal processing packages MDSPPs for Mathematica which are described in Chapters 4 8. Mathematica is a widely available standardized documented portable environment. On some platforms, Mathematica o ers a sophisticated multimedia Notebook user interface. We have written several Notebooks to accompany the MDSPPs so that the environment is selfdocumenting. Chapter 9 discusses the signal processing Notebooks.
24
CHAPTER 3 Multidimensional Multirate Signal Processing
The ultimate goal of the new multidimensional signal processing packages MDSPPs for Mathematica is to aid in the development of multidimensional multirate signal processing algorithms. At the heart of these algorithms lies the linear periodically timevarying discretetime operations of upsampling and downsampling. In one dimension, upsampling and downsampling as well as their interconnections have been thoroughly researched 49, 50 . The primary purpose of this chapter is to generalize the rules for the onedimensional multirate structures compiled in 49 to multiple dimensions. We also develop the rules for rewriting an upsampler and downsampler in cascade and an upsampler, shifter, and downsampler in cascade. In the case of the up downsampler cascade, we analyze the algorithms underlying both the time and frequencydomain conditions for commutativity. For completeness, we include rules reported by other authors such as polyphase implementations of multidimensional rational decimation systems rate changers 51, 52 . This chapter derives the algorithms underlying these rules in such a way as to be easily implementable by computer. We have encoded the rules and the underlying algorithms in the MDSPPs. In one dimension, many of the multirate rules depend on the commutativity property of multiplication over the semigroup of nonzero integers. Because of commutativity, integers can be factored into a product of primes. In multiple dimensions, resampling is characterized by the semigroup of nonsingular square integer matrices called resampling matrices. The product of resampling matrices rarely commutes, and
25
Lattice and Matrix Theory
as a consequence, factoring must be generalized using left and right matrix factors. The left and right matrix factors, however, are the same for the subset of diagonal resampling matrices. Diagonal resampling matrices correspond to processing multidimensional data separately, one dimension at a time. In simplifying and rearranging cascades of upsampling and downsampling operations, it is essential to transform multidimensional multirate operators into a separable form. In this chapter, we address several fundamental issues associated with the design and implementation of multidimensional multirate systems. Section 3.1 reviews pertinent background material including concepts from lattice theory and matrix theory especially the Smith Form decomposition and the least common right multiple. These mathematical disciplines are needed for a rigorous theory of multidimensional resampling operations. Section 3.2 investigates the onedimensional rules involving cascades of upsamplers and downsamplers that can be generalized to multiple dimensions by using only matrix multiplication and inversion. In Section 3.3, we apply the Smith Form decomposition of a resampling matrix to several problems: computing coset vectors, nding greatest common sublattices, and simplifying up downsampler cascades. Section 3.4 discusses two new sets of equivalent conditions for the commutativity of an up downsampler cascade. Section 3.5 de nes other new rules and collects them all together to form a symbolic algebra for multidimensional multirate signal processing which we have encoded as a part of the multidimensional signal processing packages. We summarize our ndings in Section 3.6.
3.1 Lattice and Matrix Theory
In this section, we present a summary of relevant concepts from the theory of lattices and the theory of matrices. Both are necessary for a rigorous understanding of uniformly sampled or resampled signals. Multidimensional data can be sampled on rectangular, hexagonal, polar, and
26
Lattice and Matrix Theory
other grids 53, 54, 55 . In this chapter, we only consider uniform sampling. Two common uniform sampling grids for image processing are rectangular and hexagonal. Each point on a uniform sampling grid is de ned by an integer vector from the origin to the point. We will refer to the collection of all such integer vectors as RI . RI is an example of a lattice, i.e. a regular arrangement of points in a space. Lattices, sublattices, resampling matrices, and other related topics are de ned in Section 3.1.1. Section 3.1.2 de nes multidimensional upsampling and downsampling. The remaining subsections discuss decomposing and factoring resampling matrices and the relationship of resampling matrices to lattices.
3.1.1 De nitions
The results in this chapter rely on the mathematical concepts of lattices, sublattices, cosets, and coset vectors as they relate to resampling matrices. A lattice is a sampling grid rectangular, hexagonal, etc., a sublattice is a subset of a lattice obtained by an integral linear mapping of a lattice, and cosets represent di erent ways to shift the sampling grid. A brief discussion of these ideas follows. For a more complete discussion of the geometry of numbers, the reader is referred to 6, 56, 57, 58 .
Lattice: A lattice is a discrete set of vectors points, ntuples in Euclidean ndimensional space Rn that forms a group under ordinary vector addition i.e.,
given two points in a lattice, their Euclidean di erence and sum are also in the lattice. A lattice, L, contains an in nite number of points
L = fn j n = u a + : : : + un ang
1 1
3:1
where ui is an integer, and the vectors faig are linearly independent real vectors in an ndimensional space. A lattice can be characterized by a sampling matrix whose columns are the faig vectors. RI will denote the integer lattice.
Resampling Matrix: A resampling matrix is a nonsingular square integer matrix that maps RI into itself. 27
Lattice and Matrix Theory Sublattice: A sublattice is a lattice which is contained within a larger lattice. In this chapter, we will only consider sublattices of the integer lattice RI . A sublattice
associated with the resampling matrix S is de ned as the range of S : sublatticeS fS n j n 2 RI g
Common Sublattice: Given three lattices A, B, and C , if A is a sublattice of B and A is also a sublattice of C , then A is called a common sublattice of B and C. Greatest Common Sublattice: The greatest common sublattice is the union of all
common sublattices.
Cosets and Coset Vectors: A coset is obtained by simply shifting a sublattice by an integer vector, k. Any vector from the origin to a point on a coset is a coset
vector. There are an in nite number of coset vectors.
Distinct Coset Vectors: The distinct coset vectors associated with the resampling matrix S are the j det S j integer vectors that lie within the fundamental parallelepiped FPD of S , denoted by FPDS = fS x 2 RI j 0 xi 1g where xi is the ith element of x. The distinct coset vectors provide a convenient way to enumerate a block of samples for resampling. For example, in onedimensional downsampling upsampling by a factor of N , the distinct coset vectors are the indices 0 : : : N , 1 of each block of input output samples being deleted inserted. The following Theorem and Corollary will be helpful in proving a later result. The proof is omitted here because it can be found in a variety of texts and papers on lattice theory 6, 57 .
28
Lattice and Matrix Theory Theorem 1: The sublattice associated with a resampling matrix A contains the sublattice associated with the resampling matrix B if and only if B = AC where C is also a resampling matrix.
Corollary 1.1: The sublattice associated with a resampling matrix A
is equivalent to the sublattice associated with the resampling matrix B if and only if B = AQ where Q is a resampling matrix such that j det Qj = 1. Corollary 1.1 is true because multiplying on the right by Q such that j det Qj = 1 does not alter the coset vectors. However, multiplying on the left by such a matrix shu es the coset vectors which in turn alters the shape of the sublattice.
3.1.2 Uniform Resampling of Lattices
The fundamental building blocks of multirate systems are the dual operations of upsampling and downsampling.
Upsampling De nition: An upsampler is characterized by a resampling matrix L called the upsampling matrix. The upsampling operation with input x n and output x n is de ned as 8 x L, n if L, n 2 RI x n= 3:2
u 1 1 u
:
and is denoted by " L or "n L or "L; n.
0
otherwise
An upsampler maps the input signal, x n , de ned on RI , to the signal x n which is also de ned on RI . The values of x n are mapped to locations in x n which are on the sublatticeL. All the values of samples of x n not on the sublatticeL are set to zero. The multidimensional Fourier transform which is de ned in 53 relationship for an upsampler is 53, 59
u u u
X ! = X LT !
u
3:3
29
Lattice and Matrix Theory
where ! is a vector of the discretetime frequency variables. In multiple dimensions, a discretetime frequency response X ! is periodic with period 2 in each discretetime frequency variable: X ! = X ! + 2 N r 8 r where r is an integer vector and N is an integer periodicity matrix. The upsampler maps the periodicity by LT , so if the upsampling matrix in nonrectangular, then the periodicity of the output frequency domain would be nonseparable. The sampling density decreases by a factor of j det Lj, the upsampling factor.
Downsampling De nition: A downsampler is characterized by a resampling matrix
M called the downsampling matrix. The downsampling operation, with input x n and output xd n , is de ned as xd n = x M n
and is denoted by M or n M or M; n. A downsampler maps the input signal x n to the output signal xd n by discarding samples of x n that do not lie on sublatticeM and then compacts sublatticeM onto RI . Downsampling introduces aliasing and increases the sampling density by a factor of j det M j, the downsampling factor. The equivalent input output Fourier transform relationship is M 1 j Xj, X
M T , ! , 2k 3:5 Xd ! = j det M j i i where ! is the vector of discretetime frequency variables and ki is a distinct coset vector of M . Note that k is the zero vector Section 3.3.1 shows how to compute distinct coset vectors. The term 2M T , ki for i 6= 0 is called an aliasing vector. Aliasing vectors are periodic with a period of 2 in each dimension. A normalized aliasing vector is an aliasing vector with the 2 term removed so it is periodic with period 1 in each dimension.
det 1 1 =0 0 1
3:4
30
Lattice and Matrix Theory
3.1.3 Regular Unimodular Matrices
In one dimension, upsampling or downsampling by a factor of one either has no e ect on the data resampling by +1 or reverses the data in time resampling by ,1. Resampling factors of one are, however, very useful in higher dimensions. They correspond to resampling matrices with determinants of 1, called regular unimodular resampling matrices 60 .
Unimodular Matrix: A unimodular matrix is a matrix with determinant of 0, +1 or ,1. Regular Unimodular Matrix: A regular unimodular matrix is a matrix with a determinant of +1 or ,1. The inverse of a regular unimodular matrix is also a
regular unimodular matrix.
Regular Unimodular Resampling Matrix: A regular unimodular resampling matrix is a regular unimodular matrix with integer components.
From the input output relationships 3.2 and 3.4 for upsamplers and downsamplers, resampling by a regular unimodular matrix does not change the sampling density but does rearrange the input samples. This permutation of input samples is an essential step in decomposing upsampling and downsampling operations into separable form. Since samples are neither deleted nor added when resampling by a regular unimodular matrix, upsampling by such a matrix is equivalent to downsampling by its inverse and viceversa.
3.1.4 Decomposing A Resampling Matrix Into Smith Form
As mentioned in the introduction, it is often desirable to express a multidimensional multirate operation in separable or decoupled form. One approach is to decompose the resampling matrix underlying the operation. The Smith form of an m n integer matrix S is the matrix product U V , where U and V are square regular unimodular
31
Lattice and Matrix Theory
resampling matrices m m and n n, respectively and is an m n diagonal matrix whose diagonal elements are nonzero integers 6, 59, 60, 61 . Figure 3.1 shows the U , , and V matrices at each step in the computation of the Smith form. In the rst iteration, the element smallest in absolute value is pivoted to the 1; 1 entry and all elements in the rst column except the diagonal entry are reduced modulo the pivot. When all the entries in the rst column other than the diagonal are zero, the process continues with the second column. The rst column is left alone, and the next pivot is found and moved to the 2; 2 entry. The process continues down the diagonal until the integer matrix has been diagonalized. For Smith forms of resampling matrices, all component matrices are square and have the same dimension. By decomposing a downsampling matrix M into its Smith form UM M VM , the downsampling operation becomes equivalent to a shu ing of the input data by UM followed by a separable downsampling by M and a reshu ing of the result by VM . A similar interpretation holds for an upsampling operation. The following properties hold for the Smith form component matrices U , , and V of a resampling matrix S : 1. If S is diagonal, then a valid choice for U and V is the identity matrix. 2. If S is diagonal, then another choice for U is a permutation matrix. V must permute the columns in the same way that U permutes the rows, so V = U , .
1
3. When U = V , , which occurs rarely for nondiagonal resampling matrices, the Smith form corresponds to an eigendecomposition yielding integer eigenvalues the diagonal elements of with eigenvector matrix U .
1
4. If S is symmetric, then U = V T . 5. j det S j = j det j. Because of the last property, the absolute values of the diagonal elements of are integer factors of the absolute value of the determinant of S . If j det S j is a prime
32
Lattice and Matrix Theory
{u, d, v} = SmithNormalForm[ resmat, Dialogue > True ] 1 0 0 1 33 27 24 21 1 0 0 1
pivot location = {2, 2} 1 0 0 1 21 24 27 33 1 0 0 1
modulo reduction by pivot (21) 1 1 1 0 21 3 6 3 1 1 1 0
pivot location = {2, 1} 1 1 1 0 3 21 3 6 1 1 1 0
modulo reduction by pivot (3) 8 7 1 1 3 0 0 15 2 1 1 0
{{{8, 1}, {7, 1}}, {{3, 0}, {0, 15}}, {{2, 1}, {1, 0}}}
Figure 3.1: Iterative Algorithm in the MDSPPs to Compute the Smith Form The original resampling matrix is the middle matrix in the rst two lines of dialogue, i.e. resmat = ff33, 24g, f27, 21gg.
33
Lattice and Matrix Theory
S
Smith Form Smith Canonical Form Alternate Smith Form
2 6 4 2 736 6 256 =4
,982 ,7 ,5 3 2 ,4 0 0 3 2 22 165 47 3 ,421 ,3 ,2 7 6 0 ,180 0 7 6 ,70 ,527 ,150 7 5 4 5 4 5 ,841 ,6 ,4 0 0 24 ,21 ,158 ,45
2 982 6 421 4 2 6 4
3060 1016 864 308 7 5 424 1068 428
3
175 12 73 5 7 5 841 146 10 14 5 3 7 ,1 11 11 5 ,7 ,3 ,2
3
3
2 4 60 4
0 0 12 0 7 5 0 0 360 0 07 24 0 5 0 0 36
3
3
2 6 4 2 6 4
22 165 47 ,1092 ,8221 ,2340 7 5 511 3847 1095
3
2 20 60 4
,884 ,2901 ,1045 3 5271 17311 6234 7 5 ,3558 ,11685 ,4208
Figure 3.2: Examples of Di erent Smith Forms number, then one diagonal element of must be equal to det S and the other diagonal elements must be equal to 1. The Smith form of a matrix is not unique. Figure 3.2 shows three di erent Smith forms for the same resampling matrix. Smith form and Smith canonical form algorithms decompose m n integer matrices in Os and Os elementary arithmetic operations, respectively, where s = m + n + log kS k1 such that kS k1 = max jSi;j j i;j 60, 62 . This algorithm analysis assumes that the intermediate integral calculations t into the machine's integer format. The Smith canonical form, in which each diagonal element of the matrix is a factor of the next diagonal element, is crucial for an e cient implementation of the multidimensional DFT 6 . For the purposes of simplifying and rewriting multidimensional multirate structures, however, the rules are independent of the Smith form algorithm used.
4 5
34
Lattice and Matrix Theory
3.1.5 Finding Common Resampling Matrix Factors
Because the product of resampling matrices does not always commute, resampling matrices have left and right matrix factors as well as left and right matrix multiples 60, 63 . If three resampling matrices A, C , and D exist such that A = CD, then
D is called a right divisor factor of A, C is called a left divisor factor of A, A is called a left multiple of D, and A is called a right multiple of C .
Note that divisor is the term used in the mathematical literature, but for this thesis, factor is usually more descriptive. For multidimensional multirate structures, a usual operation is to nd a greatest common matrix factor and a least common matrix multiple of two resampling matrices. The common right multiple of two matrices A and B is C when 1. C is a right multiple of A, i.e. C = AR for some resampling matrix R, and 2. C is a right multiple of B , i.e. C = BS for some resampling matrix S .
De nition: The least common right multiple LCRM of A and B is a
right multiple whose determinant is less than or equal to the determinants of all other right multiples of A and B in absolute value. The LCRM of two matrices is unique only to within a multiplication on the right by a regular unimodular resampling matrix. Associated with the LCRM is the greatest common left divisor GCLD. A similar de nition holds for the least common left multiple LCLM and its associated greatest common right divisor GCRD.
35
Lattice and Matrix Theory Computing Matrix Factors and Matrix Multiples
The LCRM GCLD and the LCLM GCRD of two matrices can be determined numerically. A method for constructing the LCRM of two n n nonsingular matrices A and B 62, 63 is outlined here. First, form the 2n 2n matrix
2 6 4
A B
0 0
3 7 5
where 0 denotes an n n matrix of all zero entries. Then, by using integral" Gaussian elimination operations, build a regular unimodular square integer matrix regular unimodular is de ned in Section 3.1.3 X of order 2n such that
2 6 4
A B
0 0
32 76 54
X X
11 21
X X
12 22
3 2 7=6 5 4
H 0
0 0
3 7 5
X maps the matrix form formed by concatenating A and B into Hermite normal lower triangular form. From the above matrix equation,
0 = AX + BX C = AX = ,BX
12 12 22
22
3.6 3.7
Here, C is the LCRMA; B and H is the GCLDA; B . A dual algorithm can simultaneously nd the LCLM and GCRD. For two n n matrices A and B , these algorithms require Os elementary arithmetic operations where s = n + log kAk1 + log kB k1 such that kS k1 = max jSi;j j 62 . The LCRM, GCLD, LCLM, and GCRD routines i;j in the MDSPPs bear the same name, e.g. LCRM and so forth.
3
Relationship Between Matrix Factors and Lattices
Associated with the LCRM is the greatest common left divisor. For example, the least common multiple of 2 and 3 is 6 but the greatest common divisor of 2 and 3 is one. Similarly, the LCRM of A and B could be ABV where V is a regular unimodular matrix so the greatest common left divisor is the identity matrix. This is the
36
Resampling Operations in Cascade
multidimensional analog of relative primeness, which is useful for generalizing commutativity conditions for an up downsampler in cascade because it relates to nding the greatest common sublattice of two lattices associated with the up downsampling matrices:
Theorem 2: If C = LCRMA; B , then sublatticeC equals the greatest
common sublattice of sublatticeA and sublatticeB .
Proof: Since C is the LCRMA; B , C = AR = BS for some resampling
matrices R and S . Now applying Theorem 1, we see that sublatticeC is contained in sublatticeA and sublatticeB . For every matrix G that is a common right multiple of A and B and C , sublatticeG is contained in sublatticeA, sublatticeB , and sublatticeC . Now, since the greatest common sublattice of sublatticeA and sublatticeB contains all the common sublattices associated with A and B , the proof is complete.
3.2 Resampling Operations in Cascade
In the design of multidimensional multirate lter banks, a fundamental understanding of the tandem connections of upsamplers and downsamplers is often desirable. These properties are already well known for onedimensional systems 50, 64, 65, 66 . In this section, we simply present the fundamental interconnections of interest. The remaining sections of this chapter develop conditions for simplifying and rearranging these operations. The simplest interconnection to analyze is upsampling by S followed by downsampling by S . Based on the time domain equations, this tandem processing has no e ect on the input signal since S , S = I , As mentioned in the introduction, resampling matrices do not generally commute. Upsampling by L followed by upsampling by L is equivalent to upsampling by
1 1 2
37
Resampling Operations in Cascade

" L0

"K

K

M0

a

" L0

M0

b
Figure 3.3: Equivalent Structures for " L = KL0 Followed by M = KM 0
L L because according to equation 3.2, it is the inverse of the upsampling matrix
2 1
L L , = L, L,
2 1 1 1 1 2
1
that determines how the samples from the input signal are mapped to the output signal so the order of operations is switched. It is also easy to prove that downsampling by M followed by downsampling by M is equivalent to downsampling by M M . Note that in general, the cascade of two upsamplers or two downsamplers is not commutative. They are commutative if and only if the matrix products commute, i.e. L L = L L for upsamplers or M M = M M for downsamplers. Because resampling matrices do not commute, resampling matrices have left and right matrix factors. In one dimension, the common factors in an up downsampling cascade can be removed. In higher dimensions, this condition can be de ned in terms of common left matrix factors. Consider the system shown in Figure 3.3 which downsamples by M = K M 0 after upsampling by L = K L0. Since " KL0 is equivalent to " L0 followed by " K , we see that the system in Figure 3.3a can be simpli ed to the system in Figure 3.3b. In other words, if we can factor the two resampling matrices
1 2 1 2 1 2 2 1 1 2 2 1
38
The Role of Smith Form Matrices
such that they share a common matrix factor on the left, then that common factor can be removed without changing the overall system response.
3.3 The Role of Smith Form Matrices
The Smith form decomposition of a matrix has been useful in many aspects of multidimensional signal processing including the design of multidimensional multirate lter banks 59, 61, 67 and the design of multidimensional DFT and FFT algorithms 6 . Furthermore, the Smith form decomposition has been used in other areas of electrical engineering including control theory for linear systems 39 . In this section, we present new applications of the Smith form decomposition. The rst is the computation of distinct coset vectors Section 3.3.1 which can generate the aliasing vectors associated with a downsampling matrix see Section 3.1.2. Aliasing vectors play a role in determining the commutativity of up downsampler cascades in Section 3.4. The second new application of Smith forms in this section is the nding of the greatest common sublattice GCS of two lattices Section 3.3.2. The GCS leads to rules for simplifying cascades of resampling operations Section 3.3.3. Later in Section 3.4, the GCS will be key in nding a second set of commutativity conditions for up downsamplers in cascade.
3.3.1 Computation of Coset Vectors
The distinct coset vectors of resampling matrix S are those vectors drawn from the origin to the points inside the fundamental parallelepiped of S the points always have nonnegative values. In one dimension, the parallelepiped is the block of samples being reindexed. In the separable multidimensional case, S is diagonal and the fundamental parallelepiped is a rectangular prism whose lower lefthand" corner is the origin and whose upper righthand" corner is p , 1; p , 1; : : : where pi = jSi;ij. Clearly, the ease in enumerating points inside a rectangular prism simpli es the com1 2
39
The Role of Smith Form Matrices
putation of the distinct coset vectors for diagonal resampling matrices. For more general resampling matrices, it is not as easy to compute the distinct coset vectors because it is harder to enumerate the points in the fundamental parallelepiped. The following theorem, o ered by Gardos 68 , handles the case where the resampling matrix S is the product of a regular unimodular matrix U and a diagonal matrix . Although the theorem states that all of the diagonal elements of must be positive, this actually covers the case for all integer diagonal matrices because V can always be adjusted so that has only positive diagonal elements. The proof is given in the appendix of 72 .
Theorem 3: Let U be a unimodular integer matrix and be a diagonal
matrix with all positive elements, then FPDU = U FPD. For a general resampling matrix S , the points inside of the fundamental parallelepiped can be enumerated by rst decomposing S into Smith form U V . Then, if any of the diagonal elements of are negative, multiply on the right and V on the left by the diagonal matrix D such that Di;i = sgni;i. As mentioned in the beginning of this section, it is easy to calculate the distinct coset vectors for because it is diagonal. From Theorem 3, FPDU is just the mapping of FPD by U . From Corollary 1.1, since V is unimodular, sublatticeS = sublatticeU , so FPDS = FPDU mod S = U FPD mod S Therefore, the points in the fundamental parallelepiped of S become the points in the fundamental parallelepiped of mapped by U modulo S : FPDS = fU lS j l 2 FPDg where M is vectormatrix modulo operator de ned as xM = x , M bM , xc
1
3:8
such that the bc operator takes the oor of every component of the vector argument.
40
The Role of Smith Form Matrices
For a given n n resampling matrix S , computing 3.8 requires calculating S , and U . With s = n + log kS k1, computing 3.8 requires one Smith form decomposition Os plus one matrix inversion On plus three matrixvector operations, one integer vector subtraction, and one vector oor operation per point of the fundamental parallelepiped of . Since the fundamental parallelepiped contains j det S j points, equation 3.8 requires Os + n + n j det S j multiplications plus On j det S j oor operations. When integral resampling occurs in every dimension, i.e. when none of the diagonal elements are 1, the n term can be ignored because j det S j 2n n, so the matrix inversion operation is free". The n term can only be ignored when j det S j n which is satis ed when the number of dimensions rises above 7 and integral resampling occurs in every dimension. For image and video processing, the number of dimensions n is 2 or 3 so the n term is signi cant. E cient algorithms to compute cosets that do not use the Smith form are discussed in 69 .
1 4 3 4 3 2 2 3 4 2 4
3.3.2 Computing the Greatest Common Sublattice
The least common right multiple provides one way to compute the greatest common sublattice GCS associated with two resampling matrices see Section 3.1.5. Sometimes, the Smith form decomposition can be used to nd the GCS. The GCS associated with two diagonal resampling matrices is trivial to compute. Namely, given two diagonal resampling matrices and , compute a new diagonal matrix whose ith diagonal element is the greatest common divisor of the ith diagonal element of and . It is easily veri ed that the greatest common sublattice associated with and is simply the sublattice associated with the resampling matrix 0 0 where = 0 and = 0 . For nondiagonal resampling matrices S and S , the GCS can be found when the two resampling matrices share a common left regular unimodular matrix factor U . If the decompositions share a common left factor U = U = U , then according to Corollary 1.1 nding the GCS associated with S and S is equivalent to nding
1 2 1 2 1 2 1 2 1 1 2 2 1 2 1 2 1 2
41
The Role of Smith Form Matrices
the GCS associated with S 0 = U and S 0 = U . This is easily done by noting that U is a onetoone and onto linear operator. Hence, one can simply compute the GCS associated with the resampling matrices and and then apply the inverse mapping U , to generate the GCS associated with S and S . In order to compute the GCS using this approach, we introduce the following algorithm to decompose two resampling matrices S and S into similar Smith forms:
1 1 2 2 1 2 1 1 2 1 2
1. Find a Smith form of S = U V .
1 1 1 1
2. Set U = U
2 2
1 1
3. Find and V from the relationship U , S = V
2 1 2 2
2 1
a
2
2
i;i
is the gcd of the elements of row i of the matrix U , S .
2 2 1
b V can be computed using V = , U , S .
1 1 2
h
i
1
2
4. S and S share a common U factor if V is a regular unimodular resampling matrix
1 2 2
Any Smith form algorithm can be used in the rst step. For n n resampling matrices S and S , the algorithm requires O n + log kS k1 + n elementary arithmetic operations and On gcd operations. The On term arises from the computation of U , S . Using a dual method, a Smith form decomposition can be generated for V = V . In the MDSPPs, the routines SmithFormSameU and SmithFormSameV implement these two algorithms.
1 2 4 3 2 3 1 1 2 1 2
3.3.3 Simpli cation of Resamplers in Cascade
Armed with the ability to compute Smith forms and the GCS, we address the issue of simplying and rearranging upsamplers and downsamplers in cascade. Consider the up downsampler cascade in Figure 3.4a, where we have expressed M = UM M VM and L = ULLVL in their Smith forms. If the decomposition can be performed using
42
Commutativity of Resamplers in Cascade
the same left unimodular matrix UM = UL = U see Section 3.3.2, then it is trivial to remove a factor of U from each operator and the resulting system is shown in Figure 3.4b. Now, downsampling by M follows upsampling by L so common factors between the ith diagonal elements of M and L can be removed thereby creating two new diagonal matrices 0M and 0L shown in Figure 3.4c. Since 0L and 0M are coprime, upsampling by 0L and downsampling by 0M can be switched Figure 3.4d. Figure 3.4e combines the remaining four operations into two by using the fact that downsampling by a regular unimodular matrix is equivalent to upsampling by its inverse and viceversa Section 3.1.3. Hence, Figure 3.4 shows that sometimes the operations in an up downsampler cascade can be switched by modifying the resampling matrices. In a similar way, Figure 3.5 shows that the cascade of a downsampler followed by an upsampler can sometimes be simpli ed and even rearranged.
3.4 Commutativity of Resamplers in Cascade
We now derive two equivalent sets of conditions for the commutativity of downsampling by M and upsampling by L in cascade. In the frequency domain, the conditions are derived by comparing the frequency response of an up downsampler cascade with that of a down upsampler cascade see equations 3.3 and 3.5 in Section 3.1.2. For the twodimensional case, the conditions are 70 :
FD1 The products of the resampling matrices commute, LM = ML, and FD2 The two sets of normalized aliasing vectors N
such that
up down
=
and N
down up
=
are equivalent
N N
up down
=
down up
=
ki; i = 0 : : : j det M j , 1g h i, = fLT M T ki; i = 0 : : : j det M j , 1g
1
= f MT
h
i,1
These conditions were derived without reference to the Smith form decomposition. The Smith form decomposition of the downsampling matrix M enables the enumer
43
Commutativity of Resamplers in Cascade

" VL

" L

" UL

UM

M

VM

a Cascade in Smith form

" VL

" L

M

VM

b Simpli ed cascade if

U
M
=

U
L
" VL

" 0L

0M
VM

c Simpli ed cascade from b with common factors removed

" VL

0M

" 0L

VM

d Reversing the order of operations in c since 0L and 0M are coprime

VL, 0M
1

, " VM 0L
1

e Combining operations in d so that downsampling comes rst
Figure 3.4: Five Equivalent Forms of an Up downsampler Cascade
44
Commutativity of Resamplers in Cascade

UM

M

VM

" VL

" L

" UL

a Cascade in Smith Form

UM

M

" L

" UL

b Simpli ed cascade if

V
M
=

V
L
UM

" L

M
" UL

c Reversing order of operations in b if M and L are coprime

, " LUM
1

, M UL
1

d Combining operations in c
Figure 3.5: Four Equivalent Forms of a Down upsampler Cascade
45
Commutativity of Resamplers in Cascade
ation of the coset vectors ki see Section 3.3.1; for an n n matrix M , generating the coset vectors requires O n + log kM k1 + n + n j det M j multiplications plus On j det M j oor operations. The coset vectors can then be used to calculate the two sets of normalized aliasing vectors N = and N = each having j det M j elements. If the two sets are identical, then the cascade commutes. Thus, the Smith form is the key that unlocks the commutativity conditions based on equivalence in the frequency domain for the general multidimensional case. In the time domain, we derive the conditions by comparing the time response of an up downsampler cascade with that of a down upsampler cascade see equations 3.2 and 3.4 in Section 3.1.2. We discovered the time domain conditions 71, 72 in parallel with 52 . The two resampling operations are commutative if and only if
4 3 2 2 up down down up
TD1 ML = LM , and TD2 LCRMM; L = MLV where V is a regular unimodular matrix. There are at least two ways to test the last condition TD2. The rst is to see if the
matrix product of the inverse of ML and LCRMM; L yields a regular unimodular integer matrix. A second method checks to see if the greatest common left divisor GCLD of M and L is a regular unimodular matrix, which is faster because it does not require the matrix inversion of the rst method. In either method, the LCRM of n n matrices M and L requires Os elementary arithmetic operations where s = n + log kM k1 + log kLk1 62 . A timedomain and a frequencydomain analysis leads to similar pairs of commutativity conditions for an up downsampler cascade, but the timedomain conditions require far fewer computations. Commutativity of matrix products, the rst condition, is easy to check because it only requires On additions and multiplications for n n resampling matrices. Of the second commutativity conditions TD2 and FD2, TD2 the one based on the time domain is much simpler to check because the order of computations is independent of the determinant of the downsampling matrix and varies with n instead of n .
3 3 3 4
46
A Comprehensive Set of Rules
3.5 A Comprehensive Set of Rules
In this section, we collect many rules for multidimensional multirate signal processing. The rules are presented as gures. Figures 3.6 and 3.7 are generalizations of Figures 3.8 and 3.9 in 49 , and Figure 3.8 is a generalization of the interaction between LTI lters and up downsamplers discussed in Section 3.4 of 49 . Figure 3.9 summarizes the rules for cascades of resamplers from Section 3.2 which are adaptations of the identities in Figure 3.10 of 49 . The onedimensional rules governing the commutativity of resamplers in cascade rely on the fact that integer multiplication commutes and that factoring is complete over the integers. In one dimension, two upsamplers or two downsamplers in cascade always commute, but in multiple dimensions, the operations can be switched only if the product of the two resampling matrices commutes, as shown in Figure 3.10a and 3.10b. An up downsampler in one dimension commutes if the resampling factors are relatively prime, i.e. if the least common multiple is the product of the two factors, so this can be easily extended to separable up downsampling where the resampling matrices are diagonal Figure 3.10c. For nonseparable resampling, the two conditions for commutativity derived in Section 3.4 are that the product of the resampling matrices commutes and that the least common right multiple of the two resampling matrices equals the product of the two matrices and a regular unimodular matrix Figure 3.10d. An alternate second condition is that the two sets of aliasing vectors, one for the cascade and one for the commuted cascade, are equivalent Figure 3.10e. The timedomain conditions, however, require far fewer computations.
47
A Comprehensive Set of Rules
s
M
s
z,n0
s
s s

z,M n0
s
M
s
a Delay after downsampling
z,n0
s s
M
s s
M
s
z,M ,1 n0
s
b Delay before downsampling
s
M
s
n
i
s
i s
M
s
"M
n
c Modulation after downsampling
s
n
i s
M
s
s
M
s
i s
M
d Modulation before downsampling
n
Figure 3.6: Identities for Downsamplers
48
A Comprehensive Set of Rules
s
"L
s
z,n0
s
s s

z,L,1 n0
s
"L
s
a Delay after upsampling
z,n0
s s
"L
s s
"L
s
z,L n0
s
b Delay before upsampling
s
"L
s
n
i
s
i s
"L
s
L
n
c Modulation after upsampling
s
n
i s
"L
s
s
"L
s
i s
"L
d Modulation before upsampling
n
Figure 3.7: Identities for Upsamplers
49
A Comprehensive Set of Rules
s
H zM n0
s
M
s
s s
M
s
H zn0
s
a Filtering before downsampling 52, 61
s
"L
s
H zL n0
s
H zn0
s
"L
s
b Filtering after upsampling 52, 61
Figure 3.8: Interactions between Up downsamplers and LTI Filters
In multidimensional signal processing, one typically encounters situations that do not arise in one dimension. This is the case for decomposing matrices into Smith form. The properties of resampling operations expressed in Smith form with component diagonal and regular unimodular matrices are summarized in Figure 3.5. Figure 3.12 shows how the Smith form can be used to remove redundancy if any from up downsampler cascades. Figure 3.13 identi es the conditions for which up downsampler cascades do not commute but yet the operations can be switched by altering the resampling matrices. The last gure, Figure 3.14, shows what happens when a shifter occurs between an upsampler and a downsampler. Remembering that onedimensional downsampling in terms of blocks of samples is analogous to multidimensional downsampling in terms of fundamental parallelepipeds of samples, the shifter will either shift samples o of the downsampling lattice Figure 3.14a or it will delay the incoming samples an integer number of parallelepipeds Figure 3.14b. Figure 3.14c shows a way to decompose a shift operation that exists between an upsampler and downsampler in cascade. In
50
A Comprehensive Set of Rules
s
"S
s
S
s
s s
s
a Up downsampling by the same matrix
s
S
s
"S
s
i s
, pn = 1 for S n 2 RI 0 otherwise
1
b Down upsampling by the same matrix
s
M
1
s
M
2
s
s s
M M
1
2
s
c Cascade of downsamplers
s
"L
1
s
"L
2
s s
"L L
2
1
s
d Cascade of upsamplers
s
" L = ML0
s
M
L M
s
" L0
s
e If 0 = ,1 is a nonsingular integer matrix
L
Figure 3.9: Identities for Cascades of Upsamplers and Downsamplers
51
A Comprehensive Set of Rules
s
M
1
s
M
a if
2
s
s
M
2
s
M
1
s
M1 M2
=
M2 M1
also reported in 52
s
"L
1
s
"L
2
s
s
"L
2
s
"L
1
s
b if
L1 L2
=
L2 L1
also reported in 52
s
" L
s
M
s
s s
M
s
" L
s
c if L and M are diagonal resampling matrices whose corresponding elements are relatively prime
s
"L
s
M
LM V
s
M
s
"L
s
if = and LCRM = such d that is a regular unimodular resampling matrix see Section 3.4; independently reported in 52
ML L; M LM V
s
"L
s
M
s
s
M
s
"L
s
= sets of aliasing vectors e if equivalent and the two 3.4; adapted from 70 are see Section
LM ML
Figure 3.10: Commutativity of Cascades of Upsamplers and Downsamplers
52
A Comprehensive Set of Rules
s
" L = ULLVL
s
s s
" VL
s
" L
s
" UL
s
a The Smith form of an upsampler
s
M = UM M VM
s
UM
s
M
s
VM
s
b The Smith form of a downsampler
s
"U
s
s s
U,
1
s
s
U
s
s s
" U,
1
s
c Resampling by regular unimodular matrices
s
"U
V
s
s s
"V, U
1
s
s
U, V
1
s
d Up downsampling by regular unimodular matrices
s
V
s
"U
s
V U,
1
s
" UV ,
1
s
e Down upsampling by regular unimodular matrices
Figure 3.11: Fundamental Identities Based on the Smith Form Decomposition The Smith form decomposition rewrites a resampling matrix into a matrix product U V where U and V are regular unimodular matrices and is a diagonal matrix see Section 3.3.
53
A Comprehensive Set of Rules
s
" L
s
M
s
s
0L
s
" 0M
s
If L and M are diagonal resampling matrices, then common a factors between each corresponding pair of diagonal elements can be removed.
s
UM M V
s
" ULLV
s
s
UM M
s
" UL L
s
If the resampling matrices can be decomposed into Smith form b using the same right unimodular matrix, then that resampling operation cancels see Section 3.3.3.
s
" U L VL
s
U M VM
s
s
" 0LVL
s
0M VM
s
If the resampling matrices can be decomposed into Smith form left unimodular matrix, then c using the same as well as any common factors that resampling operation cancels between each corresponding pair of diagonal elements see Section 3.3.3.
Figure 3.12: Removing Redundancy in Cascades of Upsamplers and Downsamplers
54
A Comprehensive Set of Rules
one dimension, given two nonzero integers a and b, a family of solutions for integers and in a + b = gcda; b can always be found. An e cient algorithm to nd and rst converts the equation to a Bezout identity by dividing out gcda; b to produce a + ^ = 1 and then enumerates either from 0 to j^j , 1 if j^j jaj ^ b b b ^ or from 0 to jaj , 1 otherwise until an integral solution is found for the other ^ variable. The Bezout identity, solved by the BezoutNumbers routine in the MDSPPs, plays a key role in transforming Smith forms into canonical form 60 . In multiple dimensions, given two square nonsingular integer matrices M and L, a family of integer vectors nM and nL can always be found to satisfy M nM + L nL = n . First, multiply both sides of this equation on the left by the inverse GCRDM; L to convert ^ ^ it into a multidimensional Bezout identity, M nM + L nL = GCRD, M; L n . ^ Next, enumerate either the coset vectors of L as the candidates for nM if j det Lj j det M j or the coset vectors of M as the candidates for nL otherwise until an integer vector solution is found for the other variable. Like the onedimensional Bezout identity, the answers for nM and nL are periodic. The family of solutions can be obtained by shifting nM by L times some integer vector l and by shifting nL by ,M l. In the MDSPPs, the EuclidFactors routine nds nM and nL in the relation M nM + L nL = n for relatively prime M and L. For completeness, we include the polyphase decomposition of a multidimensional rational rate changer 51, 52 . Like the onedimensional version 49 , the polyphase decomposition resamples the lter coe cients to create a more e cient implementation of the rate changer. Just as in the onedimensional case, the polyphase form exists whenever the upsampling and downsampling matrices are relatively prime 52 . The polyphase decomposition is diagramed in Figure 3.15.
0 1 0 0
55
A Comprehensive Set of Rules
s
UM M V
s
" ULLV
s
s
, " LUM
1
s
, M UL
1
s
a Switching operations in a noncommutable down upsampler cascade if the diagonal matrices M and L are relatively prime
s
" U L VL
s
U M VM
s
s
VL, 0M
1
s
, " VM 0L
1
s
b Switching operations in a noncommutable down upsampler cascade by removing the common factors of the corresponding diagonal elements of L and M to produce 0L and 0M .
Figure 3.13: Switching Operations in Noncommutable Cascades
56
A Comprehensive Set of Rules
s
"S
s
zn0
s
S
s
s
0
z
s
a Up downsampling by when the shift vector n0 62 sublatticeS or equivalently when ,1 n0 is not an integer vector
S S
s
"S
s
zn0
s
S
s
S,1 n0
s
b Up downsampling by when the shift vector n0 2 sublatticeS or equivalently when ,1 n0 is an integer vector
S S
s
"L
s
zn0
s
M "L
s s
M
s
s
zn
L
zn
L
M
s
c For resampling matrices L and M, n0 can always be rewritten as n0 = nL + nM using a form of Euclid's algorithm or equivalently by converting the equation to a Bezout identity
M
Figure 3.14: Interaction between Up downsampler Cascades and Shifters
57
A Comprehensive Set of Rules

"L

H z

M

a Multidimensional Rational Decimation System

"L

z,k0

M

E0 z


z,k0 z,k1
z,k1

M

E1 z

M
z,Lk02

z,M k01

z,Lk12 z,Lk22

z,M k11 z,M k21
z,k2
M
E2 z

z,k2
L
b Rewriting the Downsampling Stage
z,k02

c Rewriting the Shift Operations if and are relatively prime on the left

"L

M

E0z

z,k01

"L

M
z,k02

"L

M

E1z

z,k01

M

z,k02
"L
M
E2z
z,k01

"L
ML
d Moving the Shift Operations Outside of the Up Downsamplers
z,k02

e Switching the Order of the Up Downsamplers if = and and are relatively prime
LM L M

M

E0;0z

"L

zj0

z,k01


E0;1z

"L

zj1


z,k12

M

E1;0z E1;1z

"L "L

zj0

z,k11

zj1


z,k22

M

E2;0z E2;1z

"L "L

zj0

z,k21

zj1
f Final Polyphase Form
Figure 3.15: Polyphase Implementations of Rational Decimation Systems In the end, all ltering operations are performed at the lowest sampling rate in the system adapted from 51 .
58
Summary
3.6 Summary
As shown in Figures 3.6 through 3.14, this chapter presents original research that generalizes the simpli cation and rearrangement rules for onedimensional multirate signal processing in 49 to multiple dimensions. This chapter also generalizes two other important onedimensional rules to multiple dimensions 1 the commutativity of an upsampler and downsampler in cascade, and 2 pulling out the shift operation in the cascade of an upsampler, shifter, and downsampler. Vaidyanathan and Chen 52, 61 independently derived these two rules in multiple dimensions and then applied them to nd polyphase implementations of multidimensional rate changers i.e. rational decimation systems. We include their polyphase rules for completeness. This chapter also develops the algorithms necessary to evaluate the conditions of the rules in Figures 3.9, 3.10, 3.12, and 3.13. One of new algorithms computes the greatest common sublattice by computing common left matrix factors using on Smith form decompositions. Its dual, which nds common right matrix factors, is useful in removing redundant operations and switching operations in noncommutable cascades. The new algorithm underlying Figure 3.10e uses the Smith form decomposition of the downsampling matrix to compute the two sets of aliasing vectors. All of the supporting algorithms and all of the rules in Figures 3.6 through 3.14 have been encoded in the MDSPPs. LatticeTheory.m" contains the algorithms and RewriteRules.m" contains the rules. A version of LatticeTheory.m" is available that does not rely on other routines in the MDSPPs. Also encoded in our lattice theory package is a family of Smith form algorithms that nd matrices whose diagonal elements are as close to one another in absolute value as possible and or minimize the norm of either U or V 27 . Other researchers have used this new family of Smith form algorithms in our algorithm to compute common left matrix factors as a part of a procedure to design multidimensional nonuniform lter banks 73 .
59
CHAPTER 4 The Signal Processing Packages: Implementing Linear Systems Theory in
Mathematica
The multidimensional signal processing packages MDSPPs 1, 2, 3 represent a signi cant step in implementing linear systems theory in a computer algebra program Mathematica. Theory refers to the algebraic and other symbolic operations used in studying the behavior of linear systems. The MDSPPs obey Mathematica's syntax rules and follow engineering convention when possible. The MDSPPs de ne many new signals Section 4.1 and systems Section 4.2 that are missing in the core of Mathematica. Section 4.3 and 4.4 introduce some of the plotting and symbolic abilities of the MDSPPs. Chapter 5 covers more advanced features of the MDSPPs.
4.1 De ning Signals
The MDSPPs de ne the signals listed in Table 4.1. The new signals are functions that take one argument, except for Pulse and CPulse which require two the pulse length and the variable plus the o set if any. Onedimensional signals require a variable or expression to specify the domain of the signal. For example, as a function of t, the continuous step function u, would take the form of Unit 1 t or more simply CStep t which stands for Continuous Step". The expression CStep t2 would represent a step function delayed by 2 or u, t , 2. Signals can appear in algebraic expressions, such as t CStep t , which is one way to write the ramp function.
1 1
60
De ning Signals
Object
CPulse
CStep
Delta
Dirichlet
Impulse
Pulse
Sinc Step
Unit
Meaning continuous pulse: 8 0 t 0 1 t=0 2 uLt = 1 0 t L 1 t=L 2 : 0 t L continuous step function: 8 1 t 0 1 u,1 t = t=0 2 : 0 t 0 Dirac delta function t which has a value of zero for t 6= 0, a value of in nity at t = 0, and unit area: Z1 t dt = 1 ,1 normalized Dirichlet discrete kernel: sinN !=2 d N; ! = N sin!=2 Also called the ASinc function Kronecker delta function: 1 n=0 n = 0 n 6= 0 discrete pulse: 8 0 n 0 uL n = 1 0 n L , 1 : 0 nL sinct sint=t so sinc0 = 1 discrete step function: 1 n0 un = 0 n 0 family of functions which includes CStep and Delta
Table 4.1: Signals Introduced by the Signal Processing Packages
61
De ning Systems
Separable multidimensional signals are simply products of onedimensional functions. The positive quadrant for the discretetime variables n and n would be written as u n u n which is rendered as Step n1 Step n2 in Mathematica. Nonseparable multidimensional signals are easy to represent as well. For example, the line impulse through the origin running at a 45 degree angle in the continuoustime variables t and t would be written as t , t which is rendered simply as Delta t1  t2 in Mathematica.
1 2 1 2 1 2 1 2
4.2 De ning Systems
The MDSPPs de ne the operators commonly used in signal processing but missing in Mathematica. Table 4.2 lists the new operators and their parameters. Parameters are side information such as shift factors, DFT lengths, and the variables on which to operate. Operators also take arguments i.e. signals:
operator parameter , parameter , : : :
1 2
signal , signal , : : :
1 2
Cascaded systems are formed by nesting operators. For example, with Times as the builtin product operator which takes no parameters, the expression
Times Exp 2 I Pi t 10 , Shift 10,t x t
represents the modulation by exp,2jt=10 of xt , 10 which is in turn represented by the Shift operator note that a space between terms denotes multiplication. All operators in Table 4.2 represent operations in the sense that they defer evaluation of operations until TheFunction or N is applied to them. That is, when an operator from Table 4.2 is applied, no evaluation takes place instead, the resulting function is stored symbolically, until it becomes convenient to compute it explicitly. For onedimensional operators, a variable is a symbol, like t. Arguments other than variables can be numbers, symbols or formulas. To indicate that the operator is multidimensional, we use a list of variables such as t1, t2, t3 instead of a single
62
De ning Systems
Aliasby sc, w CircularShift n0, N, n CConvolve t Convolve n DFT L, n, k
DTFT n, w Difference i, n Downsample m, n
FIR n,
fh0,
h1, : : :
g
FT t, w IIR n,
fa0,
a1, : : :
g
Interleave n InvDFT L, k, n InvDTFT w, n InvFT w, t InvL s, t InvZ z, n L t, s Periodic p, v Rev v ScaleAxis sc, w Shift v0, v Summation i, ib, ie, inc Upsample l, n
Z n, z
aliases a continuous function of w giving it a period of 2 sc and divides by sc shifts a sequence by n0 + n mod N in n continuous convolution in t discrete convolution in n the Lsample in 1D and the j det Ljsample in mD discrete Fourier transform of a function in n to a function in k the discretetime Fourier transform of a sequence in n to a continuous periodic function in w the ith backward di erence in n with the downsampling factor f equal to m in 1D or j det mj in mD, operator keeps the rst sample in every block of f samples; the sampling rate decreases by a factor of f allzero digital lter with nite impulse response; in 1D, the weights taps are h0, h1, : : : and in mD, the weights become a mD volume of weights continuous Fourier transform of a function of t into a function of w allpole digital lter with in nite impulse response; in 1D, the feedback coe cients are a1, : : : , and in mD, fa0, a1, : : : g becomes a mD volume of coe cients interleaves combines samples from each input function of n into one function inverse discrete Fourier transform see DFT inverse discretetime Fourier transform see DTFT inverse continuous Fourier transform see FT inverse Laplace transform see L inverse z transform see Z Laplace transform of a function of t argument is made periodic with period p with respect to variables v; in mD, p is an integer matrix reverse the function of variable v ip the v axis scale the w axis variable by sc shift function of v by v0 points summation operator: i = ib to ie step inc with the upsampling factor f equal to l in 1D or j det lj in mD, operator inserts f , 1 samples after each input sample; sampling rate increases by a factor of f z transform of a function of n
Table 4.2: New Operators and Their Parameters in the MDSPPs
63
Plotting Signals and Systems
variable such as t. The length of the variable list indicates the dimensionality of the operator. In multiple dimensions, the rst parameter for the upsampling, downsampling, and periodic operators becomes a square matrix, whereas all other parameters become lists of expressions, one expression per dimension. For example, the downsample operator requires a downsampling constant m and a symbolic variable n. In onedimensional downsampling, m is the integer downsampling factor and n is a symbolic index, but for the multidimensional case, m is the integer downsampling matrix and n is a list of symbolic indices. Regardless of the dimensionality, the form of the operator is Downsample m, n , and its syntax is Downsample m, n f , where f is a discrete function in n.
4.3 Plotting Signals and Systems
All of the continuoustime signals introduced by the MDSPPs, except for the Dirac delta function, can be graphed by Mathematica's basic plotting command Plot. Plot, however, will experience problems if the signal contains any of the new operators introduced by the MDSPPs, or if the signal contains Dirac delta functions, or if the signal is complexvalued. Thus, we introduce a new function called SignalPlot. SignalPlot will rst convert the expression into functions which can be plotted and then graphs the real and imaginary parts of the function in solid and dashed lines, respectively. This plotting convention for analog signals follows Bracewell 74 . SignalPlot also plots Dirac delta functions as upward or downwardpointing arrows. The height of the arrows is either the area of the delta function when $DeltaFunctionScaling is set to Scaled or the same height which is determined from the values of the rest of the plot when $DeltaFunctionScaling is set to None. For discretetime signals, we have introduced the function SequencePlot that only plots realvalued functions of an integral index variable.
64
Continuous and Discrete Convolution
Both new plotting routines can plot onedimensional and twodimensional signals. For twodimensional sequences, SequencePlot generates threedimensional graphics initially sampled at discrete indices. Sometimes, a better way to visualize a twodimensional sequence is to view it at an in nite height above the plane of the two independent variables, known as a density plot. A density plot, for example, shows clearly the sampling pattern after multidimensional upsampling. For example, Figure 4.1 shows how to use Mathematica's DensityPlot command with our upsampling operator to illustrate the e 2 of upsampling a separable sequence by ect 3 1 15 the nonrectangular upsampling matrix 4 , which has an upsampling factor ,1 1 of 1 1 , ,1 1 = 2. First, the upsampling changes the basis of the sampling lattice to a nonrectangular lattice so that the function becomes nonseparable. Second, the upsampling inserts zeros in a checkerboard pattern as seen by the black tiles which represent an amplitude of zero. This checkerboard pattern corresponds to interleaving in digital video signals, and the resampling matrix involved belongs to the family of quincunx resampling matrices. In Figure 4.1, the white tiles represent an amplitude of two.
4.4 Continuous and Discrete Convolution
Convolution is a key concept in linear shiftinvariant systems as it describes ltering in the time domain. Continuous convolution of two functions, f t and gt, is de ned as 5 Z1 f t ? gt = ,1 f gt , d In discrete time, the integral is replaced by summation 75
fn ?gn =
1 X
The graphical approach to computing onedimensional convolutions ips g in time and slides it across f . The answer is obtained by performing the integration or summation over those intervals where the functions overlap.
m=,1
f m g n,m
65
Continuous and Discrete Convolution
sin2D = ( Sin[n1 Pi / 5] + 1 ) * ( Sin[n2 Pi / 5] + 1 ); DensityPlot[ sin2D, {n1, 11, 11}, {n2, 11, 11}, PlotPoints > {23, 23}, PlotRange > All ]
10
upMatrix = {{1, 1}, {1, 1}}; upSignal = Upsample[upMatrix, {n1,n2}][sin2D]; DensityPlot[ Evaluate[TheFunction[upSignal]], {n1, 11, 11}, {n2, 11, 11}, PlotPoints > {23, 23}, PlotRange > All ]
10
5
5
0
0
5
5
10 10 5 0 5 10
10 10 5 0 5 10
DensityGraphics
DensityGraphics
Figure 4.1: Visualizing TwoDimensional Sequences as Density Plots Top views of a separable function left and the upsampling of the same function by a nonrectangular matrix right.
66
Continuous and Discrete Convolution
The MDSPPs introduce PiecewiseConvolution which performs discrete convolution and convolution for piecewise continuous functions. PiecewiseConvolution accepts f and g as algebraic expressions or as lists of Fintervals. West and McClellan 26 de ne the Finterval as a list of three elements:
ffunction, leftendpoint, rightendpointg
Fintervals can be nite or in nite in extent, but when given numeric values, the right endpoint should always be greater than the left endpoint. The function can be any algebraic expression, e.g. a t^2. As a list of Fintervals, then, t CPulse 10, t + 5 Delta t  20 becomes f ft, 0, 10g, fArea 5 , 20, 20g g. Internally, PiecewiseConvolution uses the following algorithm to compute convolutions 26 : 1. convert f and g to lists of Fintervals, 2. convolve each Finterval in f with each Finterval in g, and 3. simplify the result by a reordering sorting the Fintervals and b merging overlapping Fintervals when possible. West 26 bases the computation of the convolution of two Fintervals on the Square Matrix Rule" 76 which says that convolving two nite extent intervals produces three nite extent intervals. We have generalized the rule to include situations where one or more of the endpoints is in nite as shown in Table 4.3. Since each left endpoint in an interval can either be ,1 or a number and each right endpoint in an interval can either be a number or 1, convolving two Fintervals yields 2 2 2 2 = 16 possible combinations of intervals. All of the combinations yield one, two, or three intervals. We have enumerated all possible combinations of intervals in validating Table 4.3.
67
Continuous and Discrete Convolution
Interval t 2 ,1; lf + lg t 2 lf + lg ; lf + ug Value 0 f gt , d Valid always if lf 6= ,1 and lg 6= ,1 if lg 6= ,1 and ug 6= 1 if uf 6= 1 and ug 6= 1 always Valid always if lf 6= ,1 and lg 6= ,1 if lg 6= ,1 and ug 6= 1 if uf 6= 1 and ug 6= 1 always
Z t,l Zl t,l
f
g
t 2 lf + ug ; uf + lg t 2 uf + lg ; uf + ug t 2 uf + ug ; 1
Interval n 2 ,1; lf + lg , 1 n 2 lf + lg; lf + ug , 1
g
Zt,u u
f
g
t,u
g
f gt , d f gt , d 0
Value 0 f m g n,m
a Formulas for continuoustime convolution
n,l X m=l n,l X
f
g
f
n 2 lf + ug; uf + lg n 2 uf + lg + 1; uf + ug n 2 uf + ug + 1; 1
g
m=n,u u X m=n,u
f m g n,m f m g n,m 0
g
g
b Formulas for discretetime convolution I. if lf = ,1 and uf = 1 OR lg = ,1 and ug = 1, then the convolution extends for all time and is computed by the de nition. II. if lf = ,1 and uf 6= 1 and lg 6= ,1 and ug = 1, then the convolution extends for all time with two intervals of interest for the output variable v: A. for v 2 ,1; uf , apply de nition with upper limit being v, and B. for v 2 uf ; 1, apply de nition with upper limit being uf . III. if lf 6= ,1 and uf = 1 and lg = ,1 and ug 6= 1, then apply II above after switching the two Fintervals. Table 4.3: Modi cation to the Square Matrix Rule for Convolution Convolving one Finterval ff , lf , uf g with another Finterval fg, lg, ugg generates one, two, or three intervals with nonzero functions. The tables assume that the rst Finterval is the longer one, i.e. uf , lf ug , lg . The three special cases are enumerated above.
68
Summary
Even though PiecewiseConvolution does not use the ipandslide approach to calculate outputs, the routine does have the ability to demonstrate the ipandslide approach by means of animation, e.g. convolving two pulse functions in t:
PiecewiseConvolution CPulse 1,t , CPulse 1,t , t, Dialogue All
Enabling the Dialogue option triggers the generation of the animation sequence: False for no animation, True for the ipandslide of the rst function across the second, and All for the ipandslide of the rst function across the second as well as the corresponding overlapping intervals for each ipandslide frame. The Dialogue option is used in many other routines in the MDSPPs, as is demonstrated in Chapter 5.
4.5 Summary
Mathematica possesses some inherent signal processing capabilities. It can compute the numerical multidimensional Fourier transform of a multidimensional sequence of data. It de nes many functions such as trigonometric, exponential, logarithmic, and Bessel functions commonly used in an algebraic representation of signals. However, basic signals such as step and impulse functions and basic systems such as shifters and lters are missing in the core of Mathematica. By extending Mathematica, the MDSPPs de ne many of the missing signals and systems. This chapter also introduces some of the plotting and symbolic abilities of the MDSPPs. Examples of plotting in the time" domain include signal plots, sequence plots, and density plots for the visualization of resampling operations. The lone symbolic ability introduced by this chapter is linear convolution which is implemented for both the discrete and continuous domains. Chapter 5 describes the MDSPPs in more detail.
69
CHAPTER 5 Analysis of Multidimensional Signals and Systems
Chapter 2 introduces SPLICE and its descendants which can deduce certain properties of onedimensional discretetime signals generated by various components of a system. The multidimensional signal processing packages MDSPPs extend the ability of SPLICE to analyze signal properties to multiple dimensions. As discussed in Section 5.1, several new important properties have also been added. Besides generalizing the analysis of signal properties to multiple dimensions, the MDSPPs de ne a suite of analysis tools for multidimensional discretetime and continuoustime signals based on linear transform theory. For continuoustime signals, the MDSPPs implement the multidimensional Laplace and Fourier transforms. For discretetime signals, the MDSPPs implement the multidimensional z, discretetime Fourier, and discrete Fourier transforms. The linear transform routines are symbolic because they map algebraic expressions between the time" and various frequency" domains see Section 5.2. Analysis of separable multidimensional systems reduces to analysis of each dimension separately. The onedimensional capabilities of the multidimensional z and Laplace transform routines are useful in solving linear constant coe cient di erence and di erential equations, e.g. those describing passive networks see Section 5.3.1. Coupling the symbolic linear transforms with graphical representations of signals enables new routines to automate the analysis of onedimensional signals in the time, generalized frequency z or Laplace, and Fourier frequency domains see Section 5.3.2.
70
Extending Signal Properties in ESPLICE and ADE
The linear transform routines are extensible in that users can provide their own transform pairs. Specifying transform pairs for input signals enables the transform routines to generate the inputoutput relationships of multirate systems see Section 5.4. The forward z and Laplace transform routines support multidimensional multilateral signals, so they must track the region of convergence ROC. The ROC can be used to determine algebraic conditions for the stability of multidimensional systems see Section 5.5.1. After assigning values to free parameters, stability can be determined graphically for twodimensional systems see Section 5.5.2. The multidimensional capabilities of the transform routines assist in the automation of signal analysis and the derivation of transform properties. Automated twodimensional signal analysis Section 5.6.1 is an extension of the general signal analyzers described earlier. Deriving properties of linear transforms including tracking the changes in the region of convergence is now possible in higher dimensions Section 5.6.3.
5.1 Extending Signal Properties in ESPLICE and ADE
Myers 7, 8, 9 identi ed a set of properties for onedimensional signals as shown earlier in Table 2.1. Each of these properties can be generalized to the multidimensional case. The Period property is now described by a nonsingular integer periodicity matrix which is diagonal for separable signals. We have made StartBW and EndBW more descriptive by renaming them to StartBandwidth and EndBandwidth, respectively. The properties of Start, End, StartBandwidth and EndBandwidth become sets of coordinates. Similarly, Support and Bandwidth become lists with one element per dimension. Together, these two slots pieces of data de ne the smallest nonzero rectangular region of support containing the time" and frequency"
71
Linear Transforms
domains, respectively. We add the methods procedures InDomainFunction and InBandwidthFunction to de ne the shape of the region of support. For example, a circular region of support for the general multidimensional case would be written as
Function coordinates, Apply Plus, coordinates^2
2 1 2 2 2 1 2
r^2
which computes x + x + r where x ; x ; : : : is the list of coordinates. The methods InDomainQ and InBandwidthQ check rst if the passed set of coordinates is in the rectangular region of support and then call InDomainFunction or InBandwidthFunction to determine if the set of coordinates is really in the region of support. Many types of symmetry can exist for di erent combinations of the signal variables, so the properties Symmetry and CenterOfSymmetry become lists of the types of and centers of symmetry, respectively. De ning the types of symmmetry and their centers, however, in multiple dimensions becomes very di cult. The RealorComplex property has been replaced by DataType which can take the value of Real, Complex, Integer, etc. The MDSPPs add the property Variables to keep track of the time" domain variables. Table 5.1 lists the multidimensional signal properties recognized by the MDSPPs.
5.2 Linear Transforms
Many signal processing algorithms are ultimately based on linear transforms z, discretetime Fourier, discrete Fourier for discretetime signals and Laplace and Fourier transforms for continuoustime signals. The multidimensional signal processing packages, which are introduced in Chapter 4, implement these transforms in their most general multidimensional multisided symbolic forms by applying a sequence of transformation rules. Each routine rst checks special multidimensional rules and then applies the onedimensional rule base to each dimension of the signal. Even though we use the fact that the linear transform kernels are separable, we do not make any
72
Linear Transforms
Signal Property
Bandwidth CenterOfSymmetry DataType End EndBandwidth InBandwidthFunction InBandwidthQ
InDomainFunction InDomainQ Period Start StartBandwidth Support Symmetry Variables
Meaning nonzero rectangular extent of frequency domain center of symmetry in the time domain value is Real, Complex, Integer, : : : end of nonzero rectangular extent in the time domain end of nonzero rectangular extent of bandwidth function de ning the shape of the frequency domain checks to see if the coordinates are in the region of support for the time domain function de ning the shape of the time domain checks to see if the coordinates are in the region of support for the frequency domain multidimensional period start of nonzero rectangular extent in the time domain start of nonzero rectangular extent of bandwidth nonzero rectangular extent in the time domain multivalued parameter describing timedomain symmetry; possible values are: Symmetric, AntiSymmetric, ConjugateSymmetric, and ConjugateAntiSymmetric symbols used as time domain variables
new multidimensional property unnecessary in one dimension
Table 5.1: Signal Properties Supported by the Signal Processing Packages
73
Analysis of OneDimensional Systems Using Transforms
assumptions about the signal being transformed. The rules in each onedimensional rule base are applied sequentially. At each iteration, the rst rule that applies to the expression being transformed is invoked thereby rewriting the expression. The process repeats until the expression no longer changes. Since the rules are applied in sequential order, the order of the rules is very important. We have encoded the rules in the order that humans tend to use them: transform pairs, transform properties e.g. additivity, transforms of operators e.g. shift, and transform strategies e.g. partial fraction decomposition. Figure 5.1 diagrams the general structure of the onedimensional rule bases. All of the inverse transform routines return a nal algebraic expression, whereas the forward transform routines return a data structure. The forward z and Laplace transforms track and return the region of convergence ROC. The ROC is necessary to guarantee uniqueness because the z and Laplace transforms are multilateral. All of the linear transforms can display the intermediate calculations so that designers can check the computation. Users can also specify their own transform pairs. This feature is useful in generating transfer functions for systems as shown in Section 5.4 and deriving properties of transforms as shown in Section 5.6.3.
5.3 Analysis of OneDimensional Systems Using Transforms
5.3.1 Solving Di erence and Di erential Equations
The linear di erential equation solver LSolve uses the unilateral Laplace transform to solve onedimensional linear constant coe cient di erential equations. Figure 5.2 shows the steps taken by LSolve in solving for yt in the equation
y00t + 3 y0t + 1 yt = eatu, t 2 2
1
74
Analysis of OneDimensional Systems Using Transforms
START
Break up nonseparable functions
Known rational transform pair?
NO
YES
YES
Known nonrational transform pair?
NO
Do any properties apply?
NO
YES
Known operator?
NO
YES
Do any strategies apply?
NO
YES
Apply the catchall strategy
STOP
Figure 5.1: Structure of the OneDimensional Transform Rule Bases Onedimensional transform rule bases form the backbone of the more general multidimensional, multisided transforms. Each rule base calls a onedimensional rule base once for each dimension of the transform.
75
Analysis of OneDimensional Systems Using Transforms
In[17]:=
Solving for the unknown transform yields 1 3 y[0]  +  + s y[0] + y’[0] s  a 2 1 3 s 2  +  + s 2 2 which becomes 1 10 + 4 s + s  a 1 3 s 2  +  + s 2 2 Inverse transforming this gives y[t]:
Out[17]=
LSolve[ y’’[t] + 3/2 y’[t] + 1/2 y[t] == Exp[a t] CStep[t], y[t], y[0] > 4, y’[0] > 4, Dialogue > All ] Solving for y[t] in the differential equation y[t] 3 y’[t] a t  +  + y’’[t] = E CStep[t] 2 2 subject to the initial conditions {y[0] > 4, y’[0] > 4}
The Laplace transform of the left side is: 1 3 s 2 3 y[0] L { y[t] } ( +  + s )   2 2 2 s y[0]  y’[0] ( In the general case, the unilateral Laplace transform of the nth derivative of y[t] is: (n) L {y [t] } = L { y[t] } s n +
2 t/2 t/2 (2 (5  16 a  12 a + 6 E + 22 a E + 2 t/2 t + a t 16 a E + E ) CStep[t]) / 2 t ((1 + 3 a + 2 a ) E )
In[18]:=
1 + n 2 + n s y[t0] s y’[t0] () + () + ... + s t0 s t0 E E (1 + n) y [t0] () s t0 E where t=t0 is the initial condition. )
yoft = % /. ( CStep[t] > 1 )
Out[18]=
2 t/2 t/2 (2 (5  16 a  12 a + 6 E + 22 a E + 2 t/2 t + a t 16 a E + E )) / 2 t ((1 + 3 a + 2 a ) E )
In[19]:=
Simplify[ yoft /. t > 0 ]
Out[19]=
4
In[20]:=
The Laplace transform of the right side is: 1 s  a
Simplify[ D[ yoft, t ] /. t > 0 ]
Out[20]=
4
Figure 5.2: Interaction with the Di erential Equation Solver
76
Analysis of OneDimensional Systems Using Transforms
subject to the initial conditions y0 = 4 and y00 = 4. Because full justi cation is enabled by the option Dialogue  All, the solver describes the intermediate calculations. The user has access to the answer as well as how to obtain it. Once the answer has been obtained, the user can check the answer. Command 18 de nes the variable yoft as yt for t 0; since u, t = 1 for t 0, the continuoustime step function is replaced by 1. Command 19 veri es that y0+ is indeed 4 and the last command checks the value of y00+. Like LSolve, ZSolve uses the unilateral ztransform to solve onedimensional initial value linear constant coe cient di erence equations. Figure 5.3 shows the steps taken by the di erence equation solver to compute the closedform formula for the wellknown Fibonacci sequence:
+ + 1
yn , yn, , yn, = 0
1 2
y =0
0
y =0
1
The solver computes the nonrecursive formula for yn in exact precision. Even though the solution contains irrational numbers, the formula yields integers for the Fibonacci numbers because arithmetic is performed in exact precision. The rst 11 values of the sequence are computed by the Table command in Figure 5.3.
5.3.2 Generalized Signal Analysis
The MDSPPs support new formats for plotting signals, sequences, and their transforms. New routines plot 1D and 2D signals and systems in the time domain as well as their polezero diagrams and frequency responses. They can generate several di erent styles of frequency plots, i.e. whether the frequency domain scale is linear or logarithmic and whether the magnitude plot has a linear or decibel range. They also graph root loci for one varying parameter. Two routines, ASPAnalyze and DSPAnalyze, combine these plotting abilities with the transform capabilities to perform a complete analysis of 1D and 2D signals. These analyzers display textual and graphical information. The textual dialogue
77
Analysis of OneDimensional Systems Using Transforms
fib[n_] = ZSolve[ y[n]  y[n1]  y[n2] == 0, y[n], y[0] > 0, y[1] > 1, Dialogue > All ] Solving for y[n] in the difference equation y[2 + n]  y[1 + n] + y[n] = 0 such that n > 1 and subject to the initial conditions {y[0] > 0, y[1] > 1}
The ztransform of the left side is: 2 (1  z 1  ) Z{ y[n] } z
The ztransform of the right side is: 1 z Solving for the unknown ztransform yields 1 2 1 (1  z  ) z z Inverse transforming this gives y[n]: 1 + n 1 1 + n 2 2 () Step[1 + n] 1  Sqrt[5]  + Sqrt[5] (1  Sqrt[5]) 1 + n 1 1 + n 2 2 () Step[1 + n] 1 + Sqrt[5] Sqrt[5] (1 + Sqrt[5])
The first step is to substitute the initial conditions into the difference equation for n from 0 to 1. This will give rise
to 2 impulse function(s) which will be added to the righthand side of the difference equation. Then, the unilateral ztransform will be used to solve the difference equation in the index n.
Solving for y[n] in the difference equation augmented by the initial conditions: y[2 + n]  y[1 + n] + y[n] = Impulse[1 + n]
Table[ Simplify[fib[i]], {i, 0, 10} ] {0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55}
Figure 5.3: Solving the Fibonacci Di erence Equation
78
Analysis of Multirate Systems Using Transforms
includes the Fourier and the generalized transform z or Laplace as well as the algebraic conditions for stability in terms of the free parameters. Graphics information includes the timedomain plot, the polezero diagram, and the frequency response. The output of the signal analysis routines is meant to be selfexplanatory at the level of a student studying linear systems theory. For example, Figure 5.4 shows the general analysis of the impulse response of an analog lter. Also, Figure 5.5 shows general analysis of the upsampling by 3 of a discretetime i.e. sampled sinc function.
5.4 Analysis of Multirate Systems Using Transforms
One way to understand the behavior of a system is to derive its inputoutput relationship in a transform domain. The inputoutput relationship is simply a transfer function for linear shiftinvariant LSI systems, but it becomes a transfer function plus aliasing terms for linear periodically timevarying LPTV systems 49 . When deriving inputoutput relationships on paper, engineers implicitly use transform pairs such as x n ! X z to represent the ztransform of the input signal x n . For LSI and LPTV systems, the linear transform routines can derive inputoutput relationships if all of the unknown quantities input signals, unspeci ed lters, and so forth are assigned transform pairs. This ability applies to onedimensional and multidimensional, as well as continuoustime and discretetime, systems. Discretetime LPTV systems are also known as multirate systems. Figure 5.6 shows the derivation of the inputoutput relationship for a twochannel nonuniform onedimensional lter bank. The block diagram of the lter bank is shown in Figure 5.6a. In the lter bank, the upper channel contains the frequency band , ; of the input signal, whereas the lower channel contains the frequency band ,; , ; . Both channels are resampled at their Nyquist rates. Figure 5.6b gives the de nition of the lter bank structure using the notation of the signal
2 3 2 3 2 3 2 3
79
Analysis of Multirate Systems Using Transforms
In[2]:=
ASPAnalyze[ t Exp[a t] Cos[3 Pi t / 16] CStep[t], t, 0, 3 Pi, a > 1 ] For plotting only, these symbols will will be assigned default values: {a > 1} Real part of the function is shown as solid lines. Imaginary part of the function is shown as dashed lines.
ContinuousTime Domain Analysis 0.3 0.25 0.2 0.15 0.1 0.05 2 0.05 4 6 8 t
Since the signal is stable, the frequency response can be computed directly from the Laplace transform. 3 Pi t t Cos[] CStep[t] 16 a t E has the following frequency response: 2 2 2 256 (256 a  9 Pi + 512 I a w  256 w ) 2 2 2 2 (256 a + 9 Pi + 512 I a w  256 w )
Magnitude Response (dB) 0.01 0.05 0.1 0.5 10 15 20 25 30 35 40 5 10. w
3 Pi t t Cos[] CStep[t] 16 a t E has the following Laplace transform: 2 2 2 256 (256 a  9 Pi + 512 a s + 256 s ) 2 2 2 2 (256 a + 9 Pi + 512 a s + 256 s ) The region of convergence is: Re[a] < Re(s) < Infinity The system is stable if Re[a] < 0 The default values are now being considered.
Out[2]= 0.01
Phase Response (degrees) 0.05 0.1 0.5 25 50 75 100 125 150 5 10. w
The zeroes are: The poles are:
{1.58905, 0.410951} {1. + 0.589049 I,
LTransData[ 2 2 2 256 (256 a  9 Pi + 512 a s + 256 s ) , 2 2 2 2 (256 a + 9 Pi + 512 a s + 256 s ) Rminus[Re[a]], Rplus[Infinity], LVariables[s]]
1.  0.589049 I, 1. + 0.589049 I, 1.  0.589049 I}
Im s 0.6 0.4 0.2 1.75 O 1.5 1.25 0.2 0.4 0.6 0.75 0.5
*X*
O
Re s
*X*
The signal analyzer plots the impulse response of the of the analog lter over t 2 0; 3 when the parameter a takes value of 1. The analyzer nds the Laplace transform and displays the polezero diagram with the ROC shaded. It also prints the formula describing the frequency response and plots the magnitude and phase of the frequency response.
Figure 5.4: OneDimensional Analog Signal Analysis
80
Analysis of Multirate Systems Using Transforms
In[2]:=
DSPAnalyze[ Upsample[3,n][ Sinc[n] ], n, 20, 20 ]
DiscreteTime Domain Analysis 1 0.8 0.6 0.4 0.2 20 15 10 5 0.2 5 10 15 20 n
Transform::incomplete: The rule base could not compute the forward ztransform of Sinc[n] with respect to n. Upsample 3,n has the following frequency response: Pi Periodic 2 Pi 3 in w [CPulse 1 [ + w]] 2/3 3 [Sinc[n]]
Magnitude Response 3 2.5 2 1.5 1 0.5 6 Out[2]= 4 2 2 4 6 w
Incomplete zTransform
Figure 5.5: OneDimensional DiscreteTime Analysis of a Multirate Signal The ztransform does not exist for this case, so only the timedomain signal and its Fourier transform are displayed.
81
Analysis of Multirate Systems Using Transforms

"n 2

h n

0
n 3

"n 3

g n

0
n 2
xn
+
?
6
xn ^


h n

1
n 3

"n 3

g n
1
a Flow graph of the lter bank
upperchannel = Downsample 2,n Convolve n g0 n , Upsample 3,n Downsample 3,n Convolve n h0 n , Upsample 2,n x n lowerchannel = Convolve n g1 n , Upsample 3,n Downsample 3,n Convolve n h1 n , x n
b Representation of the lter bank in the new environment
Xhat z = ZTransform lowerchannel + upperchannel, n, z, TransformLookup x n : X z , g0 n : G0 z , g1 n : G1 z , h0 n : H0 z , h1 n : H1 z
c Mathematica code to generate the inputoutput relationship
TeXForm
1 6
Collect
p p p p ^ X z =
G0 , zH0 , z + G0 zH0 z + 2G1 zH1 z X z + p p p p G0 , zH0 e 3 z + G0 z H0 e 43 z + 2G1 z H1 e 23 z X e 23 z +
p p p p G0 zH0 e 23 z + G0 , z H0 e 53 z + 2G1 z H1 e 43 z X e 43 z
1 6 1 6
i i i i i i i i
SPSimplify 6 TheFunction Xhat z , X z , X Exp 2 I Pi 3 z , X Exp 4 I Pi 3
z
d Inputoutput relationship generated by the environment
Figure 5.6: Deriving InputOutput Relationship for a NonUniform Filter Bank
82
Multidimensional Stability Analysis
processing packages. Figure 5.6c uses the forward ztransform routine to generate ^ the relationship between X z and X z in the z domain. Figure 5.6d converts the algebraic inputoutput relationship to its TEX form using the builtin abilities of Mathematica's TeXForm command. In order to make the formula look better in TEX, ^ we wrote the Mathematica command to multiply X z by six, then collect terms, and nally generate the equivalent TEX code. By hand, we added the factors to each term.
1 6
5.5 Multidimensional Stability Analysis
Multidimensional stability analysis is di cult to describe analytically. For some classes of multidimensional systems, the conditions for stability can be derived by examining the region of convergence ROC of their generalized frequency transforms Section 5.5.1. Once values have been assigned to free parameters, stability becomes much easier to analyze by visualization Section 5.5.2.
5.5.1 Symbolic Analysis of Stability
Knowing the ROC for a signal leads to a set of stability conditions involving the free parameters of the signal. A discretetime signal is stable if the ROC of its ztransform R, jzj R contains the unit circle, i.e. R, 1 R . A continuoustime signal is stable if the ROC of its Laplace transform R, es R contains the imaginary axis, i.e. R, 0 R . The Stable command uses these simple comparisons to convert ROC information to a set of stability conditions. These comparisons generalize to multiple dimensions. Figure 5.7 shows an example of the stability analysis of a twodimensional nonseparable signal taken from 53 . In multiple dimensions, nonseparable signals have nonseparable regions of convergence, which implies that one or more of the ztransform variables will appear in the ROC. When checking for stability, the ztransform variables take values on the
+ + + +
83
Multidimensional Stability Analysis
In 36 := ZTransform Out 36 = 1 ZTransData , Rminus a b 1    z1 z2 Rplus In 37 := Stable Out 37 = Abs a Abs a , Abs b 1 Abs a 1  z1 z1, z2 ,
n1 + n2! a^n1 b^n2 Step n1,n2 n1, n2 , z1, z2
n1! n2!,
Infinity, Infinity
, ZVariables
1 && Abs b
1  Abs a
In 38 := Assuming All Out 38 = Abs a 1 && Abs b 1 Abs a 1  z1 a 1 && 1   != 0 z1
The above interaction demonstrates the stability analysis of a twodimensional nonseparable sequence in the variables n1 and n2 having the form n + n ! an1 bn2 u n ; n n!n! Since the function is not separable, neither is the region of convergence ROC. Stable generates the inequality for checking if the unit bisphere jzij = 18i is in the ROC. If the inequality does not evaluate to true or false, Stable will apply a set of simpli cation rules to it with the side information jzij = 18i. The actual stability condition is jaj + jbj 1 53 which the stability checker returns as jbj 1 ,jaj note that the rst condition jaj 1 is actually redundant. By evaluating Assuming All , the MDSPPs will list the assumptions it has made about free parameters In the output of Assuming All above, the rst two inequalities are the original inequality conditions generated from the ROC. The third inequality arose during the simpli cation of the second inequality.
1 2 1 2 1 2
Figure 5.7: Stability Analysis of a NonSeparable TwoDimensional Signal
84
Multidimensional Stability Analysis
multidimensional unit circle: jz j = 1, jz j = 1, etc. In N dimensions with ztransform variables fzigN , the comparison test for stability 53 generalizes to i
1 2 =1
i
=1::: N; i=j
jz j=1
i
max6 R, 1
i
=1::: N;i=j
jz j=1
i
min6 R
+
where R, jzj j R
+
for j = 1 : : : N
Similarly for the Laplace transform of continuoustime signals, the stability test generalizes to
i
=1::: N; i=j
es =0
i
max6 R, 0
i
=1::: N;i=j
es =0
i
min6 R
+
where R,
esj R
+
for j = 1 : : : N
The stability checking routine replaces any esi terms in the Laplace transform ROC with zero.
5.5.2 Graphical Analysis of Stability
The results of transforms are sometimes best understood in terms of graphical representations. A wide array of graphical capabilities have been implemented in the MDSPPs for 1D and 2D signals in both the discretetime and continuoustime domains: continuous and discrete plots, magnitude and phase plots, and polezero diagrams. For nonseparable 2D systems, a polezero diagram can be displayed as the locus of poles and zeros of one transform variable parameterized by the other 53 . By inspection of the loci of poles and zeros, one can determine the stability of the signal by using the same rules as in one dimension: if the poles lie outside of the unit circle, then the causal implementation of the system would be unstable. Underlying the polezero root loci diagrams is a root locus plotting command for one free variable. Figure 5.8 analyzes the location of the poles and zeros for a signal. Zeros are shown as O's, poles are shown as X's, the unit circle is shown as a solid curve, and the lower bound of the ROC R, is graphed with dashed lines. The upper bound of the ROC R is not shown because it is in nite for both transform variables. The left column of Figure 5.8 shows the zero and pole diagrams of the second transform variable z2 projected onto z1 by the mapping z = expj! . The right column
+ 2 2
85
Multidimensional Stability Analysis
PoleZeroPlot ZTransform n1+n2! 4 5^n1 2 7^n2 Step n1,n2 n1, n2 , z1, z2
I w Numerator polynomial in z1: 35 E Denominator polynomial in z1: I w 28 E  10 z1 + 35 E
ZeroIm z1 Map Root 2 1.5 1 0.5
n1! n2!,
I w
z1
Numerator polynomial in z2: 35 E Denominator polynomial in z2: I w I w  28 z2 + 35 E
ZeroIm z2 Map Root 2 1.5 1 0.5
z2
I w z1 10 E
z2
2
1.5
1
0.5
O
0.5
1
1.5
2
Re z1
2
1.5
1
0.5
O
0.5
1
1.5
2
Re z2
0.5 1 1.5 2
0.5 1 1.5 2
PoleIm z1 Map Root 2 1.5 1 0.5
PoleIm z2 Map Root 2 1.5 1 0.5
2
1.5
1
0.5 0.5 1 1.5 2
0.5
1
X
1.5
2
Re z1
2
1.5
1
0.5 0.5 1 1.5 2
0.5
1
X 1.5
2
Re z2
Figure 5.8: TwoDimensional PoleZero Diagrams for an Unstable Signal The lower limit of the region of convergence ROC is plotted as dashed lines. The ROC extends to in nity in both z variables. In parts of the root loci, the poles and zeros actually move outside of the ROC because the system is unstable.
86
Analysis of Multidimensional Multirate Systems Using Transforms
of Figure 5.8 shows the zero and pole diagrams of the rst transform variable z1 projected onto z2 by the mapping z = expj! . Since the poles fall outside of the unit circle, the signal is unstable.
1 1
5.6 Analysis of Multidimensional Multirate Systems Using Transforms
This section discusses the use of the MDSPPs to analyze multidimensional multirate systems. First, Section 5.6.1 applies general twodimensional signal analysis to a downsampled signal. Second, Section 5.6.2 shows how the MDSPPs can help in the visualization of aliasing in two dimensions. Last, Section 5.6.3 performs a purely symbolic analysis of resampling.
5.6.1 Automated TwoDimensional Signal Analysis
We now elaborate on the twodimensional signal analysis abilities introduced in previous sections of this chapter. Mathematica can already plot twodimensional continuoustime functions as threedimensional mesh plots via its Plot3D command. Plot3D can also be used to generate the mesh plots for twodimensional sequences. For twodimensional sequences, however, a better alternative is to use DensityPlot to render an aerial view of the domain and range see Figure 4.1. In two dimensions, polezero diagrams for ztransforms become root maps 53 , as mentioned in Section 5.5.2. Likewise for twodimensional Laplace transforms, the polezero diagram is a root map s is plotted with s = j! for di erent values of ! and visaversa for s . Frequency responses are plotted using a parametric plotting routine. Figure 5.9 performs analysis on a discretetime signal which su ers from aliasing distortion.
1 2 1 1 2
87
Analysis of Multidimensional Multirate Systems Using Transforms
lowpass2Dtile = CPulse[Pi, w1 + Pi/2] CPulse[Pi, w2 + Pi/2]; lowpass2D = InvDTFTransform[ lowpass2Dtile, {w1, w2} ] n1 Pi n2 Pi Sinc[] Sinc[] 2 2 4 downMatrix = {{2, 1}, {1, 3}}; DSPAnalyze[ Downsample[downMatrix, {n1,n2}][lowpass2D], {n1, n2}, {5, 5}, {5, 5} ]
DiscreteTime Domain Analysis
CPulse
Pi 3 (2 Pi + w1) 4 Pi + w2 [ +   ] Pi 2 5 5 Pi 2 Pi + w1 2 (4 Pi + w2) [   + ] + Pi 2 5 5
CPulse
CPulse
Pi 3 (2 Pi + w1) 2 Pi + w2 [ +   ] Pi 2 5 5 Pi 2 Pi + w1 2 (2 Pi + w2) [   + ]) / 5] Pi 2 5 5
Magnitude Response
CPulse
0.2 2 0.1 .1 0 4 2 0 n1 2 4 w1 4 2 0 2 0.5 .5 n2 0 0 2 4 6 0 2 4 w2 4 1 6
Transform::incomplete: The rule base could not compute the forward n1 ztransform of Sinc[] with respect to n1. 2 ZTransform::notvalid: The forward ztransform could not be found. n1 Pi n2 Pi Sinc[] Sinc[] 2 2 [] n1 4 n2 Pi 3 w1 w2 [ +   ] Pi 2 5 5
0.5 5 0.25 5 0 0.25 25 0.5 .5 0
Phase Response (degrees)
Downsample 2 1 1 3  
6 4 w2 2 w1 4 6 0 2
has the following frequency response: Periodic 2 Pi 0 CPulse 0 2 Pi   w1 w2 [(CPulse
Incomplete zTransform
Pi w1 2 w2 [   + ] + Pi 2 5 5
CPulse
Pi 3 (4 Pi + w1) 6 Pi + w2 [ +   ] Pi 2 5 5 Pi 4 Pi + w1 2 (6 Pi + w2) [   + ] + Pi 2 5 5
CPulse
CPulse
Pi 3 (4 Pi + w1) 4 Pi + w2 [ +   ] Pi 2 5 5 Pi 4 Pi + w1 2 (4 Pi + w2) [   + ] + Pi 2 5 5
CPulse
The rst command de nes a 2D timedomain signal that is an ideal quarterband lowpass lter. The second command invokes the digital signal analyzer. The analyzer rst plots the time response, then reports that it cannot nd the z transform because of the twosided sinc functions, and nally computes and plots the DTFT. Downsampling by 5 quintuples the area of support in the the frequency domain which causes aliasing. Aliasing appears in the magnitude plot for amplitudes above 0.2.
Figure 5.9: Signal Analysis of Aliasing in Two Dimensions
88
Analysis of Multidimensional Multirate Systems Using Transforms
5.6.2 Visualization of Downsampling in Two Dimensions
We will now use the MDSPPs to visualize and analyze the e ects of downsampling by matrix M on a baseband signal X !. The baseband signal X ! is periodic with period of 2 in each discretetime frequency variable. In terms of the discretetime frequency variables !, downsampling by M periodically replicates the baseband j det M j times in each 2 2 frequency tile. The centers of the new frequency bands occur at the original center of the baseband shifted by the aliasing vectors associated with M see Section 3.1. The aliasing vectors, computed according to
, equation 3.5, are: !i = 2 M T ki, where i = 0; : : : ; j det M j , 1 and ki is the ith coset vector. In the MDSPPs, the DownsamplingAnalysis routine determines if the downsampled signal experiences aliasing and if it covers the frequency domain. A downsampled signal can su er from two kinds of aliasing: intraband aliasing and interband aliasing. In intraband aliasing, a frequency band overlaps with itself, whereas in interband aliasing, a frequency band overlaps with another frequency band. Clearly, if a downsampled signal does not experience aliasing and covers the frequency domain, then the signal has been resampled at its Nyquist rate. DownsamplingAnalysis represents the downsampling operation by its downsampling matrix and the baseband signals as a list of polygons. For example, Figure 5.10 shows the e ect of downsampling by a quincunx downsampling matrix on a diamondshaped passband. A quincunx resampling matrix has a determinant of 2 or ,2 and nonzero entries. The lattice sampling grid generated by quincunx matrices correspond to diagonal interleaving, which is commonly found in video signals. The vertices of the diamondshaped passband are obtained by mapping a separable lowpass lter by a quincunx matrix. In Figure 5.10, the quincunx downsampler resamples the baseband signal at its Nyquist rate. In Figure 5.11, we keep the same baseband but use a di erent downsampling matrix. The downsampling matrix has a determinant of three. Since the baseband
1
89
Analysis of Multidimensional Multirate Systems Using Transforms
DownsamplingAliasing[ {{1, 1}, {1, 1}}, Polygon[ {{Pi,0}, {0,Pi}, {Pi,0}, {0,Pi}} ], Dialogue > All ] Analyzing aliasing for the downsampling matrix 1 1 1 for a baseband whose domain is 1 described by the polygon with vertices {{Pi, 0}, {0, Pi}, {Pi, 0}, {0, Pi}} The downsampling will yield the baseband plus 1 shifted/skewed copy of the baseband. The shifting vectors are {{0, 0}, {Pi, Pi}} The area of the baseband is 2 2 Pi whose numerical value is 19.7392. Initial analysis: the downsampler resamples the baseband signal at its Nyquist rate (assuming that no aliasing is present). Band #1 of 2 is free of intraband aliasing.
Band #2 of 2 is free of intraband aliasing.
π
Band #2 of 2
0
?π
0
π
Tiling of input frequency domain for the 2 bands:
π
Tiling of Input Frequency Domain
π
Band #1 of 2
0 0
?π ?π
0
π
The downsampler does not introduce aliasing.
0
π
It resamples the baseband at the Nyquist rate.
Figure 5.10: Quincunx Downsampling Without Aliasing
90
Analysis of Multidimensional Multirate Systems Using Transforms
occupies 50 of the frequency domain, generating three copies of the baseband will produce interband aliasing. DownsamplingAliasing reports that no intraband aliasing occurs and that interband aliasing occurs in 50 of the frequency domain. The DownsamplingAliasing routine works as follows. For each new band, the routine computes the location of the vertices of the band, decomposes the band into a union of triangles, and determines how each triangle maps into the fundamental frequency tile ! 2 ,; S ! 2 ,; by calling the TriangleModWithSquare routine. The TriangleModWithSquare routine is recursive. For each triangle,
1 2
if each of the triangle vertices is either inside or on the border of the fundamental tile, then return the triangle; if one or more of the corners of the fundamental tile is inside the triangle, then split the triangle into several triangles by using the corners of the fundamental tile inside of it as the new vertices and pass each new triangle to TriangleModWithSquare; if the triangle is completely outside of the fundamental tile, then shift the triangle by the right number of periods to put at least one vertex in the fundamental tile and pass it to TriangleModWithSquare; and if none of the above cases apply, then the triangle intersects one or more of the edges of the fundamental tile, so we split the triangle at each edge intersection into smaller triangles: return the triangles inside the fundamental tile, and pass triangles outside of the fundamental tile to TriangleModWithSquare. Once the new band has been generated, the DownsamplingAliasing routine checks for intraband aliasing by checking for overlapping regions in the band via the OverlappingRegions routine. The new band is then plotted with intraband aliasing region if it exists shaded black. After all the bands have been generated, the
91
Analysis of Multidimensional Multirate Systems Using Transforms
DownsamplingAliasing[ {{1, 1}, {2, 1}}, Polygon[ {{Pi, 0}, {0, Pi}, {Pi,0}, {0,Pi}} ], Dialogue > True ] Analyzing the aliasing for the downsampling matrix 1 2 1 for a baseband whose domain 1 is described by the polygon with vertices {{Pi, 0}, {0, Pi}, {Pi, 0}, {0, Pi}} Band #1 of 3 is free of intraband aliasing. Band #3 of 3 is free of intraband aliasing.
π
Band #1 of 3
π
Band #3 of 3
0
0
?π
0
π
?π
0
π
Tiling of the input frequency domain for the 3 Band #2 of 3 is free of intraband aliasing. bands with interband aliasing shown in black:
π
Band #2 of 3
π
Tiling of Input Frequency Domain
0
0
?π
0
π
?π
0
π
The downsampler introduces interband aliasing. Aliasing occurs in 50% of the frequency domain.
Figure 5.11: Quincunx Downsampling With Aliasing
92
Analysis of Multidimensional Multirate Systems Using Transforms
routine checks for interband aliasing by checking for overlapping regions between bands. The routine then computes the total area of the domain in which aliasing occurs.
DownsamplingAliasing
5.6.3 Automatic Derivation of Transform Properties
When all the free parameters have been assigned values, graphicalbased signal analysis provides much insight. Sometimes, symbolic analysis can help in the choice of values for free parameters. For example, Section 5.4 showed that the transform routines can derive the inputoutput relationship of a system in a transform domain if abstract" transform pairs such as x n ! X z are provided. In this section, we examine the use of abstract transform pairs to derive properties of the transforms of multidimensional multirate signals. When3 upsampling a signal in two dimensions say by some upsampling matrix 2 u u 5 4 , the resulting transform is u u
11 21 12 22
X zu11 zu21 ; zu12 zu22
1 2 1 2
given that the ztransform of the original signal is X z ; z . The MDSPPs can generate this relationship by specifying the abstract transform pair x n ; n ! X z ; z :
1 2 1 2 1 2
In 2 := upMatrix = In 3 := ZTransform
u11, u12 ,
u21, u22
; x n1,n2 : X z1,z2 ,
Upsample upMatrix, n1,n2 , z1, z2 , TransformLookup 
n1, n2 x n1,n2
u11 u21 u12 u22 Out 3 = ZTransData X z1 z2 , z1 z2 , Rminus Rplus Infinity, Infinity , ZVariables z1, z2
0, 0
,
93
Summary
The MDSPPs will also track the region of convergence of the upsampled signal if we give the ROC of the input signal:
In 4 := ZTransform Upsample upMatrix, n1, n2 n1, n2 , z1, z2 , TransformLookup x n1,n2 : X z1,z2 , x n1, n2 ,
rm1,rm2 ,
rp1,rp2
u11 u21 u12 u22 Out 4 = ZTransData X z1 z2 , z1 z2 , 1 u11 1 u21 1 u12 1 u22 rm1 rm2 , rm1 rm2 1 u11 1 u21 1 u12 1 u22 rp1 rp2 , rp1 rp2 z1, z2
Rminus
,
Rplus
,
ZVariables
The resulting multidimensional ROC is not separable. For each zero element uij in the upsampling matrix, the corresponding xu and x =u terms equal 1 for any value of x. When the input signal is stable, the ROC of the ztransform of the upsampled signal will be a compressed version of the ROC of the ztransform of the input signal because of the nth root terms. The MDSPPs can track the ROC in higher dimensions.
ij
1
ij
5.7 Summary
This chapter rst identi es important properties of multidimensional signals based on the properties of onedimensional signals used in SPLICE and ADE and then discusses the highlevel capabilities of the signal processing packages. These highlevel routines, which are summarized in Table 5.2, rely on the representation of signals and systems as functions and operators, respectively, as discussed in Chapter 4. The highlevel routines cover a wide range of graphical and symbolic analyses for both
94
Summary
New Function
ASPAnalyze CTFTransform DFTransform DTFTransform DSPAnalyze InvCTFTransform InvDFTransform InvDTFTransform InvLaPlace InvZTransform LaPlace LSolve PiecewiseConvolution PiecewisePlot PoleZeroPlot RootLocus SequencePlot SignalPlot Stable ZSolve ZTransform
Dimension of Signal 1D and 2D 1D and mD 1D and mD 1D and mD 1D and 2D 1D and mD 1D and mD 1D and mD 1D and mD 1D and mD 1D and mD 1D 1D 1D 1D and 2D 1D 1D and 2D 1D and 2D 1D and mD 1D 1D and mD
Description general analog signal analyzer forward continuoustime Fourier transform forward discrete Fourier transform forward discretetime Fourier transform general digital signal analyzer inverse continuoustime Fourier transform inverse discrete Fourier transform inverse discretetime Fourier transform inverse Laplace transform inverse z transform forward Laplace transform linear constant coe cient di erential equations solver discrete and continuous convolution piecewise function plotter display of polezero diagrams root locus of one free parameter digital signal plotter analog signal plotter stability checker linear constant coe cient di erence equations solver forward z transform
Table 5.2: HighLevel Abilities of the Signal Processing Packages
95
Summary
discretetime and continuoustime signals and systems. The highlevel capabilities work on nonseparable multidimensional signals and systems, except for the convolution, root locus, and di erential di erence equation solving routines which only work for separable ones. The extended graphics routines help engineers to visualize onedimensional and twodimensional signals and systems. The transform, convolution, and di erential di erence equation DE solvers can provide closedform symbolic answers to many problems that would otherwise be di cult and prone to error when worked by hand. In addition, these routines can show how they compute the answer in a stepbystep manner for veri cation. This dialogue is in natural language for the transform and DE solvers, but it takes the form of animation for convolution. We have combined graphical and symbolic analyses for onedimensional and twodimensional signals into two general signal analyzers, one for continuoustime and one for discretetime signals.
96
CHAPTER 6 Rearranging Multidimensional Systems
This chapter discusses the rearrangement of algorithms with the goal of nding better or even optimal implementations. By rearranging components of an algorithm, ESPLICE and ADE have derived e cient implementations for three speci c onedimensional multirate systems 31 . ESPLICE derived the fact that the cascade of an upsampler and downsampler commutes if their resampling factors are relatively prime 9 . ESPLICE and ADE have generated polyphase forms for systems to change the sampling rate by rational factors 9, 11 . ADE has found an e cient multiband structure for the discretetime implementation of optimal detectors frequencychirp matched lters of FSKcoded sonar signal beams based on the pruned FFT 10, 12 . ESPLICE and ADE know how to rearrange cascade and parallel combinations of linear onedimensional multirate operators. Many of the rules are based on operator properties such as linearity. Section 6.1 extends these properties to multidimensional multirate operators and identi es new properties that became important in multiple dimensions. ESPLICE and ADE do supplement their propertybased rules to describe how to rewrite combinations of operators, especially between resamplers and other operators. We have identi ed the multidimensional versions of these rules 39, 52, 53, 71, 72 and encoded them in the multidimensional signal processing packages MDSPPs, as discussed in Chapter 3. Section 6.2 describes di erent strategies to apply a set of rearrangement rules. ESPLICE and ADE utilize cost functions to rank the e ciency of equivalent implementations of an algorithm. ESPLICE measures the number of additions and the number of multiplications, whereas ADE also takes the amount of memory needed into account. Given a cost function, rearrangement rules can be applied until a
97
Extending System Properties in ESPLICE and ADE
better or even optimal implementation is discovered. ADE uses heuristics to reduce the number of equivalent implementations called the equivalence space that it generates. The MDSPPs use heuristics as well, as discussed in Section 6.3.
6.1 Extending System Properties in ESPLICE and ADE
In order to prevent a combinatoric explosion in the number of rearrangement rules, ESPLICE and ADE base many of their rearrangement rules on system properties such as those given earlier in Table 2.2 c.f. Appendix D of 8 . These system properties extend naturally to higher dimensions. The MDSPPs add the properties DISCRETE and CONTINUOUS to distinguish the domain of the input signals. This was not an issue in ESPLICE and ADE because they only supported discretetime signals. We have also introduced the properties LINEARPHASE and DELAY because they are important parameters in image processing and communications, respectively. Another new property, SEPARABLE, indicates if a multidimensional operation can be decomposed into separable onedimensional operations. Table 6.1 lists the properties that the MDSPPs implement. To complement their propertybased rules, ESPLICE and ADE introduce rules to capture special relationships between operators. The rules containing operators that are separable in multiple dimensions can be updated easily. Most of the rules, however, do not fall into this category because they contain downsampling, upsampling, or ltering operations which are not always separable in multiple dimensions. Chapter 3 discusses in detail the rules for rewriting cascades of multidimensional multirate operations, which have been encoded in the MDSPPs.
98
Extending System Properties in ESPLICE and ADE
System Property Meaning ASSOCIATIVEy can change grouping of inputs ADDITIVEy distributes over addition y COMMUTATIVE can change order of inputs CONTINUOUS inputs are continuous signals DELAY amount of delay before output is meaningful DISCRETE inputs are discrete signals y HOMOGENEOUS scaled input gives scaled output LINEARy additive and homogeneous LINEARPHASE true if the frequency phase response is a linear function of the frequency variable y MEMORYLESS output does not depend on previous inputs or outputs; if a singleinput system, then
SHIFTINVARIANT SEPARABLE SHIFTINVARIANTy
true if separable in all dimensions, false if completely nonseparable, or a list of variables in which the operator is separable shifted input gives shifted output
y
property also included in ESPLICE 8
Table 6.1: System Properties in the MDSPPs
99
Algorithm Rearrangement
6.2 Algorithm Rearrangement
In rearranging operations in an algorithm, both simpli cation and rearrangement rules are applied to the algorithm. Simpli cation rules produce more e cient algorithms through the removal of one or more components e.g. a shift by l undoes a shift by ,l. Rearrangement rules produce di erent but not necessarily better implementations. For linear systems, some rearrangement rules can be based on system properties e.g. the order of two linear shiftinvariant operators can be switched, but most rules are expressed in terms of speci c operators e.g. the cascade of two shift operators can be rewritten as one shift operator whose shift is the sum of the two original shift factors.
6.2.1 General Procedure
The procedure to rewrite algorithms is iterative. At each iteration, the rst step is to determine which rearrangement rules apply, beginning at the outputs of the algorithm. The search for applicable rules may continue until the inputs to the algorithm is reached. Then, the algorithm is rewritten by applying one or more of the rules. Finally, at each iteration, each equivalent form of the algorithm is simpli ed. The end result is a tree of equivalent forms. In determining which rules may apply to an algorithm, rules that would produce an expression that has already been generated would not be considered. If maintaining parallelism in a part of the algorithm is important, then rules would have to apply to every branch in a parallel structure called a regularity constraint in ADE. Once a collection of candidate rules has been determined, the procedure could apply all of them to the algorithm which would likely create a tree that would grow exponentially in size. The other extreme is to apply only one rule at each iteration e.g. choose one at random. These two opposing strategies correspond to a breadth rst search and a depth rst search, respectively, through the equivalence space the tree of equivalent
100
Algorithm Rearrangement
implementations. Clearly, the second approach would be suboptimal, but the rst approach may not nd the optimal solution in a reasonable amount of time.
6.2.2 Rearranging A Multidimensional Rational Rate Changer
A simple example of algorithm rearrangement is a structure to change the two2 3 6 85 dimensional sampling rate rst by upsampling the input signal x n ; n by 4 3 3 2 3 9 12 5 and then by downsampling the result by 4 . In the MDSPPs, the structure 7 7 would be written as a nested algebraic expression:
1 2
Downsample Upsample
9, 12 ,
7, 7
,
n1,n2 , n1,n2 x n1,n2
6, 8 ,
3, 3
The rearrangement rule shown in Figure 3.12a will seek to decompose up downsampling cascades by rst factoring the up downsampling matrices into their Smith forms and then by removing redundant resampling operations. In this case, 2 3 2 32 3 2 3 2 32 3 6 85 42 0543 45 9 12 5 4 3 0 5 4 3 4 5 4 4 = = 3 3 0 3 1 1 7 7 0 7 1 1 so the common right unimodular matrix can be removed, thereby making the resampling separable:
Downsample Upsample 3, 0 , 0, 7 , n1,n2 , n1,n2 x n1,n2 2, 0 , 0, 3
Another rearrangement rule decomposes separable multidimensional resampling operations into onedimensional resampling operations:
Downsample 3, n1 Downsample 7, n2 Upsample 3, n2 x n1,n2 Upsample 2, n1
Applying the rearrangement rule that commutes downsamplers in cascade and another that commutes upsamplers in cascade, we ultimately get
101
Finding Optimal Algorithms
Upsample 2, n1 Downsample 3, n1 Downsample 7, n2 x n1, n2 Upsample 3, n2
The overall rate change is the product of the original downsampling matrix and the inverse of the upsampling matrix:
2 9 4
12 5 4 6 8 5, 4 0 5 = 0 7 7 3 3
1 2 7 3
32
3
23
3
6.3 Finding Optimal Algorithms
Adding the ability to measure the performance of equivalent forms of algorithms can help guide the rearrangement process to search for better implementations. At each iteration in the depth rst search strategy for applying rearrangement rules, the lone rule could be chosen as the one among the candidates that produces the best new implementation. When no more rules apply, the best implementation would be chosen from all of the equivalent forms generated. For the breadth rst search strategy, candidate rules could be pruned by only selecting those that will produce an algorithm with a lower cost or increase the cost over the current or original algorithm by no more than some factor. Both ESPLICE and ADE implement cost functions to rank equivalent algorithms 77 . ESPLICE counts only the number of additions and multiplications times the sampling rate. ADE also counts the number of memory elements which is independent of the sampling rate. The MDSPPs measure all three costs. In multiple dimensions, however, the sampling rate is no longer a real number but instead a realvalued matrix. The MDSPPs measure the number of additions and multiplications times the average output sampling rate. Unlike ESPLICE and ADE, the MDSPPs can associate a range of values with a free parameter. The implementation cost of linear operators for each free parameter either monotonically increases or monotonically decreases for each free parameter.
102
Summary
When one free parameter is allowed to take values over an interval, the associated implementation cost for an operator depending on that free parameter also takes values over an interval. Even when an operator depends on two or more such free parameters, the cost function still takes values over an interval. As a consequence, the implementation cost of two algorithms can be compared in terms of intervals.
6.4 Summary
The primary goal behind rearranging the form of an algorithm is to nd a better implementation of algorithm. To nd better implementations, a computer program needs a set of simpli cation and rearrangement rules as well as a function to measure the implementation cost of an algorithm. When applied, simpli cation rules will produce a better implementation because they remove inverse operations that appear in cascade. Rearrangement rules, on the other hand, simply replace a part of an algorithm by an equivalent form that may or may not be more e cient. ESPLICE and ADE express some rearrangement rules for onedimensional systems in terms of system properties, but the majority are based on speci c operators. Since system properties such as linearity extend to higher dimensions, we have recast the rearrangement rules based on system properties for multidimensional systems. We also add new properties and base new rules on them. In Chapter 3, we derive the special case rules in multiple dimensions. We also extend the simpli cation rules present in ESPLICE and ADE to multiple dimensions. In computing computational cost for onedimensional multirate algorithms, the number of additions and multiplications is scaled by the sampling rate. For multiple dimensions, the computational cost depends instead on the sampling matrix instead of the sampling rate. The other costs such as the amount of memory elements are independent of the sampling rate and extend naturally to higher dimensions.
103
CHAPTER 7 Generating Equivalent Code for Algorithms
As the previous two chapters show, symbolic reasoning can assist the engineer in designing algorithms with a relatively small number of free parameters and interconnections see Figure 5.6, for example. Because we implement reasoning on mathematical formulas, we do not have a convenient way to represent the many connections in complex systems. The layout of complex systems is better captured graphically by tools that can organize block diagrams hierarchically. Generating equivalent code for developed algorithms is essential to transport the algorithm to another environment for layout, typesetting, further simulation, and so forth. Speci cally, the MDSPPs translate signal processing expressions into TEX document processing commands 78 and Ptolemy block diagram representations 13, 14, 15, 16 . We have chosen TEX because it is a typesetting tool commonly used by engineers to document designs, write reports, etc. We have chosen the Ptolemy environment because it can represent complex systems and simulate algorithms as well as generate C programs, DSP assembly language code, and VHDL descriptions. The ability to generate Ptolemy code places the MDSPPs in the rapid prototyping process shown earlier in Table 1.1.
7.1 Generating TEX Code
TEX, a typesetting language developed by Donald Knuth, formats text and equations in a standard way for publication. When formatting equations, TEX displays mathe
104
Generating TEX Code
matical formulas according to the conventions de ned by the American Mathematical Society. When typesetting text, TEX follows standard rules followed by printers. We have chosen TEX as a target environment because it is commonly used by engineers to document projects, write technical reports, etc. An inherent ability of Mathematica is to translate a single algebraic expression into a TEX equation. These equations can then be spliced into existing documents. Ideally, this automatic translation eliminates the possibility of errors that can occur while retyping a formula. For example, TeXForm t^2 gives t^2 which is interpreted as t once the expression is put into a TEX math environment i.e. surrounded by dollar signs. A more complicated example is the inputoutput relationship of a twochannel nonuniform lter bank shown in Figure 5.6c. We rst had Mathematica factor the expression in a meaningful way and then called TeXForm to generate the equation. In the extended Mathematica environment, the multidimensional signal processing packages contain the TEX de nitions for signals and systems they introduce. For example,
2
TeXForm
Downsample 2,n
x n
generates one line of TEX code downarrow_ 2,n x n which, after TEX processing, gives ;n x n . For a more complicated example, we will use the time domain description of the twochannel nonuniform lter bank shown in Figure 5.6b:
2
TeXForm
upperchannel + lowerchannel
gives
g1 n?n "3;n 3;n h1 n ?n xn+ 2;n g0 n?n "3;n 3;n h0 n?n "2;n xn
Note that ?n means convolution in the variable n. The impulse response of the lters are not subscripted here as we did not encode a rule to print g0 as g .
0
105
Generating Complete Ptolemy Simulations
7.2 Generating Complete Ptolemy Simulations
We now consider issues involved in converting a set of mathematical formulas into a complete working program. In our case, we are ultimately concerned with converting a set of mathematical formulas into its Ptolemy block diagram speci cation with the option of generating additional code to direct the running of a complete simulation, as is discussed in Section 7.2.2. First, Section 7.2.1 examines the general problem of converting formulas to a complete highlevel program.
7.2.1 Program Synthesis
Like its ability to generate TEX code, Mathematica can convert a formula into its equivalent C or Fortran code via the CForm and FortranForm commands, respectively. With a highlevel programming language comes the additional problem of linking routines to libraries. Code generated by Mathematica will have the correct syntax, but it may not contain valid function calls. For example, CForm t^2 produces Powert,2. In order for such a statement to compile properly, one must either de ne a C macro to preprocess Power into the proper function say pow or provide a library function called Power. Alternately, one can program Mathematica to generate the proper function call. In our example, we want CForm to translate the power function in Mathematica to the pow function in the C math library:
Unprotect Power, pow ; Clear pow ; Format Power, CForm Protect Power, pow ; := pow;
Now, CForm t^2 would return powt,2. Even once we provide all of these missing hooks" to the target language, we can only generate lines of code that will compile properly. We still do not have the ability to generate subroutines, and therefore, we cannot create complete programs from
106
Generating Complete Ptolemy Simulations
scratch. SINAPSE, developed at the Schlumberger Laboratory for Computer Science, is an example of a system that uses Mathematica to generate complete programs 79, 80, 81 . SINAPSE translates a set of partial di erential wave equations with boundary conditions into nite di erence approximations. From the nite di erence approximations, SINAPSE generates optimized Fortran77, Connection MachineTM Fortran, or C source code. Complete source code generation takes on the order of 10 minutes on a Sun SparcStation 2. SINAPSE directs the user in choosing the wave phenomena, the boundary conditions, the nite di erence technique, the target language, and the target machine either a sequential computer or the parallel Connection Machine.
7.2.2 Converting Algebraic Formulas to Working Ptolemy Simulations
As previously mentioned, we would like to translate signal processing algorithms expressed as mathematical formulas to the Ptolemy environment so that the MDSPPs can t into a rapid prototyping process. Ptolemy has both a graphical user interface and a command line interpreter. The command line interpreter is the actual target for code generation. Ptolemy represents systems using block diagrams. Unlike a highlevel language, Ptolemy does not support nested function calls. Therefore, we must unravel nested calls in an algebraic formula to decompose it into a sequence of simple block operations. A simple block operation consists of a set of inputs, one function call, and a set of outputs. In Mathematica, the PtolemyProgram command converts an algebraic formula to a complete Ptolemy simulation; its algorithm is given in Figure 7.1. The order and syntax of the information given to the PtolemyProgram routine mimics the calling sequence of the various plotting commands in Mathematica:
routine expression, fv , vmin , vmaxg, options
107
Generating Complete Ptolemy Simulations
The preprocessing converts impulse responses of lters that appear in convolution operators into ltering operations and maintains the names of the original constants so that they will appear as the same name in the Ptolemy code. The Ptolemy code can then be edited at a future time to insert new values for the constants. Since the generated Ptolemy program will ultimately be interpreted and run, we must specify a header for the program that will de ne the signals and systems commonly used in the MDSPPs but missing in Ptolemy's library. We provide a default header simply called header.pt" which will either be loaded by the Ptolemy program as its rst step or be copied verbatim to the beginning of the Ptolemy program the default. In the rst case, the global parameter $PtolemyProlog must be set to
StringJoin "load "", "SignalProcessing`ObjectOriented`header.pt", $Path " "" , FindFile
The default behavior allows the output of the PtolemyProgram routine to be piped to a computer possibly on a separate network that runs Ptolemy. Ptolemy runs simulations by incrementing a oatingpoint counter over a range of values. As a part of the code generation, we de ne the counter and connect it to the variables over which to run the simulation. The primary part of the conversion process decomposes the algebraic expression into simple block operations and then translates simple block operations into a sequence of Ptolemy commands. Once the block operations have been converted into Ptolemy form, we add Ptolemy code to run the simulation. The example to demonstrate the code generation ability is the twochannel nonuniform lter bank of Figure 5.6. Figure 7.2 draws the schematic of the lter bank, gives the Mathematica code representing the system, and adds the Mathematica commands necessary to invoke the simulation. In order for the simulation to run, we must specify the analysis synthesis lters and the input signal x n . We use the lters designed in 82 so that the lter bank achieves near perfect reconstruction of
108
Generating Complete Ptolemy Simulations
Given an algebraic expression parsed into tree form with variables as leaves and operators as nodes a list of each fvariable, minimum value, maximum value, incrementg over which to run the simulation the increment if not speci ed defaults to 1 optionally, a list of values to assign to constants that appear in the algebraic expression the assignment will take place in the generated Ptolemy code Preprocess the algebraic expression by converting LSI convolution operations to ltering operations extracting constants from the algebraic expression constants are symbols that the user has not speci ed as variables Generate a complete Ptolemy simulation by determining the Ptolemy code to use as the prolog generating code that will allocate constant terms appearing in the expression, assigning to them the value given by the user see above or one if a value was not speci ed generating code that will allocate variables generating code that will de ne the algorithm by converting the algebraic expression to block diagram form converting the block diagram form to Ptolemy code generating code to invoke the simulation Figure 7.1: Algorithm to Convert Algebraic Expressions to Ptolemy Simulations
109
Generating Complete Ptolemy Simulations

"n 2

h n

0
n 3

"n 3

g n

0
n 2
xn
+
?
6
xn ^


h n

1
n 3

"n 3

g n
1
a Flow graph
upperchannel = Downsample 2,n Convolve n g0 n , Upsample 3,n Downsample 3,n Convolve n h0 n , Upsample 2,n x n lowerchannel = Convolve n g1 n , Upsample 3,n Downsample 3,n Convolve n h1 n , x n
b Representation in the new environment
h0 n h1 n g0 n g1 n x n = FIR n, Hold ReadList "ptolemy h0", Number = FIR n, Hold ReadList "ptolemy h1", Number = FIR n, Hold ReadList "ptolemy g0", Number = FIR n, Hold ReadList "ptolemy g1", Number = Cos 2 Pi n 3 Sinc Pi n 6 3; ; ; ; ;
PtolemyProgram
upperchannel + lowerchannel, n, 1, 100 , Dialogue  All "!interpreter"
c Speci cation of lters and input signal
Figure 7.2: Ptolemy Simulation of Filter Bank Run From Within Mathematica
110
Generating Complete Ptolemy Simulations
the input signal. The lter coe cients will be read by Ptolemy when it runs the simulation the Hold command prevents Mathematica from evaluating the ReadList command which would read the les in before code generation. We want to choose the input signal x n to test the reconstruction properties of the lter bank. The lter bank decomposes the input frequency band into , ; for the upper channel and ,; , ; for the lower channel. So, one suitable choice for x n is a bandpass signal having frequency content in the range of frequencies , ; , ; . By setting the amplitude to be 1 over these bands, we can ask the MDSPPs to generate the time response of the bandpass signal. as shown in Figure 7.3. The intermediate block diagram form, which unravels the nested operations in the algorithm, is shown in Figure 7.4. The generated Ptolemy code is listed in Appendix A.1.
2 3 2 3 2 3 2 3 5 6 1 2 1 2 5 6
7.2.3 Code Generation After Algorithm Rearrangement
If the analysis synthesis lters are implemented as either FIR or IIR lters, then the MDSPPs can nd better implementations of the twochannel lter bank. Using the default cost function, which counts the number of additions, multiplication, and memory elements, a much better implementation results after rewriting each analysis and synthesis channel in its polyphase form:
upperchannel = PolyphaseResample 2, FIR n, Hold ReadList "ptolemy g0",Number PolyphaseResample 3, FIR n, Hold ReadList "ptolemy h0",Number x n ; , 3, n , 2, n
lowerchannel = PolyphaseUpsample 3, FIR n, Hold ReadList "ptolemy g1",Number PolyphaseDownsample 2, FIR n, Hold ReadList "ptolemy h1",Number x n ; , n , n
By speci ng the same input
x n
as in Figure 7.2c, the Ptolemy simulation is
111
Generating Complete Ptolemy Simulations
In 14 := freqResponse = CPulse Pi 3, w + 5 Pi 6 Pi Out 14 = CPulse  + w Pi 3 2
+ CPulse Pi 3,
w  Pi 2
5 Pi + CPulse  + w Pi 3 6 freqResponse, w, n
In 15 := timeResponse = InvDTFTransform
2 I 3 n Pi n Pi 2 I 3 n Pi n Pi E Sinc E Sinc 6 6 Out 15 =  + 6 6 In 16 := SPSimplify timeResponse, Variables n
2 n Pi n Pi Cos  Sinc 3 6 Out 16 = 3
Figure 7.3: Deriving the Input Bandpass Signal for the Filter Bank Simulation
112
Generating Complete Ptolemy Simulations
rationalconst1 := 1 3 timesconst1 := 2 Pi 3 timesconst2 := Pi 6 timesbyconstant1 := n timesconst1 cos1 := Cos timesbyconstant1 timesbyconstant2 := n timesconst2 sinc1 := Sinc timesbyconstant2 times1 := cos1 sinc1 timesbyconstant3 := rationalconst1 times1 upsample1 := Upsample timesbyconstant3 2,n fir1 := FIR n, Hold ReadList ptolemy h0, Number downsample1 := Downsample fir1 3,n upsample2 := Upsample downsample1 3,n fir2 := FIR n, Hold ReadList ptolemy g0, Number downsample2 := Downsample fir2 2,n fir3 := FIR n, Hold ReadList ptolemy h1, Number timesbyconstant3 downsample3 := Downsample fir3 3,n upsample3 := Upsample downsample3 3,n fir4 := FIR n, Hold ReadList ptolemy g1, Number plus1 := downsample2 + fir4
upsample1
upsample2
upsample3
Figure 7.4: Block Diagram Form of the TwoChannel Filter Bank
113
Summary
invoked by
PtolemyProgram upperchannel + lowerchannel, Dialogue All n, 1, 100 ,
"!interpreter"
The Ptolemy code for this more e cient lter bank structure is given in Appendix A.2.
7.3 Summary
This chapter demonstrates the ability of the MDSPPs to convert an algorithm composed of a set of mathematical formulas into either a set of equivalent lines of TEX code for typesetting or a complete Ptolemy program. Ptolemy program generation is necessary for the environment to take part in rapid prototyping because the algorithm can easily be embedded into a complex system. The complex system, once simulated in Ptolemy, can then be translated into a C program or DSP assembly language program. Until recently, Ptolemy has only o ered limited abilities to process images 15 . Now, it supports multidimensional scheduling 16 so generating code for multidimensional multirate algorithms should be possible.
114
CHAPTER 8 Interactive Design of TwoDimensional Decimation Systems
This chapter demonstrates the ability of the multidimensional signal processing p