**Part 2**

## **Hardware and Photonics Applications**

94 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

Beckwith, Buck & Marangoni, (1982), Mechanical Measurements, Third Edition, 1982. Daly, James; Riley, William; McConnell, Kennet, (1984), Instrumentation for Engineering

Doebelin, E., (2003), Meaurement Systems: Application and Design, 5th Edition, McGraw

Herceg, Ed, (2006), Factors to Consider in Selecting and Specifying LVDT for applications,

Herceg, Edward, (1972), Handbook of Measurement and Control: An authoritative treatise on the theory and application of the LVDT, Schaevitz Engineering, 1972. Mishra, SK & Panda, G, (2006), A novel method for designing LVDT and its comparison

Mishra, SK; Panda, G; Das, DP; Pattanaik, SK & Meher, MR, (2005), A novel method of

Popović, Dobrivoje; Vlacic, Ljubo, (1999), Mechatronics in Engineering Design and Product

Syulski, J.K., Sykulska, E. and Hughes, S.T., (1992), Applications of Finite Element

Wu, Shang-The; Mo, Szu-Chieh; Wu, Bo-Siou, (2008), An LVDT-based self-actuating

with conventional design, Proceedings of the 2006 IEEE Sensors Applications

designing LVDT using artificial neural network, 2005 International Conference on Intelligent Sensing and Information Processing Proceedings, page 223-227, 2005. Morris, Alan S., (2001), Measurement & Instrumentation Principles, Elsevier Butterworth

Modelling in LVDT Design, The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, Vol. 11, No. 1, 73-76, James &

displacement transducer, Science Direct, Sensors and Actuators A 141 (2008) 558–

Measurements, 2nd Edition, 1984, John Wiley & Sons.

• "Mohammad Hussam ALDein" Omar Ali Bolad

Nikkei Electronics Asia, March 2006.

Symposium, pages: 129-134, 2006.

Development, CRC Press, 1999.

James Science Publishers Ltd.

Heinemann, 2001.

564.

• Anas Imad Ahmad Farouqa

• Yazan Khalid Talib Al-diqs

Hill, 2003.

**11. References** 

• Anas Abdullah Ahmad Al-shoubaki • Baha' Abd-Elrazzaq Ahmad Alsuradi

**0**

**5**

Peter Malík

*Slovak Republic*

**Computational Models Designed in MATLAB to**

**Improve Parameters and Cost of Modern Chips**

Methods and techniques to design integrated circuits made a huge progress since German engineer Werner Jacobi filed a first patent for an integrated-circuit-like semiconductor amplifying device in 1949 or since Jack Kilby successfully demonstrated the first working integrated circuit on September 12, 1958. The first integrated circuits were composed of a small number of transistors and their functionality was very limited. They were designed manually without any effective calculation tools and therefore the developing phase took long time. Each design error was difficult to discover and the whole process had to be repeated when uncovered during physical prototype testing. This further increased time necessary to develop a working chip. An important milestone in integrated circuits design is the first microprocessor Intel 4004 produced in April 1971. The significance lies in the beginning of mutual bond between produced microprocessors and methods to design them. It was discovered very soon that the high computation power of microprocessors can be utilized directly by the design process to improve future microprocessors. High performance microprocessors started to be used for calculations of an estimated behavior of future prototypes and their parameters. The simple models were transformed to more complex models that closer represented the reality. More precise models increased a chance to uncover a design error in early design phases and accelerated its correction. The continual upgrading process of design methods and tools matures into computer aided design (CAD) software platforms that are able to carry out many tasks automatically which had to be done

Constant refinement of design methods, tools and technology secured the steady grow of integrated circuits internal complexity. Intel co-founder Gordon E. Moore described the trend in his 1965 paper that the number of components in integrated circuits had doubled every year from the invention of the integrated circuit in 1958 until 1965. He predicted that the trend would continue for at least ten years Moore (1965). His prediction has proved to be uncannily accurate and is now used in semiconductor industry to guide long-term planning and to set targets for research and development ITRS (2009). This long-term trend is also

The digital design process uses highly automated CAD tools that significantly reduce the time necessary to design a prototype when the design is closely specified by full set of parameters. The original specification is usually general and many parameters have to be defined later. The process specifying missing parameters is an inseparable part of feasibility study or the first design phases. It is very important to set suitable parameters because they determine the resulting prototype. More complex integrated circuits require more parameters to be defined

**1. Introduction**

manually in the past.

known as Moore's law.

*Institute of Informatics, Slovak Academy of Sciences*

### **Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips**

#### Peter Malík

*Institute of Informatics, Slovak Academy of Sciences Slovak Republic*

#### **1. Introduction**

Methods and techniques to design integrated circuits made a huge progress since German engineer Werner Jacobi filed a first patent for an integrated-circuit-like semiconductor amplifying device in 1949 or since Jack Kilby successfully demonstrated the first working integrated circuit on September 12, 1958. The first integrated circuits were composed of a small number of transistors and their functionality was very limited. They were designed manually without any effective calculation tools and therefore the developing phase took long time. Each design error was difficult to discover and the whole process had to be repeated when uncovered during physical prototype testing. This further increased time necessary to develop a working chip. An important milestone in integrated circuits design is the first microprocessor Intel 4004 produced in April 1971. The significance lies in the beginning of mutual bond between produced microprocessors and methods to design them. It was discovered very soon that the high computation power of microprocessors can be utilized directly by the design process to improve future microprocessors. High performance microprocessors started to be used for calculations of an estimated behavior of future prototypes and their parameters. The simple models were transformed to more complex models that closer represented the reality. More precise models increased a chance to uncover a design error in early design phases and accelerated its correction. The continual upgrading process of design methods and tools matures into computer aided design (CAD) software platforms that are able to carry out many tasks automatically which had to be done manually in the past.

Constant refinement of design methods, tools and technology secured the steady grow of integrated circuits internal complexity. Intel co-founder Gordon E. Moore described the trend in his 1965 paper that the number of components in integrated circuits had doubled every year from the invention of the integrated circuit in 1958 until 1965. He predicted that the trend would continue for at least ten years Moore (1965). His prediction has proved to be uncannily accurate and is now used in semiconductor industry to guide long-term planning and to set targets for research and development ITRS (2009). This long-term trend is also known as Moore's law.

The digital design process uses highly automated CAD tools that significantly reduce the time necessary to design a prototype when the design is closely specified by full set of parameters. The original specification is usually general and many parameters have to be defined later. The process specifying missing parameters is an inseparable part of feasibility study or the first design phases. It is very important to set suitable parameters because they determine the resulting prototype. More complex integrated circuits require more parameters to be defined

Digital circuit design is composed of many steps that continually shape the final integrated circuit. The technology advancement paved the road to very complex integrated circuits called system on a chip (SoC). The more complex the chip is the more design steps it requires. Considerable interest has been oriented toward computational models with lower abstraction levels. This resulted in the development of precise computational models that are integrated in modern professional design tools. Easy access to these models is a great advantage. They are utilized mostly during moderate and finalizing design steps. This chapter is oriented to computational models on higher abstraction levels that are used in feasibility studies and first

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 99

At the beginning of design process, the computational models depend considerably on a selection of design characteristics. In this design phase, the designers create more detailed specifications and develop the guidelines to reach them. The original specification is usually too general and its parametric goals can be reached by different means. Each selected solution has to be analyzed and evaluated. Simple computational models are very important in this process. These models can differ from each other significantly. They use limited amount of parameters and many unnecessary features are omitted. Models are used in many situations, e.g., a comparison analysis of hardware implementation with software solution, analysis of different algorithms, evaluation of different hardware architectures or components. These are

Each hardware implementation performs one or several algorithms. An algorithm is evaluated by a number of operations that has to be carried out to produce a set of output data. Long computation process can produce a significant computation error that distorts the results. Each algorithm has different sensitivity to accumulation of errors. Therefore, an analysis of computation errors by models is very important and the results usually dictate minimal precision requirements of hardware components. These models are usually very simple with minimal utilization of hardware parameters. In general, all operations are transformed to basic mathematical operations that the selected hardware is capable to perform. Each transformed operation uses rounding or truncating operation to model limited

Mathematical operations can be performed with data in different data representations. Selecting an appropriate data representation can improve the resulted parameters. This effect can be modeled and analyzed. Data can be represented by floating point or fixed point. Fixed point representations are also called integer representations. Floating point has better precision for wide data intervals but if data variations are small an integer is a better solution Inacio & Ombres (1996). An integer does not contain information about decimal point placement and therefore it has to be managed externally. Floating point computation components are more complex and thus larger area of a final chip is taken Kwon et al. (2005). The increased area means higher cost of integrated circuits; hence, integer components utilization can reduce the price of a final chip. These are universal guidelines but the optimal

Models that analyze data variations are the good starting point to decide between floating point and integer. The level of specialization is the key information in the decision process. More universal applications usually use floating point due to more difficult estimation of data variations. In highly specialized applications integer can be a better solution because internal data structure can be highly optimized. The optimization results depend considerably on the

Integer is more susceptible to data overflow than floating point Beauchamp et al. (2008). Data overflow significantly distorts actual values and therefore digital designs have to include some form of protection mechanism. Before it can be designed, the knowledge where and

few examples with the most significant influence on the following design steps.

design steps.

precision of hardware components.

results can be achieved only by deep analyses.

quality of used models and on the optimization extent.

and therefore more analysis and simulation work is necessary. Usually this process requires more time in comparison to the physical design that is mostly done by modern CAD tools.

The modern CAD tools include models that closely represent the resulting prototype; however, their calculation is too much time consuming. The feasibility study and first design stages require different models. The models have to be relatively simple with the acceptable correlation to reality. An excellent software environment to create and evaluate these models is MATLAB.

Users can create new functions and procedures in MATLAB similarly as in a common programing language environment. An advantage is that MATLAB has a more simple and intuitive syntax and thorough help with many practical examples. Additionally, it has an extensive database of internal functions that can be used directly in the source code of new functions. Models usually describe a simulated behavior by mathematical equations or some iterative process. Both types of model description can be written in MATLAB intuitively with the short source code. However, the biggest advantage is very fast and simple work with calculated results and the ability to visualize the resulted data with internal MATLAB functions.

A hardware implementation of an algorithm or computation scheme implies new challenges that have to be overcome. The general aim is to design an integrated circuit with minimal area that produces results with a pre-specified overall precision. Deep analysis and lot of simulations are needed to attain this. MATLAB can be used as a part of this optimization process. The advantages are extensive mathematical support with the option to create own optimization code and easy way to sort and evaluate excessive calculated results. A simple software computational model can be created in MATLAB and later used within optimization routines. An example of this optimization process is presented in the following sections with an integrated circuit for acceleration of time consuming computations of forward and backward modified discrete cosine transform (MDCT) that is used in audio compression and coding. MDCT is used in many audio coding standards for time-to-frequency transformation of digital signals. It is perfect reconstruction cosine-modulated filter bank based on the concept of time domain aliasing cancellation Princen et al. (1987). MDCT is the basic processing block for high-quality audio compression in the international audio coding standards and commercial audio compression algorithms in consumers electronics Bosi & Goldberg (2003); Spanias et al. (2007). Forward and backward MDCT computations are the most computationally intensive operations in the well-known audio coding standards. Thus, an efficient implementation of the MDCT has become the key technology to realize real-time and low-cost audio decoders.

#### **2. Computational models in digital design process**

A computational model is a set of equations that describe some selected attribute or feature in a more simplified way compared to its true physical nature. The idea behind creating a computational model is to be able to run complex analyses and simulations of specific scenarios on PC. The fast computational speed and the ability to calculate with a wide spectrum of input data variations is very important. A computation power of common PC has been increased significantly in recent years, which enabled to run usual simulations on ordinary PC that required big server clusters in the past. The ability to create a computational model easily and in short time accelerates the design process in general. The result is that a designer can verify his/her designs in early developing phases and use the saved time to further improvements.

2 Will-be-set-by-IN-TECH

and therefore more analysis and simulation work is necessary. Usually this process requires more time in comparison to the physical design that is mostly done by modern CAD tools. The modern CAD tools include models that closely represent the resulting prototype; however, their calculation is too much time consuming. The feasibility study and first design stages require different models. The models have to be relatively simple with the acceptable correlation to reality. An excellent software environment to create and evaluate these models

Users can create new functions and procedures in MATLAB similarly as in a common programing language environment. An advantage is that MATLAB has a more simple and intuitive syntax and thorough help with many practical examples. Additionally, it has an extensive database of internal functions that can be used directly in the source code of new functions. Models usually describe a simulated behavior by mathematical equations or some iterative process. Both types of model description can be written in MATLAB intuitively with the short source code. However, the biggest advantage is very fast and simple work with calculated results and the ability to visualize the resulted data with internal MATLAB

A hardware implementation of an algorithm or computation scheme implies new challenges that have to be overcome. The general aim is to design an integrated circuit with minimal area that produces results with a pre-specified overall precision. Deep analysis and lot of simulations are needed to attain this. MATLAB can be used as a part of this optimization process. The advantages are extensive mathematical support with the option to create own optimization code and easy way to sort and evaluate excessive calculated results. A simple software computational model can be created in MATLAB and later used within optimization routines. An example of this optimization process is presented in the following sections with an integrated circuit for acceleration of time consuming computations of forward and backward modified discrete cosine transform (MDCT) that is used in audio compression and coding. MDCT is used in many audio coding standards for time-to-frequency transformation of digital signals. It is perfect reconstruction cosine-modulated filter bank based on the concept of time domain aliasing cancellation Princen et al. (1987). MDCT is the basic processing block for high-quality audio compression in the international audio coding standards and commercial audio compression algorithms in consumers electronics Bosi & Goldberg (2003); Spanias et al. (2007). Forward and backward MDCT computations are the most computationally intensive operations in the well-known audio coding standards. Thus, an efficient implementation of the MDCT has

become the key technology to realize real-time and low-cost audio decoders.

A computational model is a set of equations that describe some selected attribute or feature in a more simplified way compared to its true physical nature. The idea behind creating a computational model is to be able to run complex analyses and simulations of specific scenarios on PC. The fast computational speed and the ability to calculate with a wide spectrum of input data variations is very important. A computation power of common PC has been increased significantly in recent years, which enabled to run usual simulations on ordinary PC that required big server clusters in the past. The ability to create a computational model easily and in short time accelerates the design process in general. The result is that a designer can verify his/her designs in early developing phases and use the saved time to

**2. Computational models in digital design process**

further improvements.

is MATLAB.

functions.

Digital circuit design is composed of many steps that continually shape the final integrated circuit. The technology advancement paved the road to very complex integrated circuits called system on a chip (SoC). The more complex the chip is the more design steps it requires. Considerable interest has been oriented toward computational models with lower abstraction levels. This resulted in the development of precise computational models that are integrated in modern professional design tools. Easy access to these models is a great advantage. They are utilized mostly during moderate and finalizing design steps. This chapter is oriented to computational models on higher abstraction levels that are used in feasibility studies and first design steps.

At the beginning of design process, the computational models depend considerably on a selection of design characteristics. In this design phase, the designers create more detailed specifications and develop the guidelines to reach them. The original specification is usually too general and its parametric goals can be reached by different means. Each selected solution has to be analyzed and evaluated. Simple computational models are very important in this process. These models can differ from each other significantly. They use limited amount of parameters and many unnecessary features are omitted. Models are used in many situations, e.g., a comparison analysis of hardware implementation with software solution, analysis of different algorithms, evaluation of different hardware architectures or components. These are few examples with the most significant influence on the following design steps.

Each hardware implementation performs one or several algorithms. An algorithm is evaluated by a number of operations that has to be carried out to produce a set of output data. Long computation process can produce a significant computation error that distorts the results. Each algorithm has different sensitivity to accumulation of errors. Therefore, an analysis of computation errors by models is very important and the results usually dictate minimal precision requirements of hardware components. These models are usually very simple with minimal utilization of hardware parameters. In general, all operations are transformed to basic mathematical operations that the selected hardware is capable to perform. Each transformed operation uses rounding or truncating operation to model limited precision of hardware components.

Mathematical operations can be performed with data in different data representations. Selecting an appropriate data representation can improve the resulted parameters. This effect can be modeled and analyzed. Data can be represented by floating point or fixed point. Fixed point representations are also called integer representations. Floating point has better precision for wide data intervals but if data variations are small an integer is a better solution Inacio & Ombres (1996). An integer does not contain information about decimal point placement and therefore it has to be managed externally. Floating point computation components are more complex and thus larger area of a final chip is taken Kwon et al. (2005). The increased area means higher cost of integrated circuits; hence, integer components utilization can reduce the price of a final chip. These are universal guidelines but the optimal results can be achieved only by deep analyses.

Models that analyze data variations are the good starting point to decide between floating point and integer. The level of specialization is the key information in the decision process. More universal applications usually use floating point due to more difficult estimation of data variations. In highly specialized applications integer can be a better solution because internal data structure can be highly optimized. The optimization results depend considerably on the quality of used models and on the optimization extent.

Integer is more susceptible to data overflow than floating point Beauchamp et al. (2008). Data overflow significantly distorts actual values and therefore digital designs have to include some form of protection mechanism. Before it can be designed, the knowledge where and

The simple model is an advantage when used under optimization routines with many computation loops. Many mutually dependent parameters can easily rise the number of optimization cycles over the acceptable level. In these cases, the very simple models can produce more data on parameters' dependency that is later used to reduce the number of optimization dimensions and widths of input data intervals. This procedure can be repeated

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 101

Hardware implementation has special preferences for computational algorithms. The number and complexity of computation operations are most important. An algorithm with smaller number of operations computes results in shorter time which represents improved computational speed. The number of computation operations reduced only by single operation can significantly improve overall performance if the computation is run many times under a loop or in more complex algorithm. Not all operations are equal. Addition and subtraction operations are similar in their complexity. Multiplication operation is much more complex and therefore consumes much more area on a chip. There is an exception. Multiplication by a constant that equals to power of 2 (2, 4, 8, 16, 32, 64, ...) is the same as shifting the original binary number to the left by the number of bits that equals to base 2 logarithm of the constant. It is similar to a constant that equals to negative power of 2 (1/2, 1/4, 1/8, ...) with the difference that shifting is to the right. Hence, one removed multiplication is more significant than a removed addition or subtraction, and multiplication by power of 2

The algorithm selected for hardware implementation is presented in Britanak & Rao (2002). The algorithm computes MDCT that is used in audio compression. It has low number of mathematical operations. Additionally, it has a symmetric computational structure that is an advantage in digital design. The symmetric structure can be implemented more efficiently. The resulting hardware is dedicated to accelerate the time consuming MDCT computations in

Correct verifications during the whole design process save time. Frequent verifications facilitate a process of locating errors, because only few modifications were made between two following inspections. An algorithm verification can be done by another algorithm that is widely accepted and well proven. The best choice for verification of MDCT algorithm is the

original MDCT algorithm Princen & Bradley (1986). Its mathematical equations are:

*N* 2 

> *N* 2

where {*ck*} are MDCT coefficients and {*x*ˆ*n*} represents the time domain aliased data sequence recovered by backward MDCT. The verification model was designed straightforward by simple computation loops. The source code is shown in the following MATLAB example.

(2*k* + 1)

(2*k* + 1)

, *<sup>k</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*

<sup>2</sup> <sup>−</sup> 1, (1)

, *n* = 0, 1, . . . , *N* − 1, (2)

2*n* + 1 +

2*n* + 1 +

later with more robust and more precise models.

**3. MDCT computational models in MATLAB**

can be neglected.

audio compression.

*ck* =

*<sup>x</sup>*ˆ*<sup>n</sup>* <sup>=</sup> <sup>2</sup> *N*

Z = zeros(1, N/2); for k = 0:((N/2) - 1) sum = 0.0;

for n = 0:(N - 1)

*N*−1 ∑ *n*=0

> *N* <sup>2</sup> −1 ∑ *k*=0

function Z = MDCT\_analytic(X, N)

*xncos <sup>π</sup>* 2*N*

> *ckcos <sup>π</sup>* 2*N*

when the data overflow happens is necessary. This is usually analyzed with simple models that work with wide data intervals and the point of interest is oriented to data that produce highest peak values. It is also important to monitor the highest negative peaks, because they can cause data underflow which has similar effect as data overflow. The data producing these peaks are later used as test vectors for verifications of created designs.

An important part of modern integrated circuit is memory and memory controller. The cost of memory has been falling in last years, which changed the focus from memory size to memory controller bandwidth. The bandwidth depends on the speed of access to memory, on memory delays and on memory interface bit-width. Access to memory can be improved by rearranging of data in the memory. Reading data from memory that are placed next to each other is much faster than reading data that are randomly distributed through the whole memory. Grouping data that are often accessed at the same time and placing them in the shared memory block as close as possible can significantly reduce the access time and increase the bandwidth. Optimal data rearranging can be modeled and analyzed.

Another improvement of the memory bandwidth can be achieved by analysis of data values. If the data values are from a narrow data range, the universal exponent can be transferred separately and following data are represented only with a mantissa. This methodology is mostly used in serial data transfers and with integer representations. In this case, the better bandwidth utilization is increasing complexity of memory controller and therefore it has to be evaluated more thoroughly. A similar principle can be applied to optimization of integer representation itself.

Modern algorithms often require large memory. Lot of data are stored for a short time only. If these time periods are not overlapped, the data can be rewritten and memory size can be reduced. In many cases, the time periods are overlapped only partially and small changes within an order of operations can eliminate overlapping. Of course, this requires more complex models and deeper analysis. The result is represented by a significant reduction of memory size.

Digital designs can utilize different computation flows. A design with few universal computation blocks depends more on a memory subsystem. The advantage is a simple design and support for many different algorithms. The optimization objectives are mostly oriented to a memory subsystem and to the number of universal computational blocks. Deeper analysis can evaluate the necessary complexity level of used computational blocks. An exchange of several computational blocks with simplified blocks can reduce the final area of a chip.

A highly specialized digital design has usually a custom hardware architecture that gone through a long optimization process. Its significant part is done at a higher abstraction level. The important subject of analysis is the effect of increased parallel structures or utilization of computational pipeline. Precision and overflow analysis is more the necessity than an option. Each computational block can use unique interconnection, special features and attached latch or register. The aim is to utilize these options efficiently. Then, the resulting design has excellent parameters. The time consuming optimizations are costly and still economical in large production quantities.

During the computational model design, it is important to select an area of interest and to be able to verify precision of results. A model that produces results that differ from reality too significantly has to be corrected immediately. In many situations, the simple model composed of few mathematical equations can be sufficient, e.g., to verify an algorithm or some independent feature. Mutually dependent features have to be modeled by more robust models. The models that mimic hardware architectures with high precision give designers a chance to analyze the hardware behavior before it is actually produced. This saves time and reduces the overall cost.

The simple model is an advantage when used under optimization routines with many computation loops. Many mutually dependent parameters can easily rise the number of optimization cycles over the acceptable level. In these cases, the very simple models can produce more data on parameters' dependency that is later used to reduce the number of optimization dimensions and widths of input data intervals. This procedure can be repeated later with more robust and more precise models.

#### **3. MDCT computational models in MATLAB**

4 Will-be-set-by-IN-TECH

when the data overflow happens is necessary. This is usually analyzed with simple models that work with wide data intervals and the point of interest is oriented to data that produce highest peak values. It is also important to monitor the highest negative peaks, because they can cause data underflow which has similar effect as data overflow. The data producing these

An important part of modern integrated circuit is memory and memory controller. The cost of memory has been falling in last years, which changed the focus from memory size to memory controller bandwidth. The bandwidth depends on the speed of access to memory, on memory delays and on memory interface bit-width. Access to memory can be improved by rearranging of data in the memory. Reading data from memory that are placed next to each other is much faster than reading data that are randomly distributed through the whole memory. Grouping data that are often accessed at the same time and placing them in the shared memory block as close as possible can significantly reduce the access time and increase the bandwidth. Optimal

Another improvement of the memory bandwidth can be achieved by analysis of data values. If the data values are from a narrow data range, the universal exponent can be transferred separately and following data are represented only with a mantissa. This methodology is mostly used in serial data transfers and with integer representations. In this case, the better bandwidth utilization is increasing complexity of memory controller and therefore it has to be evaluated more thoroughly. A similar principle can be applied to optimization of integer

Modern algorithms often require large memory. Lot of data are stored for a short time only. If these time periods are not overlapped, the data can be rewritten and memory size can be reduced. In many cases, the time periods are overlapped only partially and small changes within an order of operations can eliminate overlapping. Of course, this requires more complex models and deeper analysis. The result is represented by a significant reduction

Digital designs can utilize different computation flows. A design with few universal computation blocks depends more on a memory subsystem. The advantage is a simple design and support for many different algorithms. The optimization objectives are mostly oriented to a memory subsystem and to the number of universal computational blocks. Deeper analysis can evaluate the necessary complexity level of used computational blocks. An exchange of several computational blocks with simplified blocks can reduce the final area of a chip. A highly specialized digital design has usually a custom hardware architecture that gone through a long optimization process. Its significant part is done at a higher abstraction level. The important subject of analysis is the effect of increased parallel structures or utilization of computational pipeline. Precision and overflow analysis is more the necessity than an option. Each computational block can use unique interconnection, special features and attached latch or register. The aim is to utilize these options efficiently. Then, the resulting design has excellent parameters. The time consuming optimizations are costly and still economical in

During the computational model design, it is important to select an area of interest and to be able to verify precision of results. A model that produces results that differ from reality too significantly has to be corrected immediately. In many situations, the simple model composed of few mathematical equations can be sufficient, e.g., to verify an algorithm or some independent feature. Mutually dependent features have to be modeled by more robust models. The models that mimic hardware architectures with high precision give designers a chance to analyze the hardware behavior before it is actually produced. This saves time and

peaks are later used as test vectors for verifications of created designs.

data rearranging can be modeled and analyzed.

representation itself.

of memory size.

large production quantities.

reduces the overall cost.

Hardware implementation has special preferences for computational algorithms. The number and complexity of computation operations are most important. An algorithm with smaller number of operations computes results in shorter time which represents improved computational speed. The number of computation operations reduced only by single operation can significantly improve overall performance if the computation is run many times under a loop or in more complex algorithm. Not all operations are equal. Addition and subtraction operations are similar in their complexity. Multiplication operation is much more complex and therefore consumes much more area on a chip. There is an exception. Multiplication by a constant that equals to power of 2 (2, 4, 8, 16, 32, 64, ...) is the same as shifting the original binary number to the left by the number of bits that equals to base 2 logarithm of the constant. It is similar to a constant that equals to negative power of 2 (1/2, 1/4, 1/8, ...) with the difference that shifting is to the right. Hence, one removed multiplication is more significant than a removed addition or subtraction, and multiplication by power of 2 can be neglected.

The algorithm selected for hardware implementation is presented in Britanak & Rao (2002). The algorithm computes MDCT that is used in audio compression. It has low number of mathematical operations. Additionally, it has a symmetric computational structure that is an advantage in digital design. The symmetric structure can be implemented more efficiently. The resulting hardware is dedicated to accelerate the time consuming MDCT computations in audio compression.

Correct verifications during the whole design process save time. Frequent verifications facilitate a process of locating errors, because only few modifications were made between two following inspections. An algorithm verification can be done by another algorithm that is widely accepted and well proven. The best choice for verification of MDCT algorithm is the original MDCT algorithm Princen & Bradley (1986). Its mathematical equations are:

$$\mathbf{c}\_{k} = \sum\_{n=0}^{N-1} \mathbf{x}\_{n} \cos \left[ \frac{\pi}{2N} \left( 2n + 1 + \frac{N}{2} \right) (2k+1) \right], \quad k = 0, 1, \ldots, \frac{N}{2} - 1,\tag{1}$$

$$\pounds\_{n} = \frac{2}{N} \sum\_{k=0}^{\frac{N}{2}-1} c\_k \cos \left[ \frac{\pi}{2N} \left( 2n + 1 + \frac{N}{2} \right) (2k+1) \right], \quad n = 0, 1, \ldots, N-1,\tag{2}$$

where {*ck*} are MDCT coefficients and {*x*ˆ*n*} represents the time domain aliased data sequence recovered by backward MDCT. The verification model was designed straightforward by simple computation loops. The source code is shown in the following MATLAB example.

```
function Z = MDCT_analytic(X, N)
Z = zeros(1, N/2);
for k = 0:((N/2) - 1)
    sum = 0.0;
    for n = 0:(N - 1)
```
where **I***N*/2 is the identity matrix, **J***N*/2 is the reflection matrix, all of order *N*/2. The matrix

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 103

Equation (4) can be decomposed into several stand-alone operations that are easily transformed to spare matrices. The first matrix **X2***N*/2×*<sup>N</sup>* represents subtraction in the brackets of equation (4). The second matrix **<sup>G</sup>***N*×*N*/2 represents multiplication by cosines and sines constants. The third matrix **G2***N*/2×*<sup>N</sup>* represents the remaining subtraction and addition. The matrix **X2***N*/2×*<sup>N</sup>* models a set of subtracters, matrix **<sup>G</sup>***N*×*N*/2 models a set of multipliers and

, **G2***N*/2×*<sup>N</sup>* =

*cos <sup>π</sup>* 2*N*

*sin <sup>π</sup>* 2*N*

*cos <sup>π</sup>* 2*N*

*sin <sup>π</sup>* 2*N*

where **0***N*/4 is the zero matrix, **I***N*/4 is the identity matrix, **J***N*/4 is the reflection matrix, all of order *<sup>N</sup>*/4. The matrices **X2***N*/2×*N*, **<sup>G</sup>***N*×*N*/2 and **G2***N*/2×*<sup>N</sup>* written in MATLAB are shown in

�

*cos* <sup>3</sup>*<sup>π</sup>* 2*N* ...

*sin* <sup>3</sup>*<sup>π</sup>* 2*N* ...

**I***N*/4 **0***N*/4 −**I***N*/4 **0***N*/4 **<sup>0</sup>***N*/4 **<sup>I</sup>***N*/4 **<sup>0</sup>***N*/4 **<sup>I</sup>***N*/4 �

> *cos*(*N*/2−1)*<sup>π</sup>* 2*N sin* (*N*/2−1)*<sup>π</sup>* 2*N*

, (9)

, (10)

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

**X1***N*×*<sup>N</sup>* written in MATLAB is shown in the following example.

fliplr(eye(N/2)) -eye(N/2) ];

**G2***N*/2×*<sup>N</sup>* models a set of adders and subtracters. They are given by

...

*cos* <sup>3</sup>*<sup>π</sup>* 2*N*

*sin* <sup>3</sup>*<sup>π</sup>* 2*N*

eye(N/4) zeros(N/4) -eye(N/4) zeros(N/4)];

zeros(N/4) eye(N/4) zeros(N/4) eye(N/4) ];

G = [ diag(cos((((N/2)-(1:2:N/2))\*pi)/(2\*N))) zeros(N/4);

zeros(N/4) diag(fliplr(cos((((N/2)-(1:2:N/2))\*pi)/(2\*N)))); zeros(N/4) fliplr(diag(sin((((N/2)-(1:2:N/2))\*pi)/(2\*N)))); flipud(diag(sin((((N/2)-(1:2:N/2))\*pi)/(2\*N)))) zeros(N/4) ];

...

X2 = [zeros(N/4) -eye(N/4) zeros(N/4) eye(N/4);

G2 = [ eye(N/4) zeros(N/4) -eye(N/4) zeros(N/4);

*sin* (*N*/2−1)*<sup>π</sup>* 2*N*

**0***N*/4 −**I***N*/4 **0***N*/4 **I***N*/4 **<sup>I</sup>***N*/4 **<sup>0</sup>***N*/4 <sup>−</sup>**I***N*/4 **<sup>0</sup>***N*/4 �

> *cos*(*N*/2−1)*<sup>π</sup>* 2*N*

X1 = [ eye(N/2) fliplr(eye(N/2));

**X2***N*/2×*<sup>N</sup>* =

�

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

**<sup>G</sup>***N*×*N*/2 =

the following example.

arg = pi / (2 \* N) \* ((2 \* n + 1 + N/2) \* (2 \* k + 1)); sum = sum + (X(n + 1) \* cos(arg)); end Z(k + 1) = sum;

Complete formulas constituting the selected MDCT algorithm are as follows:

$$\begin{split} z\_{2k} &= (-1)^k \frac{\sqrt{2}}{2} \sum\_{n=0}^{N-1} \left( a\_n \cos \left[ \frac{\pi (2n+1)k}{2(N/4)} \right] - b\_n \sin \left[ \frac{\pi (2n+1)k}{2(N/4)} \right] \right), \\ z\_{2k + \frac{N}{2}} &= (-1)^{k + \frac{N}{2}} \frac{\sqrt{2}}{2} \sum\_{n=0}^{N-1} (-1)^{n+1} \left( a\_n \sin \left[ \frac{\pi (2n+1)k}{2(N/4)} \right] + b\_n \cos \left[ \frac{\pi (2n+1)k}{2(N/4)} \right] \right), \\ k &= 0, 1, \ldots, \frac{N}{4} - 1, \end{split} \tag{3}$$

where

$$\begin{aligned} a\_n &= \left(\mathbf{x}\_n' - \mathbf{x}\_{\frac{n}{2}-1-n}'\right) \cos \frac{\pi(2n+1)}{2N} - \left(\mathbf{x}\_n'' - \mathbf{x}\_{\frac{n}{2}-1-n}'\right) \sin \frac{\pi(2n+1)}{2N}, \\ b\_n &= \left(\mathbf{x}\_n' - \mathbf{x}\_{\frac{n}{2}-1-n}'\right) \sin \frac{\pi(2n+1)}{2N} + \left(\mathbf{x}\_n'' - \mathbf{x}\_{\frac{n}{2}-1-n}'\right) \cos \frac{\pi(2n+1)}{2N}, \quad n = 0, 1, \dots, \frac{N}{4} - 1, \end{aligned} \tag{4}$$

and

$$\mathbf{x}\_{\mathrm{n}}^{\prime} = \mathbf{x}\_{\mathrm{n}} - \mathbf{x}\_{\mathrm{N}-1-\mathrm{n}\prime} \quad \mathbf{x}\_{\mathrm{n}}^{\prime} = \mathbf{x}\_{\mathrm{n}} + \mathbf{x}\_{\mathrm{N}-1-\mathrm{n}\prime} \quad n = 0, 1, \ldots, \frac{\mathrm{N}}{2} - 1. \tag{5}$$

Final MDCT coefficients are obtained as

$$c\_{2k} = z\_{2k}, \quad c\_{2k+1} = -z\_{N-2-2k}, \quad k = 0, 1, \ldots, \frac{N}{4} - 1. \tag{6}$$

The model of selected algorithm is used for verification and optimization of an internal computational structure and it is also used under optimization routines with many loops. Therefore, the model has to be designed more efficiently in comparison to the algorithm verification model. MATLAB can calculate matrix equations more efficiently than equations written in linear code. A model described by matrix equations can reduce computation time significantly. This is very important when the model is a part of optimization process with many loops, because the saved computation time is multiplied by a number of loops. Equations (5) can be transformed to matrix form

$$\mathbf{x}' = \mathbf{X}'\_{N/2 \times N} \mathbf{x}, \quad \mathbf{x}'' = \mathbf{X}''\_{N/2 \times N} \mathbf{x}, \tag{7}$$

where matrices **X**� *<sup>N</sup>*/2×*<sup>N</sup>* and **<sup>X</sup>**�� *<sup>N</sup>*/2×*<sup>N</sup>* are spare matrices that compute addition or subtraction of two input data. They can be merged to single matrix **X1***N*×*N*. These matrices model a set of adders and subtracters and can be used as a part of hardware architecture model. They are given by

$$\mathbf{X}'\_{N/2 \times N} = \begin{bmatrix} \mathbf{I}\_{N/2} \ -\mathbf{J}\_{N/2} \end{bmatrix}, \quad \mathbf{X}''\_{N/2 \times N} = \begin{bmatrix} \mathbf{I}\_{N/2} \ \mathbf{J}\_{N/2} \end{bmatrix}, \quad \mathbf{X}1\_{N \times N} = \begin{bmatrix} \mathbf{I}\_{N/2} \ \mathbf{J}\_{N/2} \\ \mathbf{J}\_{N/2} - \mathbf{I}\_{N/2} \end{bmatrix}, \tag{8}$$

where **I***N*/2 is the identity matrix, **J***N*/2 is the reflection matrix, all of order *N*/2. The matrix **X1***N*×*<sup>N</sup>* written in MATLAB is shown in the following example.

$$\begin{array}{rcl} \mathtt{X1} &=& [\mathtt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\langle\}}}}}}}}}}}}}}}}}}}]}}\mathtt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\langle\beta}}}}}}}}}}}}}}}}\new{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}}}}}}\new{\texttt{\texttt{\texttt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}}}}\new{\texttt{\texttt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}}\new{\texttt{\texttt{\texttt{\vdots}}}}\mathtt{X}}\end{array}}{\mathtt{\texttt{\texttt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}}}\mathtt{\texttt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}}\mathtt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\texttt{\vdots}}}}\mathtt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\langle\beta\rangle}}}\mathtt{\texttt{\texttt{\langle\beta\rangle}}}\mathtt{\texttt{\langle\beta\rangle}}\mathtt{\texttt{\langle\beta\rangle}}\end{array}}} ]\mathtt{\texttt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\langle\beta\rangle}}}\mathtt{\texttt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\langle\beta\rangle}}}}\mathtt{\texttt{\texttt{\$$

6 Will-be-set-by-IN-TECH

sum = sum + (X(n + 1) \* cos(arg));

Complete formulas constituting the selected MDCT algorithm are as follows:

(−1)*n*+<sup>1</sup>

*π*(2*n* + 1)*k* 2(*N*/4)

> *ansin*

 *x* �� *<sup>n</sup>* − *x* � *N* <sup>2</sup> −1−*n* 

 *x* �� *<sup>n</sup>* − *x* � *N* <sup>2</sup> −1−*n cos*

��

*<sup>c</sup>*2*<sup>k</sup>* <sup>=</sup> *<sup>z</sup>*2*k*, *<sup>c</sup>*2*k*+<sup>1</sup> <sup>=</sup> <sup>−</sup>*zN*−2−2*k*, *<sup>k</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*

many loops, because the saved computation time is multiplied by a number of loops.

*<sup>N</sup>*/2×*N***x**, **<sup>x</sup>**

�� = **X**��

*<sup>N</sup>*/2×*<sup>N</sup>* <sup>=</sup> [**I***N*/2 **<sup>J</sup>***N*/2] , **X1***N*×*<sup>N</sup>* <sup>=</sup>

of two input data. They can be merged to single matrix **X1***N*×*N*. These matrices model a set of adders and subtracters and can be used as a part of hardware architecture model. They are

The model of selected algorithm is used for verification and optimization of an internal computational structure and it is also used under optimization routines with many loops. Therefore, the model has to be designed more efficiently in comparison to the algorithm verification model. MATLAB can calculate matrix equations more efficiently than equations written in linear code. A model described by matrix equations can reduce computation time significantly. This is very important when the model is a part of optimization process with

end Z(k + 1) = sum;

> *z*2*k*<sup>+</sup> *<sup>N</sup>* 2

where

*an* = *x* � *<sup>n</sup>* − *x* �� *N* <sup>2</sup> −1−*n cos*

*bn* = *x* � *<sup>n</sup>* − *x* �� *N* <sup>2</sup> −1−*n* 

and

*<sup>z</sup>*2*<sup>k</sup>* = (−1)*<sup>k</sup>*

= (−1)*k*<sup>+</sup> *<sup>N</sup>*

*<sup>k</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*

*x* �

where matrices **X**�

*<sup>N</sup>*/2×*<sup>N</sup>* <sup>=</sup>

given by

**X** �

Final MDCT coefficients are obtained as

√2 2

> 4 √2 2

*N* <sup>4</sup> −1 ∑ *n*=0

 *ancos* 

*N* <sup>4</sup> −1 ∑ *n*=0

*π*(2*n* + 1) 2*N* −

<sup>2</sup>*<sup>N</sup>* <sup>+</sup>

*sinπ*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>)

*<sup>n</sup>* = *xn* − *xN*−1−*n*, *<sup>x</sup>*

Equations (5) can be transformed to matrix form

*<sup>N</sup>*/2×*<sup>N</sup>* and **<sup>X</sup>**��

**I***N*/2 −**J***N*/2

**x** � = **X** �

 , **X** ��

arg = pi / (2 \* N) \* ((2 \* n + 1 + N/2) \* (2 \* k + 1));

− *bnsin*

*π*(2*n* + 1)*k* 2(*N*/4)

*<sup>n</sup>* <sup>=</sup> *xn* <sup>+</sup> *xN*−1−*n*, *<sup>n</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*

<sup>4</sup> <sup>−</sup> 1, (3)

*sinπ*(2*<sup>n</sup>* <sup>+</sup> <sup>1</sup>) <sup>2</sup>*<sup>N</sup>* ,

*π*(2*n* + 1)

*π*(2*n* + 1)*k* 2(*N*/4)

+ *bncos*

 ,

<sup>2</sup>*<sup>N</sup>* , *<sup>n</sup>* <sup>=</sup> 0, 1, . . . , *<sup>N</sup>*

*<sup>N</sup>*/2×*N***x**, (7)

**I***N*/2 **J***N*/2 **J***N*/2 −**I***N*/2 , (8)

*<sup>N</sup>*/2×*<sup>N</sup>* are spare matrices that compute addition or subtraction

*π*(2*n* + 1)*k* 2(*N*/4)

 ,

> <sup>4</sup> <sup>−</sup> 1, (4)

<sup>2</sup> <sup>−</sup> 1. (5)

<sup>4</sup> <sup>−</sup> 1. (6)

Equation (4) can be decomposed into several stand-alone operations that are easily transformed to spare matrices. The first matrix **X2***N*/2×*<sup>N</sup>* represents subtraction in the brackets of equation (4). The second matrix **<sup>G</sup>***N*×*N*/2 represents multiplication by cosines and sines constants. The third matrix **G2***N*/2×*<sup>N</sup>* represents the remaining subtraction and addition. The matrix **X2***N*/2×*<sup>N</sup>* models a set of subtracters, matrix **<sup>G</sup>***N*×*N*/2 models a set of multipliers and **G2***N*/2×*<sup>N</sup>* models a set of adders and subtracters. They are given by

where **0***N*/4 is the zero matrix, **I***N*/4 is the identity matrix, **J***N*/4 is the reflection matrix, all of order *<sup>N</sup>*/4. The matrices **X2***N*/2×*N*, **<sup>G</sup>***N*×*N*/2 and **G2***N*/2×*<sup>N</sup>* written in MATLAB are shown in the following example.

```
X2 = [zeros(N/4) -eye(N/4) zeros(N/4) eye(N/4);
      eye(N/4) zeros(N/4) -eye(N/4) zeros(N/4)];
G = [ diag(cos((((N/2)-(1:2:N/2))*pi)/(2*N))) zeros(N/4);
    zeros(N/4) diag(fliplr(cos((((N/2)-(1:2:N/2))*pi)/(2*N))));
    zeros(N/4) fliplr(diag(sin((((N/2)-(1:2:N/2))*pi)/(2*N))));
   flipud(diag(sin((((N/2)-(1:2:N/2))*pi)/(2*N)))) zeros(N/4) ];
G2 = [ eye(N/4) zeros(N/4) -eye(N/4) zeros(N/4);
     zeros(N/4) eye(N/4) zeros(N/4) eye(N/4) ];
```
d4 -d2 d5 -(d4+d5) 1 -(d4+d5) d5 -d2 d4];

It can be seen from the 9-point DCT-II example that the matrix **DC**9×<sup>9</sup> is not a spare matrix. It is composed of many nonzero constants. Most of them are not equal to power of 2; hence, shift operations are used rarely. When the computation of this matrix is carried out directly, the number of multiplications is equal to the number of matrix elements that are not equal to ±1 or 0. The number of multiplications is equal to 60. 12 of them are multiplication by 0.5 and

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 105

Many scientific papers oriented to reduction of mathematical operations within specific algorithms were published. The research in this area of interest is highly valued. The paper Britanak & Rao (2001) presents linear code that requires 10 multiplications for a computation of 9-point DCT-II and 2 of them are multiplication by 0.5. The optimized linear code needs 6-times less multiplications. This code can be transformed to a set of spare matrices that calculate results with minimized number of mathematical operations when multiplied in correct order. One version of these spare matrices is shown in the following MATLAB

therefore the number of nontrivial multiplication is 48.

DC1 = [ 0 0 0 1 0 1 0 0 0;

0 0 0 1 0 -1 0 0 0; 0 0 1 0 0 0 1 0 0; 0 0 -1 0 0 0 1 0 0; 0 1 0 0 0 0 0 1 0; 0 1 0 0 0 0 0 -1 0; 1 0 0 0 0 0 0 0 1; -1 0 0 0 0 0 0 0 1; 0 0 0 0 1 0 0 0 0];

DC2 = [ 0 0 0 0 1 0 0 0 1;

1 0 1 0 0 0 0 0 0; 0 0 1 0 0 0 -1 0 0; 1 0 0 0 0 0 -1 0 0; 1 0 -1 0 0 0 0 0 0; 0 1 0 -1 0 0 0 0 0; 0 0 0 1 0 0 0 1 0; 0 1 0 0 0 0 0 -1 0; 0 1 0 1 0 0 0 0 0; 0 0 0 0 0 -d1 0 0 0; 0 0 0 0 d2 0 0 0 0; 0 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 1];

DC3 = [ 0 1 0 0 0 0 0 0 0 0 0 1 0 0;

0 0 0 0 0 1 0 0 0 0 0 0 1 0; 0 0 -d3 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -d4 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 -d5 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 -d6 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 -d7 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 -d8 0 0 0 0 0;

example.

Equation (6) represents permutation operations that change signs and order of data. It can be transformed to the matrix form. The source code of the resulting matrix **<sup>Q</sup>***N*/2×*N*/2 is shown in the following MATLAB example.

```
function Q = permatrix(N);
Q = zeros(N/2);
for k = 1:(N/4)
Q(2*k-1, k) = 1;
Q(2*k, N/2-k+1) = -1;
end
```
The most complex part of the selected algorithm is equation (3). It represents N/4-point DCT-II and N/4-point DST-II combined by a butterfly operation Britanak & Rao (2002). DCT-II and DST-II can be represented by a matrix but it cannot be generated universally for any N/4 due to their not symmetrical internal computational structure. The MDCT computational model is optimized for audio standard MP3 that uses 36-point MDCT for slowly changing audio signal and 12-point MDCT for rapidly changing audio signal. Hence, the model requires 9-point DCT-II/DST-II and 3-point DCT-II/DST-II matrices. The matrices DCT-II and DST-II are similar and it is possible to transform one to the other with simple modifications. The procedure to transforming DST-II to DCT-II is described in Britanak (2002). The transformation is composed of extra operations that invert each odd input data and reverse order of all input data. It results in one DCT-II that is used two times with different sets of data. This is advantageous when dedicated hardware architecture is to be designed. The single hardware block can be optimized within shorter time and overall hardware architecture can be reduced easily by multiplexing input data from two datasets to single hardware block. The 9-point DCT-II matrix can be derived directly from the mathematic equation. The result is a matrix that represents many operations of addition and multiplication. The number of these operations can be reduced significantly by matrix decomposition to several spare matrices or by optimized linear code. The both principles utilize computation with partial results that influence several output data. It is more efficient to use already computed partial results during the following computations than to calculate all results from input data only. The 9-point DCT-II matrix **DC**9×<sup>9</sup> written in MATLAB is shown in the following example.


8 Will-be-set-by-IN-TECH

Equation (6) represents permutation operations that change signs and order of data. It can be transformed to the matrix form. The source code of the resulting matrix **<sup>Q</sup>***N*/2×*N*/2 is shown

The most complex part of the selected algorithm is equation (3). It represents N/4-point DCT-II and N/4-point DST-II combined by a butterfly operation Britanak & Rao (2002). DCT-II and DST-II can be represented by a matrix but it cannot be generated universally for any N/4 due to their not symmetrical internal computational structure. The MDCT computational model is optimized for audio standard MP3 that uses 36-point MDCT for slowly changing audio signal and 12-point MDCT for rapidly changing audio signal. Hence, the model requires 9-point DCT-II/DST-II and 3-point DCT-II/DST-II matrices. The matrices DCT-II and DST-II are similar and it is possible to transform one to the other with simple modifications. The procedure to transforming DST-II to DCT-II is described in Britanak (2002). The transformation is composed of extra operations that invert each odd input data and reverse order of all input data. It results in one DCT-II that is used two times with different sets of data. This is advantageous when dedicated hardware architecture is to be designed. The single hardware block can be optimized within shorter time and overall hardware architecture can be reduced easily by multiplexing input data from two datasets to single hardware block. The 9-point DCT-II matrix can be derived directly from the mathematic equation. The result is a matrix that represents many operations of addition and multiplication. The number of these operations can be reduced significantly by matrix decomposition to several spare matrices or by optimized linear code. The both principles utilize computation with partial results that influence several output data. It is more efficient to use already computed partial results during the following computations than to calculate all results from input data only. The

9-point DCT-II matrix **DC**9×<sup>9</sup> written in MATLAB is shown in the following example.

DC = [ 1 1 1 1 1 1 1 1 1;

d7 d1 d8 (d7-d8) 0 (d8-d7) -d8 -d1 -d7; -d3 d2 (d3+d5) -d5 -1 -d5 (d3+d5) d2 -d3; d1 0 -d1 -d1 0 d1 d1 0 -d1; -(d3+d4) -d2 d3 d4 1 d4 d3 -d2 -(d3+d4); (d7-d6) -d1 -d6 d7 0 -d7 d6 d1 (d6-d7); d2 -1 d2 d2 -1 d2 d2 -1 d2; d6 -d1 (d6+d8) -d8 0 d8 -(d6+d8) d1 -d6;

in the following MATLAB example.

function Q = permatrix(N);

Q = zeros(N/2); for k = 1:(N/4) Q(2\*k-1, k) = 1; Q(2\*k, N/2-k+1) = -1;

u=2 \* pi / 9; d1 = sqrt(3)/2;

d2 = 0.5; d3 = cos(4\*u); d4 = cos(2\*u); d5 = cos(u); d6 = sin(4\*u); d7 = sin(2\*u); d8 = sin(u);

end

$$\begin{array}{ccccccccc} \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ } & \text{\$ }$$

It can be seen from the 9-point DCT-II example that the matrix **DC**9×<sup>9</sup> is not a spare matrix. It is composed of many nonzero constants. Most of them are not equal to power of 2; hence, shift operations are used rarely. When the computation of this matrix is carried out directly, the number of multiplications is equal to the number of matrix elements that are not equal to ±1 or 0. The number of multiplications is equal to 60. 12 of them are multiplication by 0.5 and therefore the number of nontrivial multiplication is 48.

Many scientific papers oriented to reduction of mathematical operations within specific algorithms were published. The research in this area of interest is highly valued. The paper Britanak & Rao (2001) presents linear code that requires 10 multiplications for a computation of 9-point DCT-II and 2 of them are multiplication by 0.5. The optimized linear code needs 6-times less multiplications. This code can be transformed to a set of spare matrices that calculate results with minimized number of mathematical operations when multiplied in correct order. One version of these spare matrices is shown in the following MATLAB example.


zeros(14,9) DC2 ];

zeros(11,14) DC3 ];

zeros(14,11) DC4 ];

zeros(9,14) DC5 ];

matrix written in MATLAB is shown in the following example.

D = [ 1, zeros(1,N/2-1);

C = nc\*Q\*D\*DCS5\*DCS4\*DCS3\*DCS2\*DCS1\*G2\*G\*X2\*X1\*x;

these matrices also improves the resulting hardware architecture.

process which is based on solving these mathematical equations. The equation

The results from DCT-II and DST-II are combined by a butterfly operation. It is a symmetrical operation and therefore can be described universally for any *<sup>N</sup>*/4 by a matrix **<sup>D</sup>***N*/2×*<sup>N</sup>*/2. The

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 107

zeros(N/4-1,1) eye(N/4-1) zeros(N/4-1,1) -fliplr(eye(N/4-1));

zeros(N/4-1,1) -fliplr(eye(N/4-1)) zeros(N/4-1,1) -eye(N/4-1)];

The presented matrices are the main elements of the MDCT model. The model computation core is composed of multiplications of these matrices in a correct order. The MDCT model is

All matrices used in the model are spare matrices with low number of mathematical operations which increases computational speed in MATLAB. They also correspond to basic hardware computational blocks as an adder, subtracter or multiplier. Further optimization of

The MDCT computational model is created. The next step is analysis of data variations. Its results represent actual intervals of data inside the modeled computational structure. The positive and negative extreme input values can cause data overflow or data underflow and therefore it is important to understand how the extreme values are propagated through the computational structure. This knowledge facilitates a selection of appropriate data representations and is necessary when integer representations are selected. Extreme values of partial results can be paired with input data that caused them. These input data can be included into a verification data set and used later to test if overflowing or underflowing can happen. A computational model described by a matrix equation simplifies the analysis

zeros(1,N/4) -1 zeros(1,N/4-1);

**<sup>y</sup>***X*<sup>1</sup> <sup>=</sup> **X1***N*×*N***x**, (11)

DCS3 = [ DC3 zeros(11,14);

DCS4 = [ DC4 zeros(14,11);

DCS5 = [ DC5 zeros(9,14);

shown in the following MATLAB example.

function C = MDCT36(x);

N = 36;

d2 = 0.5; d3 = cos(4\*u); d4 = cos(2\*u); d5 = cos(u); d6 = sin(4\*u); d7 = sin(2\*u); d8 = sin(u);

nc = sqrt(2)/2; u=2 \* pi / 9; d1 = sqrt(3)/2;


The matrices **DC1**9×9, **DC2**14×9, **DC3**11×14, **DC4**14×<sup>11</sup> and **DC5**9×<sup>14</sup> use the same constants as the matrix **DC**9×<sup>9</sup> from the previous example. They are designed to model hardware architecture. Addition or subtraction operations have only two inputs which is characteristic for hardware implementations. Mutual addition of more input data is distributed through several matrices. The last few lines of the several matrices perform no mathematical operation. These lines provide access to input data or partial results for the following matrices.

Computation of 9-point DST-II with 9-point DCT-II requires two small modifications. Inversion of odd input data can be modeled by a modified identity matrix that has all odd diagonal elements equal to −1. Reverse order of input data can be modeled by a reflection matrix. These two matrices are used together with **DC1**9×<sup>9</sup> during computation of 9-point DST-II. The computations of 9-point DCT-II and 9-point DST-II are done independently with clearly separated input data and results. The model with both 9-point DCT-II and 9-point DST-II is shown in the following MATLAB example.

```
M1 = [ fliplr(eye(9)), zeros(9);
          zeros(9), eye(9)];
M2 = diag([ones(1,9), repmat([1 -1],1,4), 1 ]);
DCC1 = [ DC1 zeros(9);
       zeros(9) DC1 ];
DCS1 = DCC1*M2*M1;
DCS2 = [ DC2 zeros(14,9);
```

```
zeros(14,9) DC2 ];
DCS3 = [ DC3 zeros(11,14);
      zeros(11,14) DC3 ];
DCS4 = [ DC4 zeros(14,11);
      zeros(14,11) DC4 ];
DCS5 = [ DC5 zeros(9,14);
      zeros(9,14) DC5 ];
```
10 Will-be-set-by-IN-TECH

0 0 0 0 0 0 0 0 0 0 -1 0 0 1; 1 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 1 0 0 0 0];

DC4 = [ d2 0 0 0 0 0 0 0 0 0 0;

0 -d1 0 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 1 0 0; 0 0 -1 0 0 0 0 0 1 0 0; 0 0 0 1 0 0 0 0 1 0 0; 0 0 0 0 0 1 0 0 0 0 1; 0 0 0 0 0 -1 0 0 0 0 1; 0 0 0 0 0 0 1 0 0 0 1; 1 0 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 0 1 0; 0 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 1 0 0 0 0 0 0; 0 0 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 1 0 0 0];

DC5 = [ 0 0 0 0 0 0 0 0 1 0 0 0 0 0;

DST-II is shown in the following MATLAB example.

zeros(9), eye(9)]; M2 = diag([ones(1,9), repmat([1 -1],1,4), 1 ]);

M1 = [ fliplr(eye(9)), zeros(9);

zeros(9) DC1 ];

DCS2 = [ DC2 zeros(14,9);

DCC1 = [ DC1 zeros(9);

DCS1 = DCC1\*M2\*M1;

0 0 0 0 0 0 0 -1 0 0 0 0 0 1; 0 0 -1 0 0 0 0 0 0 0 0 1 0 0; 0 1 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0 0 -1 0 0 0; 0 0 0 0 0 0 1 0 0 0 0 0 -1 0; 1 0 0 0 0 0 0 0 0 -1 0 0 0 0; 0 0 0 0 0 1 0 0 0 0 0 0 0 1; 0 0 0 0 1 0 0 0 0 0 0 1 0 0];

The matrices **DC1**9×9, **DC2**14×9, **DC3**11×14, **DC4**14×<sup>11</sup> and **DC5**9×<sup>14</sup> use the same constants as the matrix **DC**9×<sup>9</sup> from the previous example. They are designed to model hardware architecture. Addition or subtraction operations have only two inputs which is characteristic for hardware implementations. Mutual addition of more input data is distributed through several matrices. The last few lines of the several matrices perform no mathematical operation.

Computation of 9-point DST-II with 9-point DCT-II requires two small modifications. Inversion of odd input data can be modeled by a modified identity matrix that has all odd diagonal elements equal to −1. Reverse order of input data can be modeled by a reflection matrix. These two matrices are used together with **DC1**9×<sup>9</sup> during computation of 9-point DST-II. The computations of 9-point DCT-II and 9-point DST-II are done independently with clearly separated input data and results. The model with both 9-point DCT-II and 9-point

These lines provide access to input data or partial results for the following matrices.

The results from DCT-II and DST-II are combined by a butterfly operation. It is a symmetrical operation and therefore can be described universally for any *<sup>N</sup>*/4 by a matrix **<sup>D</sup>***N*/2×*<sup>N</sup>*/2. The matrix written in MATLAB is shown in the following example.

```
D = [ 1, zeros(1,N/2-1);
   zeros(N/4-1,1) eye(N/4-1) zeros(N/4-1,1) -fliplr(eye(N/4-1));
                 zeros(1,N/4) -1 zeros(1,N/4-1);
   zeros(N/4-1,1) -fliplr(eye(N/4-1)) zeros(N/4-1,1) -eye(N/4-1)];
```
The presented matrices are the main elements of the MDCT model. The model computation core is composed of multiplications of these matrices in a correct order. The MDCT model is shown in the following MATLAB example.

```
function C = MDCT36(x);
N = 36;
nc = sqrt(2)/2;
u=2 * pi / 9;
d1 = sqrt(3)/2;
d2 = 0.5;
d3 = cos(4*u);
d4 = cos(2*u);
d5 = cos(u);
d6 = sin(4*u);
d7 = sin(2*u);
d8 = sin(u);
C = nc*Q*D*DCS5*DCS4*DCS3*DCS2*DCS1*G2*G*X2*X1*x;
```
All matrices used in the model are spare matrices with low number of mathematical operations which increases computational speed in MATLAB. They also correspond to basic hardware computational blocks as an adder, subtracter or multiplier. Further optimization of these matrices also improves the resulting hardware architecture.

The MDCT computational model is created. The next step is analysis of data variations. Its results represent actual intervals of data inside the modeled computational structure. The positive and negative extreme input values can cause data overflow or data underflow and therefore it is important to understand how the extreme values are propagated through the computational structure. This knowledge facilitates a selection of appropriate data representations and is necessary when integer representations are selected. Extreme values of partial results can be paired with input data that caused them. These input data can be included into a verification data set and used later to test if overflowing or underflowing can happen. A computational model described by a matrix equation simplifies the analysis process which is based on solving these mathematical equations. The equation

$$\mathbf{y}^{X1} = \mathbf{X1}\_{N \times N} \mathbf{x}\_{\prime} \tag{11}$$

The value of the least significant bit (LSB) decreases exponentially with linear increasing of integer interface bit-width. When the LSB value is too small the information that it carries can be discarded without increasing the overall computation error. This can be modeled by introduction of rounding or truncating operations into the model. MATLAB contains internal functions that transform data to binary integer and back. This is shown in the following MATLAB example where b1 is an integer variable and a2 is recovered value of a1 with finite

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 109

The functions used in the previous example take too much time to compute. A faster method to model finite precision of data is the utilization of truncating operations that directly represents discarding of several LSBs. The rounding operation is more complex and requires additional computations. The rounding operation is usually transformed to a combination of addition and truncating operation where the added value is equal to the half of new LSB after rounding. Data in hardware are usually represented in two's-complement arithmetic and discarding of several LSBs in this system represents rounding toward zero. The function that performs rounding toward zero in MATLAB is called **fix()**. The operation of the **fix()** function is confined to the placement of decimal point. The shift of decimal point modifies the amount of data that is going to be discarded. The new data has to be in original form and this implies a second shift of the decimal point in opposite direction. The combination of shift, rounding toward zero and opposite shift can model the truncating operation in hardware. A shift is represented by multiplication of a constant. A constant in the binary shift is equal to power of 2. The positive power represents data shift to the left and negative power stands for data shift to the right. This is shown in the following MATLAB example that is a modification of

precision set by q.

a1 = 0.54321;

q = quantizer([16,15]);

the MDCT computational model.

fyX1 = fix(2^s1 \* yX1) \* 2^(-s1);

fyX2 = fix(2^s2 \* yX2) \* 2^(-s2);

fyG2 = fix(2^s4 \* yG2) \* 2^(-s4);

fyDCS1 = fix(2^s5 \* yDCS1) \* 2^(-s5);

fyDCS2 = fix(2^s6 \* yDCS2) \* 2^(-s6);

fyDCS3 = fix(2^s7 \* yDCS3) \* 2^(-s7);

fyDCS4 = fix(2^s8 \* yDCS4) \* 2^(-s8);

fyDCS5 = fix(2^s9 \* yDCS5) \* 2^(-s9);

fyD = fix(2^s10 \* yD) \* 2^(-s10);

fyG = fix(2^s3 \* yG) \* 2^(-s3);

yX1 = X1\*x;

yX2 = X2\*fyX1;

yG = G\*fyX2;

yG2 = G2\*fyG;

yDCS1 = DCS1\*fyG2;

yDCS2 = DCS2\*fyDCS1;

yDCS3 = DCS3\*fyDCS2;

yDCS4 = DCS4\*fyDCS3;

yDCS5 = DCS5\*fyDCS4;

yD = D\*fyDCS5;

b1 = num2bin(q, a1); a2 = bin2num(q, b1);

can be transformed to a set of equations

$$\begin{aligned} y\_1^{\mathbf{X}1} &= \mathbf{X1}(1, \cdot)\mathbf{x} = \mathbf{X1}(1, 1)\mathbf{x}\_1 + \mathbf{X1}(1, 2)\mathbf{x}\_2 + \dots \mathbf{X1}(1, N)\mathbf{x}\_{N\prime} \\ y\_2^{\mathbf{X}1} &= \mathbf{X1}(2, \cdot)\mathbf{x} = \mathbf{X1}(2, 1)\mathbf{x}\_1 + \mathbf{X1}(2, 2)\mathbf{x}\_2 + \dots \mathbf{X1}(2, N)\mathbf{x}\_{N\prime} \end{aligned}$$

$$\vdots$$

$$y\_N^{\mathbf{X1}} = \mathbf{X1}(N, \cdot)\mathbf{x} = \mathbf{X1}(N, 1)\mathbf{x}\_1 + \mathbf{X1}(N, 2)\mathbf{x}\_2 + \dots \mathbf{X1}(N, N)\mathbf{x}\_{N\prime} \tag{12}$$

where **X1**(*i*, :) is the i-th line of the matrix **X1***N*×*N*. The positive extreme value of *<sup>y</sup>X*<sup>1</sup> *<sup>i</sup>* is produced by addition of partial results **X1**(*i*, *j*)*xi* with maximal values that are all positive. This means that if **X1**(*i*, *j*) is positive then *xi* should be maximal value within the interval of input data. If it is negative then *xi* should be maximal negative value within the interval of input data. Similarly, the negative extreme value of *yX*<sup>1</sup> *<sup>i</sup>* is produced by addition of partial results **X1**(*i*, *j*)*xi* with maximal values that are all negative. The matrix **X1***N*×*<sup>N</sup>* is very simple with only additions and subtractions and therefore it is intuitive that maximal value of results is double the maximal input value. This simple conclusion is valid only at the beginning of computation. In the middle, it is influenced by the internal constants and computational structure.

The computation of extreme values that are produced by other matrices is done by the same principle. The matrix **X1***N*×*<sup>N</sup>* in equation (11) is substituted by product of the matrix under investigation and other matrices that precede it. This principle is shown in equation (13). The product of necessary matrices can be calculated easily in MATLAB. The whole process can be also automated by new function created in MATLAB.

$$\begin{aligned} \mathbf{y}^{\text{XI}} &= \mathbf{M}\mathbf{X}\mathbf{2} = \mathbf{X}\mathbf{2} \star \mathbf{X}\mathbf{1} \star \mathbf{x} \\ \mathbf{y}^{\text{C}} &= \mathbf{M}\mathbf{G} = \mathbf{G} \star \mathbf{X2} \star \mathbf{X1} \star \mathbf{x} \\ \mathbf{y}^{\text{C}} &= \mathbf{M}\mathbf{G2} = \mathbf{G2} \star \mathbf{G} \star \mathbf{X2} \star \mathbf{X1} \star \mathbf{x} \\ \mathbf{y}^{\text{DCS1}} &= \mathbf{M}\mathbf{DCS1} = \mathbf{DCS1} \star \mathbf{G2} \star \mathbf{G} \star \mathbf{X2} \star \mathbf{X1} \star \mathbf{x} \\ \mathbf{y}^{\text{DCS2}} &= \mathbf{M}\mathbf{DCS2} = \mathbf{DCS2} \star \mathbf{DCS1} \star \mathbf{G2} \star \mathbf{G \star X2} \star \mathbf{X1} \star \mathbf{x} \\ \mathbf{y}^{\text{DCS3}} &= \mathbf{M}\mathbf{DCS3} = \mathbf{DCS4} \star \mathbf{DCS3} \star \mathbf{DCS1} \star \mathbf{G2} \star \mathbf{G \star X2} \star \mathbf{X1} \star \mathbf{x} \\ \mathbf{y}^{\text{DCS4}} &= \mathbf{M}\mathbf{DCS4} = \mathbf{DCS4} \star \mathbf{DCS3} \star \mathbf{DCS3} \star \mathbf{DCS2} \star \mathbf{DCS1} \star \mathbf{G2} \star \mathbf{G \star X2} \star \mathbf{X1} \star \mathbf{x} \\ \mathbf{y}^{\text{DCS5}} &= \mathbf{M}\mathbf{D} \mathbf{D} \mathbf{S}5 \star \mathbf{DCS4} \star \mathbf{DCS3} \star \mathbf{DCS2} \star \mathbf{DCS1} \star \mathbf{G2} \star \mathbf{G \star X2} \star$$

The next step is optimization of all data representations within the whole computational structure. Integer representations are selected for the MDCT model. The extreme values are known, therefore data overflow is no longer a problem. The aim of the following optimization is to set an optimal ratio between overall computation error, input data precision and output data precision. The MDCT algorithm is composed of many mathematical operations and most of them increase the interval of their results. This means widening of integer interface. An addition usually increases the interface by one bit and a multiplication increases it to double. The value of the least significant bit (LSB) decreases exponentially with linear increasing of integer interface bit-width. When the LSB value is too small the information that it carries can be discarded without increasing the overall computation error. This can be modeled by introduction of rounding or truncating operations into the model. MATLAB contains internal functions that transform data to binary integer and back. This is shown in the following MATLAB example where b1 is an integer variable and a2 is recovered value of a1 with finite precision set by q.

```
q = quantizer([16,15]);
a1 = 0.54321;
b1 = num2bin(q, a1);
a2 = bin2num(q, b1);
```
12 Will-be-set-by-IN-TECH

<sup>1</sup> = **X1**(1, :)**x** = **X1**(1, 1)*x*<sup>1</sup> + **X1**(1, 2)*x*<sup>2</sup> + ... **X1**(1, *N*)*xN*,

<sup>2</sup> = **X1**(2, :)**x** = **X1**(2, 1)*x*<sup>1</sup> + **X1**(2, 2)*x*<sup>2</sup> + ... **X1**(2, *N*)*xN*, . . .

where **X1**(*i*, :) is the i-th line of the matrix **X1***N*×*N*. The positive extreme value of *<sup>y</sup>X*<sup>1</sup>

produced by addition of partial results **X1**(*i*, *j*)*xi* with maximal values that are all positive. This means that if **X1**(*i*, *j*) is positive then *xi* should be maximal value within the interval of input data. If it is negative then *xi* should be maximal negative value within the interval of

results **X1**(*i*, *j*)*xi* with maximal values that are all negative. The matrix **X1***N*×*<sup>N</sup>* is very simple with only additions and subtractions and therefore it is intuitive that maximal value of results is double the maximal input value. This simple conclusion is valid only at the beginning of computation. In the middle, it is influenced by the internal constants and computational

The computation of extreme values that are produced by other matrices is done by the same principle. The matrix **X1***N*×*<sup>N</sup>* in equation (11) is substituted by product of the matrix under investigation and other matrices that precede it. This principle is shown in equation (13). The product of necessary matrices can be calculated easily in MATLAB. The whole process can be

*<sup>N</sup>* = **X1**(*N*, :)**x** = **X1**(*N*, 1)*x*<sup>1</sup> + **X1**(*N*, 2)*x*<sup>2</sup> + ... **X1**(*N*, *N*)*xN*, (12)

*<sup>i</sup>* is

*<sup>i</sup>* is produced by addition of partial

can be transformed to a set of equations

*yX*<sup>1</sup>

*yX*<sup>1</sup>

*yX*<sup>1</sup>

**<sup>y</sup>***X*<sup>1</sup> <sup>=</sup> **MX2** <sup>=</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup> <sup>y</sup>***<sup>G</sup>* <sup>=</sup> **MG** <sup>=</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>**

**<sup>y</sup>***G*<sup>2</sup> <sup>=</sup> **MG2** <sup>=</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>**

structure.

input data. Similarly, the negative extreme value of *yX*<sup>1</sup>

also automated by new function created in MATLAB.

**<sup>y</sup>***DCS*<sup>1</sup> <sup>=</sup> **MDCS1** <sup>=</sup> **DCS1** <sup>∗</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>**

**<sup>y</sup>***DCS*<sup>2</sup> <sup>=</sup> **MDCS2** <sup>=</sup> **DCS2** <sup>∗</sup> **DCS1** <sup>∗</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>**

**<sup>y</sup>***DCS*<sup>3</sup> <sup>=</sup> **MDCS3** <sup>=</sup> **DCS3** <sup>∗</sup> **DCS2** <sup>∗</sup> **DCS1** <sup>∗</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>**

**<sup>y</sup>***DCS*<sup>4</sup> <sup>=</sup> **MDCS4** <sup>=</sup> **DCS4** <sup>∗</sup> **DCS3** <sup>∗</sup> **DCS2** <sup>∗</sup> **DCS1** <sup>∗</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>**

**<sup>y</sup>***DCS*<sup>5</sup> <sup>=</sup> **MDCS5** <sup>=</sup> **DCS5** <sup>∗</sup> **DCS4** <sup>∗</sup> **DCS3** <sup>∗</sup> **DCS2** <sup>∗</sup> **DCS1** <sup>∗</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup> <sup>y</sup>***<sup>D</sup>* <sup>=</sup> **MD** <sup>=</sup> **<sup>D</sup>** <sup>∗</sup> **DCS5** <sup>∗</sup> **DCS4** <sup>∗</sup> **DCS3** <sup>∗</sup> **DCS2** <sup>∗</sup> **DCS1** <sup>∗</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>**

**<sup>y</sup>***<sup>Q</sup>* <sup>=</sup> **MQ** <sup>=</sup> **<sup>Q</sup>** <sup>∗</sup> **<sup>D</sup>** <sup>∗</sup> **DCS5** <sup>∗</sup> **DCS4** <sup>∗</sup> **DCS3** <sup>∗</sup> **DCS2** <sup>∗</sup> **DCS1** <sup>∗</sup> **G2** <sup>∗</sup> **<sup>G</sup>** <sup>∗</sup> **X2** <sup>∗</sup> **X1** <sup>∗</sup> **<sup>x</sup>** (13)

The next step is optimization of all data representations within the whole computational structure. Integer representations are selected for the MDCT model. The extreme values are known, therefore data overflow is no longer a problem. The aim of the following optimization is to set an optimal ratio between overall computation error, input data precision and output data precision. The MDCT algorithm is composed of many mathematical operations and most of them increase the interval of their results. This means widening of integer interface. An addition usually increases the interface by one bit and a multiplication increases it to double. The functions used in the previous example take too much time to compute. A faster method to model finite precision of data is the utilization of truncating operations that directly represents discarding of several LSBs. The rounding operation is more complex and requires additional computations. The rounding operation is usually transformed to a combination of addition and truncating operation where the added value is equal to the half of new LSB after rounding. Data in hardware are usually represented in two's-complement arithmetic and discarding of several LSBs in this system represents rounding toward zero. The function that performs rounding toward zero in MATLAB is called **fix()**. The operation of the **fix()** function is confined to the placement of decimal point. The shift of decimal point modifies the amount of data that is going to be discarded. The new data has to be in original form and this implies a second shift of the decimal point in opposite direction. The combination of shift, rounding toward zero and opposite shift can model the truncating operation in hardware. A shift is represented by multiplication of a constant. A constant in the binary shift is equal to power of 2. The positive power represents data shift to the left and negative power stands for data shift to the right. This is shown in the following MATLAB example that is a modification of the MDCT computational model.

```
yX1 = X1*x;
fyX1 = fix(2^s1 * yX1) * 2^(-s1);
yX2 = X2*fyX1;
fyX2 = fix(2^s2 * yX2) * 2^(-s2);
yG = G*fyX2;
fyG = fix(2^s3 * yG) * 2^(-s3);
yG2 = G2*fyG;
fyG2 = fix(2^s4 * yG2) * 2^(-s4);
yDCS1 = DCS1*fyG2;
fyDCS1 = fix(2^s5 * yDCS1) * 2^(-s5);
yDCS2 = DCS2*fyDCS1;
fyDCS2 = fix(2^s6 * yDCS2) * 2^(-s6);
yDCS3 = DCS3*fyDCS2;
fyDCS3 = fix(2^s7 * yDCS3) * 2^(-s7);
yDCS4 = DCS4*fyDCS3;
fyDCS4 = fix(2^s8 * yDCS4) * 2^(-s8);
yDCS5 = DCS5*fyDCS4;
fyDCS5 = fix(2^s9 * yDCS5) * 2^(-s9);
yD = D*fyDCS5;
fyD = fix(2^s10 * yD) * 2^(-s10);
```
yQ = Q\*fyD; fyQ = fix(2^s11 \* yQ) \* 2^(-s11); C = nc\*fyQ;

The variables s1, s2, . . . , s11 from the previous example represent the precision of partial results. The values of variables are natural numbers and they depend on the input data precision, input data interval and intervals of partial results. They all can be used as inputs of an optimization process where low values of the variables mean low overall computation precision and high values of the variables imply large area of a resulted hardware implementation. Some of them are more important than others. Multiplication doubles the output precision. It means that a multiplier output has two time more bits in comparison to a single multiplier input and therefore the precision of multiplication results has to be modified. Overall computation precision is influenced significantly by precision of internal constants. The MDCT algorithm uses constants three times during the whole computation. The first set of constants is represented by the matrix **<sup>G</sup>***N*×*<sup>N</sup>*/2. The second set of constants is used during computation of DCT-II and DST-II and the last multiplications by constant represent the normalization of all results by the scaling constant nc. Precision of the internal constants can be optimized together with other variables.

Suitable values of the variables and precision of constants can be set by a simple analysis. Several variables can be set intuitively. One of them is s1 that stands for the precision of matrix **X1***N*×*N*. This matrix performs addition and subtraction operations and adding one extra bit is required. The truncating operation at the beginning of computation affects overall computation error more significantly. This truncating operation can be omitted. A similar principle can be done also for matrix **X2***N*/2×*N*. The other variables can be set to the output precision that is increased by one or two bits. The effect of higher precision can be easily visualized up to three dimensions in MATLAB. In direct approach, more dimensions are not possible to show clearly. Therefore, it is useful to modify only two variables at a time when the visualization is the main optimization tool. An example of the described visualization is shown in Figures 1 and 2. As can be seen, the multiplication results with bit-width wider than 20 bits does not improve overall precision in the specific situation shown in Figures 1 and 2. The combination of visualizations and intuitive approach usually result in a set of suitable variable and parameter values within short time. However, complex optimization routines are necessary to calculate the best set of variables and parameters.

<sup>16</sup> <sup>18</sup> <sup>20</sup> <sup>22</sup> <sup>24</sup> <sup>26</sup> <sup>28</sup> <sup>30</sup> <sup>32</sup>

Bit−width of G matrix results

<sup>16</sup> <sup>18</sup> <sup>20</sup> <sup>22</sup> <sup>24</sup> <sup>26</sup> <sup>28</sup> <sup>30</sup> <sup>32</sup>

Bit−width of DCS3 matrix results

cost of increased precision cannot be assessed.

Maximal computational error

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 111

Fig. 1. Maximal computational error dependence on the bit-width of multiplication results.

The most important parameters of the hardware chip are the area, speed or frequency and delay. These parameters depend significantly on implementation technology. They cannot be calculated directly or at least without extremely complex models. Low abstraction level models are capable to compute them with high precision but they are too time demanding. They are used usually only for verification or fine-tuning during the final steps of the design process. The problem is that these parameters are necessary within the optimization process. They are important factors of final product specification. At the same time, they represent feasibility boundaries. The chip with the area larger than few hundreds square millimeters is not possible to produce and it is the same with too high frequency. Many other parameters influence their values, e.g., higher precision stands for larger area. Precision can be easily modeled and accurately calculated. However, this is not true for the area. An optimization process has to be able to evaluate the effect of increased precision to the area. Otherwise, a

Time consuming calculations can be replaced by an estimation. The advantage is in significantly reduced computation complexity. This allows to use an estimation as a part of the optimization process and to optimize other parameters more efficiently. A disadvantage is in limited accuracy. In general, estimation of some complex parameter, e.g., the area, can be implemented into computational models with any abstraction levels. Estimation capabilities of higher abstraction models produce more coarse results. However, the accuracy can be improved by increasing the amount of data that the estimation is based on. The estimation data can be taken from many diverse situations and can be represented by many different forms. The most common procedure is to take data from well known previous situations or

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x 10−3

Maximal errors

#### **4. Refinement of computational models and their application**

Computational models are used in analyses, simulations and optimization processes. They are the key elements to further improve the parameters of future products. Optimized computational models represent hardware architectures that are later implemented as hardware prototypes. Verification of hardware prototypes and their following modifications are necessary. Hardware prototypes are designed in different platforms; therefore, direct verification is not possible. Verification models have been created to test designed computational models. The mechanism to generate input stimuli and pair them with correct results can be utilized. The necessary modifications have to be made to transform generated data to the specific data representations that are used in hardware prototypes. The models can be upgraded to write the transformed data to an external file that is directly used as input of a verification process. The modification of the final computational model allows to include not only input output pairs but also data that represent partial results of the internal computational structure.

14 Will-be-set-by-IN-TECH

The variables s1, s2, . . . , s11 from the previous example represent the precision of partial results. The values of variables are natural numbers and they depend on the input data precision, input data interval and intervals of partial results. They all can be used as inputs of an optimization process where low values of the variables mean low overall computation precision and high values of the variables imply large area of a resulted hardware implementation. Some of them are more important than others. Multiplication doubles the output precision. It means that a multiplier output has two time more bits in comparison to a single multiplier input and therefore the precision of multiplication results has to be modified. Overall computation precision is influenced significantly by precision of internal constants. The MDCT algorithm uses constants three times during the whole computation. The first set of constants is represented by the matrix **<sup>G</sup>***N*×*<sup>N</sup>*/2. The second set of constants is used during computation of DCT-II and DST-II and the last multiplications by constant represent the normalization of all results by the scaling constant nc. Precision of the internal constants

Suitable values of the variables and precision of constants can be set by a simple analysis. Several variables can be set intuitively. One of them is s1 that stands for the precision of matrix **X1***N*×*N*. This matrix performs addition and subtraction operations and adding one extra bit is required. The truncating operation at the beginning of computation affects overall computation error more significantly. This truncating operation can be omitted. A similar principle can be done also for matrix **X2***N*/2×*N*. The other variables can be set to the output precision that is increased by one or two bits. The effect of higher precision can be easily visualized up to three dimensions in MATLAB. In direct approach, more dimensions are not possible to show clearly. Therefore, it is useful to modify only two variables at a time when the visualization is the main optimization tool. An example of the described visualization is shown in Figures 1 and 2. As can be seen, the multiplication results with bit-width wider than 20 bits does not improve overall precision in the specific situation shown in Figures 1 and 2. The combination of visualizations and intuitive approach usually result in a set of suitable variable and parameter values within short time. However, complex optimization routines

Computational models are used in analyses, simulations and optimization processes. They are the key elements to further improve the parameters of future products. Optimized computational models represent hardware architectures that are later implemented as hardware prototypes. Verification of hardware prototypes and their following modifications are necessary. Hardware prototypes are designed in different platforms; therefore, direct verification is not possible. Verification models have been created to test designed computational models. The mechanism to generate input stimuli and pair them with correct results can be utilized. The necessary modifications have to be made to transform generated data to the specific data representations that are used in hardware prototypes. The models can be upgraded to write the transformed data to an external file that is directly used as input of a verification process. The modification of the final computational model allows to include not only input output pairs but also data that represent partial results of the internal

yQ = Q\*fyD;

C = nc\*fyQ;

fyQ = fix(2^s11 \* yQ) \* 2^(-s11);

can be optimized together with other variables.

computational structure.

are necessary to calculate the best set of variables and parameters.

**4. Refinement of computational models and their application**

The most important parameters of the hardware chip are the area, speed or frequency and delay. These parameters depend significantly on implementation technology. They cannot be calculated directly or at least without extremely complex models. Low abstraction level models are capable to compute them with high precision but they are too time demanding. They are used usually only for verification or fine-tuning during the final steps of the design process. The problem is that these parameters are necessary within the optimization process. They are important factors of final product specification. At the same time, they represent feasibility boundaries. The chip with the area larger than few hundreds square millimeters is not possible to produce and it is the same with too high frequency. Many other parameters influence their values, e.g., higher precision stands for larger area. Precision can be easily modeled and accurately calculated. However, this is not true for the area. An optimization process has to be able to evaluate the effect of increased precision to the area. Otherwise, a cost of increased precision cannot be assessed.

Time consuming calculations can be replaced by an estimation. The advantage is in significantly reduced computation complexity. This allows to use an estimation as a part of the optimization process and to optimize other parameters more efficiently. A disadvantage is in limited accuracy. In general, estimation of some complex parameter, e.g., the area, can be implemented into computational models with any abstraction levels. Estimation capabilities of higher abstraction models produce more coarse results. However, the accuracy can be improved by increasing the amount of data that the estimation is based on. The estimation data can be taken from many diverse situations and can be represented by many different forms. The most common procedure is to take data from well known previous situations or

with a new graphical interface. MATLAB is also able to transform their source code to

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 113

The most important goal of computational models is to improve the resulting hardware implementation or some parameters of final chip. The original MDCT computational model had been continually upgraded. Modifications were oriented to simplify the computational structure, reduce mathematical operations and improve generality. Changes were done to the algorithm, order, type and number of mathematical operations, computational pipeline and

The first major modification is oriented to a simplification of the computational structure located between inputs and first multiplications. This part is described by mathematical equations (4) and (5) or by matrices **X1***N*×*N*, **X2***N*/2×*N*, **<sup>G</sup>***N*×*<sup>N</sup>*/2, **G2***N*/2×*<sup>N</sup>* of the MDCT computational model. Basic algebraic manipulations applied to these equations transform several additions of input data to additions of internal constants. The constants are not changed during the whole computation process and therefore the addition of internal constants can be calculated in advance. This produces new constants that are used as replacement of the original constants and the overall number of addition operations is

Another modification is associated with a scaling factor. The MDCT algorithm uses the scaling constant nc that is applied to all results. Algebraic manipulation can change the position of this multiplication and combine it with any matrix of the MDCT computational model. The matrix **DCS3**22×<sup>28</sup> represents multiplication by constants; however, few of them are equal to ±0.5 that can be implemented by shift operations. Combination of the scaling factor with these constants nullifies this advantage and then a standard multiplier has to be used. The best option is to use the matrix **<sup>G</sup>***N*×*N*/2 again. In this situation, the new constants are calculated as a product of scaling factor and original constants. This results in a reduced number of

The original MDCT computational structure does not differentiate between addition of both positive, combination of positive and negative and both negative input variables. Addition of both positive input variables is implemented by an adder. Addition of one positive and one negative is implemented by a subtracter. Addition of two negative input variables is implemented by a combination of inverter and subtracter. The inverter changes the sign of the first input and then the subtracter calculates the result. The problem is that this operation is not usually implemented into standard computation blocks and extra inverters increase the final area. The solution lies in further optimization of the MDCT computational structure. The inversion can be separated from the subtraction and combined with preceding mathematical operations by applying basic algebraic manipulations. The separated operation of inversion can be completely absorbed and then the extra inverter is omitted. An example of this

The sign change of the variable e in equation (15) is done by exchanging the input variables a and b before the variable e is calculated. Internal constants can also be used to completely

*e* = *a* − *b*, *f* = *c* + *d*, *g* = −*e* − *f* , (14)

*e* = −(*a* − *b*) = −*a* + *b* = *b* − *a*, *f* = *c* + *d*, *g* = *e* − *f* . (15)

**5. Modifications of the MDCT computational model to improve hardware**

independent application that can be run on a common PC.

reduced. This process is presented in Šimlaštík et al. (2006).

situation is described by equations (14) and (15).

multipliers used within MDCT computational structure Malík et al. (2006).

**architecture parameters**

level of parallel computations.

product versions and to interpolate between them. Higher number of situations or different versions imply a better interpolation and also more accurate estimation.

An estimated parameter is dependent on many other parameters. Each parameter included into estimation function increases the number of dimensions that the interpolation has to work with. An estimation can be simplified by reducing this number with selecting only the most important parameters. A simplified estimation can be used at the beginning of optimization process together with higher abstraction computational models. When more data are accumulated, additional parameters can be used to increase accuracy. Later, the new data can be taken from early prototypes which further improves accuracy.

More refined computational models can be transformed to independent applications. This can be useful when a product is a part of a more complex unit. Teams that design other parts can test mutual features and integration capability much sooner. Another example is the hardware-software co-design. The created applications can allow to test new software before the hardware is actually produced. Computational models with improved estimation capabilities can be transformed to applications that are offered together with final products as an extra tool. This is useful for universal products with several setting values or many parameters that can be modified by customers.

An independent application requires to design an input/output interface for users. The best choice is a graphical interface that is easy to understand and can be operated intuitively. MATLAB is a great tool to design computational models. The models can also be upgraded with a new graphical interface. MATLAB is also able to transform their source code to independent application that can be run on a common PC.

#### **5. Modifications of the MDCT computational model to improve hardware architecture parameters**

16 Will-be-set-by-IN-TECH

Maximal computational errors

Fig. 2. Maximal computational error dependence on the bit-width of multiplication results

product versions and to interpolate between them. Higher number of situations or different

An estimated parameter is dependent on many other parameters. Each parameter included into estimation function increases the number of dimensions that the interpolation has to work with. An estimation can be simplified by reducing this number with selecting only the most important parameters. A simplified estimation can be used at the beginning of optimization process together with higher abstraction computational models. When more data are accumulated, additional parameters can be used to increase accuracy. Later, the new

More refined computational models can be transformed to independent applications. This can be useful when a product is a part of a more complex unit. Teams that design other parts can test mutual features and integration capability much sooner. Another example is the hardware-software co-design. The created applications can allow to test new software before the hardware is actually produced. Computational models with improved estimation capabilities can be transformed to applications that are offered together with final products as an extra tool. This is useful for universal products with several setting values or many

An independent application requires to design an input/output interface for users. The best choice is a graphical interface that is easy to understand and can be operated intuitively. MATLAB is a great tool to design computational models. The models can also be upgraded

versions imply a better interpolation and also more accurate estimation.

data can be taken from early prototypes which further improves accuracy.

<sup>16</sup> <sup>18</sup> <sup>20</sup> <sup>22</sup> <sup>24</sup> <sup>26</sup> <sup>28</sup> <sup>30</sup> <sup>32</sup>

Bit−width of DCS3 matrix results

displayed in base 10 logarithmic scale.

parameters that can be modified by customers.

−5 −4.8 −4.6 −4.4 −4.2 −4 −3.8 −3.6 −3.4 −3.2 −3

Base 10 logarthm of maximal errors

<sup>16</sup> <sup>18</sup> <sup>20</sup> <sup>22</sup> <sup>24</sup> <sup>26</sup> <sup>28</sup> <sup>30</sup> <sup>32</sup>

Bit−width of G matrix results

The most important goal of computational models is to improve the resulting hardware implementation or some parameters of final chip. The original MDCT computational model had been continually upgraded. Modifications were oriented to simplify the computational structure, reduce mathematical operations and improve generality. Changes were done to the algorithm, order, type and number of mathematical operations, computational pipeline and level of parallel computations.

The first major modification is oriented to a simplification of the computational structure located between inputs and first multiplications. This part is described by mathematical equations (4) and (5) or by matrices **X1***N*×*N*, **X2***N*/2×*N*, **<sup>G</sup>***N*×*<sup>N</sup>*/2, **G2***N*/2×*<sup>N</sup>* of the MDCT computational model. Basic algebraic manipulations applied to these equations transform several additions of input data to additions of internal constants. The constants are not changed during the whole computation process and therefore the addition of internal constants can be calculated in advance. This produces new constants that are used as replacement of the original constants and the overall number of addition operations is reduced. This process is presented in Šimlaštík et al. (2006).

Another modification is associated with a scaling factor. The MDCT algorithm uses the scaling constant nc that is applied to all results. Algebraic manipulation can change the position of this multiplication and combine it with any matrix of the MDCT computational model. The matrix **DCS3**22×<sup>28</sup> represents multiplication by constants; however, few of them are equal to ±0.5 that can be implemented by shift operations. Combination of the scaling factor with these constants nullifies this advantage and then a standard multiplier has to be used. The best option is to use the matrix **<sup>G</sup>***N*×*N*/2 again. In this situation, the new constants are calculated as a product of scaling factor and original constants. This results in a reduced number of multipliers used within MDCT computational structure Malík et al. (2006).

The original MDCT computational structure does not differentiate between addition of both positive, combination of positive and negative and both negative input variables. Addition of both positive input variables is implemented by an adder. Addition of one positive and one negative is implemented by a subtracter. Addition of two negative input variables is implemented by a combination of inverter and subtracter. The inverter changes the sign of the first input and then the subtracter calculates the result. The problem is that this operation is not usually implemented into standard computation blocks and extra inverters increase the final area. The solution lies in further optimization of the MDCT computational structure. The inversion can be separated from the subtraction and combined with preceding mathematical operations by applying basic algebraic manipulations. The separated operation of inversion can be completely absorbed and then the extra inverter is omitted. An example of this situation is described by equations (14) and (15).

$$e = a - b, \quad f = c + d, \quad g = -e - f,\tag{14}$$

$$e = -(a - b) = -a + b = b - a, \quad f = c + d, \quad g = e - f. \tag{15}$$

The sign change of the variable e in equation (15) is done by exchanging the input variables a and b before the variable e is calculated. Internal constants can also be used to completely

**Complete Computation 4-input LUTs Flip flops Frequency Latency**

**Complete 4-input LUTs Flip flops Frequency Latency**

Difference -16.7 % -84.5 % 29.2 % -43.7 %

No optimization 100 % 100 % 100 % 100 %

After optimization 83.3 % 15.5 % 129.2 % 56.3 %

estimated results were significantly different in comparison to prototypes. A further analysis showed that the generated set of input parameters represents hardware architecture with

The MDCT computational model was used to calculate optimal input parameters. These parameters represent the precision of internal constants and several partial results. The necessary input variables are input precision and overall computation precision. The

Several hardware prototypes were implemented into FPGA technology. The first prototype represents implementation of MDCT hardware architecture without any optimization. It is implementation of a complete MDCT computational structure based on the selected MDCT

The second line of Table 1 represents parameters of the optimized MDCT computational structure with no reduction of parallel computational level. This prototype is able to calculate both forward and backward MDCT. The comparison of these two prototypes from Table 1 is

As can be seen in Table 2, all parameters were improved. The highest improvement is in reduction of internal memory elements, represented by flip flops, which are reduced by nearly 85 %. The latency is reduced by 44 % which is also induced by operational frequency increased by 29 %. The area of the combination logic, represented by 4-input LUTs, was reduced by

The MDCT implementations with a complete computational structure consume too large area and therefore several prototypes with a reduced level of parallel computation were implemented. The parameters of two such prototypes are shown in Table 3 and Table 4. The first lines of Table 3 and Table 4 represent parameters of the first reduced MDCT computational structure that was designed directly from the complete MDCT computational structure shown in the first line of Table 1. The reduced implementation prototype has only minimal extra optimizations. The second lines of Table 3 and Table 4 represent parameters of the reduced MDCT computational structure that was designed from the optimized MDCT computational structure with no reduction of parallel computational level shown in the

Table 2. Comparison of 36-point MDCT prototypes implemented into FPGA

acceptable overall parameters and the estimation was not upgraded further.

algorithm. The resulted parameters are shown in the first line of Table 1.

No optimization MDCT 18661 8174 32.26 MHz 341 ns

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 115

After optimization MDCT / IMDCT 15536 1269 41.67 MHz 192 ns

Table 1. 36-point MDCT prototypes implemented into FPGA

optimization process is presented in Malík et al. (2006).

**6. Evaluation of the optimization process**

shown in Table 2.

nearly 17 %.

absorbed extra inversions. In this case, new constants are equal to product of −1 and original constants which is the same as inversion of the original constants. The MDCT computational model described by matrix equations is a great advantage, because the inversion separation can be written by matrix decomposition. The new matrix that depicts separated inversions is spare diagonal matrix with most diagonal elements equal to 1 and only the elements that represent inversions under investigation equal to −1. MATLAB can easily calculate the product of this matrix with a preceding matrix. If the result does not represent a situation when extra inversion can be absorbed completely, MATLAB can calculate new matrix decomposition and then the process is repeated. Optimized MDCT hardware architecture with reduced extra inverters is presented in Malík et al. (2006).

The next major modification allows computation of backward MDCT by the computational structure optimized for forward MDCT. This means that computation of forward and backward MDCT can be calculated by the same MDCT kernel. Complexity of the required changes within MDCT computational structure is not high which was presented in Cho et al. (2004) with a different MDCT algorithm. The changes are related to the inversion of several signals and modification of the input and output parts of the MDCT computational structure. Backward MDCT has twice as many outputs than inputs which is exactly opposite in comparison to forward MDCT. The matrices at the beginning and at the end of the MDCT computational structure were modified to include this. Backward MDCT uses extra scaling factor equal to 2/*N* that was absorbed to internal constants. Modified MDCT computational structure and resulted hardware architecture is presented in Malík (2007).

The MDCT computational structure is composed of many mathematical operations. When they are implemented directly the resulted hardware consumes a large area. The advantage is in high computational speed. However, higher computational speed can be achieved by increasing operational frequency and therefore the large area is a disadvantage. The area can be scaled down by separating mathematical operations to different time periods. It is actually a reduction of parallel computational level. This process requires to add memory elements into architecture to store partial results and to keep them unchanged for several computational cycles. The symmetrical computational structure, e.g., the butterfly structure, can be scaled down intuitively. In this case, all results are described by the same equation with different input variables. This means that computational structure of a single result is able to compute other results if a mechanism that selects different input variables is attached to it. The situation is more difficult when a computational structure is not symmetrical, because the universal computational structure dedicated for several results has to be designed manually. This is usually a challenging and time consuming process.

The memory elements added into computational structure can be used to introduce a pipeline computation. A pipelined computational structure represents a computational structure that is divided into several parts and all these parts compute their partial results at the same time. The aim is to design these parts so the whole computation process is continual and none part has to wait for its input data. This results in significantly increased computational speed without increasing operational frequency. A minor disadvantage is in induced computational delay. The MDCT computational structure with reduced parallel computational level and pipelined computation is presented in Malík (2007).

The MDCT computational model was upgraded with estimation capabilities. The estimated parameter is the area described in 4-input Look-Up Tables (LUTs) which are basic computational blocks used in FPGA technology. The FPGA chips were used for hardware prototyping. Estimation data were taken from elementary prototypes. The estimation uses a linear interpolation to calculate estimated results. The estimated results of this simple model had been used within the optimization process. However, later evaluation showed that the


Table 1. 36-point MDCT prototypes implemented into FPGA

18 Will-be-set-by-IN-TECH

absorbed extra inversions. In this case, new constants are equal to product of −1 and original constants which is the same as inversion of the original constants. The MDCT computational model described by matrix equations is a great advantage, because the inversion separation can be written by matrix decomposition. The new matrix that depicts separated inversions is spare diagonal matrix with most diagonal elements equal to 1 and only the elements that represent inversions under investigation equal to −1. MATLAB can easily calculate the product of this matrix with a preceding matrix. If the result does not represent a situation when extra inversion can be absorbed completely, MATLAB can calculate new matrix decomposition and then the process is repeated. Optimized MDCT hardware architecture

The next major modification allows computation of backward MDCT by the computational structure optimized for forward MDCT. This means that computation of forward and backward MDCT can be calculated by the same MDCT kernel. Complexity of the required changes within MDCT computational structure is not high which was presented in Cho et al. (2004) with a different MDCT algorithm. The changes are related to the inversion of several signals and modification of the input and output parts of the MDCT computational structure. Backward MDCT has twice as many outputs than inputs which is exactly opposite in comparison to forward MDCT. The matrices at the beginning and at the end of the MDCT computational structure were modified to include this. Backward MDCT uses extra scaling factor equal to 2/*N* that was absorbed to internal constants. Modified MDCT computational

The MDCT computational structure is composed of many mathematical operations. When they are implemented directly the resulted hardware consumes a large area. The advantage is in high computational speed. However, higher computational speed can be achieved by increasing operational frequency and therefore the large area is a disadvantage. The area can be scaled down by separating mathematical operations to different time periods. It is actually a reduction of parallel computational level. This process requires to add memory elements into architecture to store partial results and to keep them unchanged for several computational cycles. The symmetrical computational structure, e.g., the butterfly structure, can be scaled down intuitively. In this case, all results are described by the same equation with different input variables. This means that computational structure of a single result is able to compute other results if a mechanism that selects different input variables is attached to it. The situation is more difficult when a computational structure is not symmetrical, because the universal computational structure dedicated for several results has to be designed manually.

The memory elements added into computational structure can be used to introduce a pipeline computation. A pipelined computational structure represents a computational structure that is divided into several parts and all these parts compute their partial results at the same time. The aim is to design these parts so the whole computation process is continual and none part has to wait for its input data. This results in significantly increased computational speed without increasing operational frequency. A minor disadvantage is in induced computational delay. The MDCT computational structure with reduced parallel computational level and

The MDCT computational model was upgraded with estimation capabilities. The estimated parameter is the area described in 4-input Look-Up Tables (LUTs) which are basic computational blocks used in FPGA technology. The FPGA chips were used for hardware prototyping. Estimation data were taken from elementary prototypes. The estimation uses a linear interpolation to calculate estimated results. The estimated results of this simple model had been used within the optimization process. However, later evaluation showed that the

with reduced extra inverters is presented in Malík et al. (2006).

structure and resulted hardware architecture is presented in Malík (2007).

This is usually a challenging and time consuming process.

pipelined computation is presented in Malík (2007).


Table 2. Comparison of 36-point MDCT prototypes implemented into FPGA

estimated results were significantly different in comparison to prototypes. A further analysis showed that the generated set of input parameters represents hardware architecture with acceptable overall parameters and the estimation was not upgraded further.

The MDCT computational model was used to calculate optimal input parameters. These parameters represent the precision of internal constants and several partial results. The necessary input variables are input precision and overall computation precision. The optimization process is presented in Malík et al. (2006).

#### **6. Evaluation of the optimization process**

Several hardware prototypes were implemented into FPGA technology. The first prototype represents implementation of MDCT hardware architecture without any optimization. It is implementation of a complete MDCT computational structure based on the selected MDCT algorithm. The resulted parameters are shown in the first line of Table 1.

The second line of Table 1 represents parameters of the optimized MDCT computational structure with no reduction of parallel computational level. This prototype is able to calculate both forward and backward MDCT. The comparison of these two prototypes from Table 1 is shown in Table 2.

As can be seen in Table 2, all parameters were improved. The highest improvement is in reduction of internal memory elements, represented by flip flops, which are reduced by nearly 85 %. The latency is reduced by 44 % which is also induced by operational frequency increased by 29 %. The area of the combination logic, represented by 4-input LUTs, was reduced by nearly 17 %.

The MDCT implementations with a complete computational structure consume too large area and therefore several prototypes with a reduced level of parallel computation were implemented. The parameters of two such prototypes are shown in Table 3 and Table 4.

The first lines of Table 3 and Table 4 represent parameters of the first reduced MDCT computational structure that was designed directly from the complete MDCT computational structure shown in the first line of Table 1. The reduced implementation prototype has only minimal extra optimizations. The second lines of Table 3 and Table 4 represent parameters of the reduced MDCT computational structure that was designed from the optimized MDCT computational structure with no reduction of parallel computational level shown in the

**7. Conclusion**

from 17 % to 93 %.

Bratislava.

**8. References**

final chip with better parameters and lower cost.

with calculated data further improves the effect of visualizations.

Science + Business Media, Inc., New York, NY.

The technology of integrated circuit production has been significantly improved which resulted in very high integration of currently produced chips. The advantage is in lower cost and higher computational speed of produced chips. Higher functionality represents increased complexity of the internal computational structure that has to be designed. The designed hardware architecture has to be optimized. Each small improvement that reduces the chip final area, represents lower cost of the produced chips. This is very important in high volume productions where the savings can grow to high values. Computer aided design systems represent powerful tools that facilitate and speed up the design process. Most of these tools include complex low abstraction level models that are able to optimize parameters of the final chip. However, they are not suitable for a general optimization, because they are too computation power demanding. An increase of modern chip complexity resulted in more optimization steps that use models with higher abstraction level. MATLAB is an excellent software platform that can be used to design, verify and highly optimize these computational models. The optimization with higher abstraction level models is faster and therefore wider intervals of input parameters and more suitable solutions can be evaluated. This implies the

Computational Models Designed in MATLAB to Improve Parameters and Cost of Modern Chips 117

Many computational models have to be created during the whole design process. The feasibility study and early design steps are examples where the computational models are created anew. MATLAB includes extensive database of internal functions that can be directly used into computational models which reduces the time necessary to design them. The created models have to be verified and therefore the ability to generate input test data and to analyze and visualize the calculated results is necessary. This is supported by embedded MATLAB functions and the representing functions are easy to use. Additional manipulation

The design process of computational models created in MATLAB was presented by MDCT computational models. The presented models are oriented to different tasks with a common goal to improve parameters of the resulting hardware architecture and its implementation to a chip. The effect of the optimization process was shown by the comparison of several hardware prototypes implemented into FPGA technology during different optimization phase. The parameters of the optimized prototypes are significantly improved. The improvements range

Beauchamp, M. J., Hauck, S., Underwood, K. D. & Hemmert, K. S. (2008). Architectural

Britanak, V. (2002). The refined efficient implementation of the mdct in mp3 and comparison

Britanak, V. & Rao, K. (2002). A New Fast Algorithm for the Unified Forward and Inverse

Britanak, V. & Rao, K. R. (2001). An Efficient Implementation of the Forward and Inverse

MDCT/MDST Computation, *Signal Processing* 82(3): 433–459.

IEEE Signal Processing Letters, Vol. 8, No. 10, Oct. 2001, p. 279.

*Transactions on Very Large Scale Integration (VLSI) Systems* 16(2): 177–187. Bosi, M. & Goldberg, R. E. (2003). *Introduction to Digital Audio Coding and Standards*, Springer

Modifications to Enhance the Floating-Point Performance of FPGAs, *IEEE*

with other methods, *Technical Report IISAS-2002-2*, Institute of Informatics SAS,

MDCT in MPEG Audio Coding, *IEEE Signal Processing Letters* 8(2): 48–51. Erratum:


Table 3. 36-point MDCT prototypes with reduced level of parallel computation implemented into FPGA


Table 4. 36-point MDCT prototypes with reduced level of parallel computation implemented into FPGA


Table 5. Comparison of 36-point MDCT prototypes with reduced level of parallel computation implemented into FPGA

second line of Table 1. The comparison of the two implemented prototypes with a reduced level of parallel computation is shown in Table 5.

As can be seen in Table 5, all parameters were improved. The highest improvement is in increasing operational frequency by nearly 93 %. This is partially caused by introducing a computational pipeline with 3 stages. As a result the latency is reduced by 16 % even though the number of clock cycles that represent latency is increased (see the fifth column of Table 4). Internal memory elements are reduced by 29 %. The area of the combinational logic is reduced by 40 % and additionally no tree-state buffers (TBUFs) were used. The optimized reduced prototype is able to calculate both forward and backward MDCT.

The MDCT hardware prototypes were also implemented into two ASIC technologies. Comparison between implementation into AMS 350nm CMOS standard cell library and implementation into UMC 90nm CMOS low power digital library technology with clock gating technique is presented in Malík et al. (2009).

The all MDCT computational models were designed in MATLAB. Most optimization work and time consuming calculations were also done in MATLAB, which substantially contributed to overall improvements. Additionally, many internal MATLAB functions were used to generate random input values and to visualize results, which accelerated the whole optimization process significantly.

#### **7. Conclusion**

20 Will-be-set-by-IN-TECH

No optimization MDCT 5472 2467 2070 0

After optimization MDCT / IMDCT 3275 1752 0 0

into FPGA

into FPGA

computation implemented into FPGA

level of parallel computation is shown in Table 5.

gating technique is presented in Malík et al. (2009).

optimization process significantly.

prototype is able to calculate both forward and backward MDCT.

Table 3. 36-point MDCT prototypes with reduced level of parallel computation implemented

**Reduced Frequency Clock period Computation Latency Latency**

No optimization 18.52 MHz 54 ns 18 x CLK 26 x CLK 1 404 ns

After optimization 35.71 MHz 28 ns 18 x CLK 42 x CLK 1 176 ns

Table 4. 36-point MDCT prototypes with reduced level of parallel computation implemented

No optimization 100 % 100 % 100 % 100 % 100 %

After optimization 59.9 % 71.0 % 0.0 % 192.9 % 83.8 %

second line of Table 1. The comparison of the two implemented prototypes with a reduced

As can be seen in Table 5, all parameters were improved. The highest improvement is in increasing operational frequency by nearly 93 %. This is partially caused by introducing a computational pipeline with 3 stages. As a result the latency is reduced by 16 % even though the number of clock cycles that represent latency is increased (see the fifth column of Table 4). Internal memory elements are reduced by 29 %. The area of the combinational logic is reduced by 40 % and additionally no tree-state buffers (TBUFs) were used. The optimized reduced

The MDCT hardware prototypes were also implemented into two ASIC technologies. Comparison between implementation into AMS 350nm CMOS standard cell library and implementation into UMC 90nm CMOS low power digital library technology with clock

The all MDCT computational models were designed in MATLAB. Most optimization work and time consuming calculations were also done in MATLAB, which substantially contributed to overall improvements. Additionally, many internal MATLAB functions were used to generate random input values and to visualize results, which accelerated the whole

Table 5. Comparison of 36-point MDCT prototypes with reduced level of parallel

Difference -40.1 % -29.0 % -100 % 92.9 % -16.2 %

**Reduced 4-input LUTs Flip flops TBUFs Frequency Latency**

**Reduced Computation 4-input LUTs Flip flops TBUFs Block RAM**

The technology of integrated circuit production has been significantly improved which resulted in very high integration of currently produced chips. The advantage is in lower cost and higher computational speed of produced chips. Higher functionality represents increased complexity of the internal computational structure that has to be designed. The designed hardware architecture has to be optimized. Each small improvement that reduces the chip final area, represents lower cost of the produced chips. This is very important in high volume productions where the savings can grow to high values. Computer aided design systems represent powerful tools that facilitate and speed up the design process. Most of these tools include complex low abstraction level models that are able to optimize parameters of the final chip. However, they are not suitable for a general optimization, because they are too computation power demanding. An increase of modern chip complexity resulted in more optimization steps that use models with higher abstraction level. MATLAB is an excellent software platform that can be used to design, verify and highly optimize these computational models. The optimization with higher abstraction level models is faster and therefore wider intervals of input parameters and more suitable solutions can be evaluated. This implies the final chip with better parameters and lower cost.

Many computational models have to be created during the whole design process. The feasibility study and early design steps are examples where the computational models are created anew. MATLAB includes extensive database of internal functions that can be directly used into computational models which reduces the time necessary to design them. The created models have to be verified and therefore the ability to generate input test data and to analyze and visualize the calculated results is necessary. This is supported by embedded MATLAB functions and the representing functions are easy to use. Additional manipulation with calculated data further improves the effect of visualizations.

The design process of computational models created in MATLAB was presented by MDCT computational models. The presented models are oriented to different tasks with a common goal to improve parameters of the resulting hardware architecture and its implementation to a chip. The effect of the optimization process was shown by the comparison of several hardware prototypes implemented into FPGA technology during different optimization phase. The parameters of the optimized prototypes are significantly improved. The improvements range from 17 % to 93 %.

#### **8. References**


**0**

**6**

*Mexico*

**Results Processing in MATLAB**

I.V. Guryev, I.A. Sukhoivanov, N.S. Gurieva, J.A. Andrade Lucio

*University of Guanajuato, Campus Irapuato-Salamanca, Division of Engineering*

The chapter is intended to provide the reader with powerful and flexible tools based on MATLAB and its open-source analogs for the processing and analyzing the results obtained

Particularly, in the chapter we give a brief overview of free and open-source as well as shareware software for computation of the photonic crystals characteristics. We concentrate attention mostly to their advantages and drawbacks and give short description of the data

The next parts of the chapter are dedicated to processing of the results obtained within the

Firstly, we consider the most basic case of plane results represented by a functional dependences. In this part, we are talking about the plane data interpolation, approximation

Data interpolation may be used as a tool intended to avoid high resolution computation (when scanning, say, through the wavelength or frequency) and, therefore, to reduce computation

Data approximation allows highly efficient data analysis by means of finding the parameters of the target function. For this reason, the numerically obtained data is approximated (within the certain accuracy) by the analytical functions including variable coefficients which are usually posses physical meaning. As a result of the approximation, the unknown coefficients

Correct and obvious data representation is very important when writing the scientific paper.

The next part of the chapter shows reader how to process and represent three-dimensional and multi-dimensional data. Basically, the operations described for plane data is extrapolated

In many cases, the results of the complex numerical research is represented in the form of multiple files with similar format. In the next part of the chapter we show how to use MATLAB for merging the data into a single multidimensional array and for representing it

time which may achieve weeks in the modern numerical problems.

Therefore, in the end of current part we give an idea of figure formatting.

**1. Introduction**

files they give results in.

specific software.

and representation.

to multidimensional data arrays.

in convenient form.

are found.

by means of highly specialized software.

**for Photonics Applications**

and O. Ibarra-Manzano


### **Results Processing in MATLAB for Photonics Applications**

I.V. Guryev, I.A. Sukhoivanov, N.S. Gurieva, J.A. Andrade Lucio and O. Ibarra-Manzano *University of Guanajuato, Campus Irapuato-Salamanca, Division of Engineering Mexico*

#### **1. Introduction**

22 Will-be-set-by-IN-TECH

118 MATLAB for Engineers – Applications in Control, Electrical Engineering, IT and Robotics

Cho, Y.-K., Song, T.-H. & Kim, H.-S. (2004). An optimized algorithm for computing the

Inacio, C. & Ombres, D. (1996). The DSP Decision fixed point or Floating?, *IEEE Spectrum*

ITRS (2009). International technology roadmap for semiconductors, 2009 edition

*Circuits and Systems*, USC Inf. Sci. Inst., Marina del Rey, CA, pp. 3331–3334. Malík, P. (2007). A generic IP Core of Indetical Forward and Inverse 12/36 Point MDCT

Malík, P., Ufnal, M., Łuczyk, A. W., Baláž, M. & Pleskacz, W. (2009). MDCT / IMDCT Low

Moore, G. E. (1965). Cramming more components onto integrated circuits, *Electronics*

Princen, J. P. & Bradley, A. B. (1986). Analysis/Synthesis Filter Bank Design Based on

Princen, J. P., Johnson, A. W. & Bradley, A. B. (1987). Sub-band/transform coding using

Spanias, A., Painer, T. & Atti, V. (2007). *Audio Signal Processing and Coding*, John Wiley & Sons,

Šimlaštík, M., Malík, P., Pikula, T. & Baláž, M. (2006). FPGA Implementation of a Fast MDCT

Vol. A, Chiang Mai, Thailand, pp. 626–628.

*Magazine of Technology Insiders* 33(9): 72–74.

*Scale Integration*, Nice, France, pp. 18–23.

Republic, pp. 144–147.

*Processing* 34(5): 1153–1161.

*ICASSP'87*, Dallas, TX, pp. 2161–2164.

38(8): 114–117.

Inc., Hoboken, NJ.

modified discrete cosine transform and its inverse transform, *Proc. of the IEEE Region 10 International Conference on Analog and Digital Techniques in Electrical Engineering*,

executive summary, *Technical report*. http://www.itrs.net/Links/2009ITRS/ 2009Chapters\_2009Tables/2009\_ExecSum.pdf [Accessed: April 2011]. Kwon, T. J., Sondeen, J. & Draper, J. (2005). Design Trade-Offs in Floating-Point Unit

Implementation for Embedded and Processing-In-Memory Systems, *Proc. of the IEEE*

Architecture and an Architectural Model Simulation Toolbox, *Proc. of. the 14th IEEE International Conference on Electronics, Cirsuits and Systems*, Marrakech, Maroko. Malík, P., Baláž, M., Pikula, T. & Šimlaštík, M. (2006). MDCT IP Core Generator with

Architectural Model Simulation, *Proc. of the IFIP International Conference on Very Large*

Power Implementations in 90 nm CMOS Technology for MP3 Audio, *Proc. of the IEEE Symposium on Design and Diagnostics of Electronic Circuits and Systems*, Liberec, Czech

Time Domain Aliasing Cancellation, *IEEE Transaction on Acoustic, Speech and Signal*

filter bank designs based on time-domain aliasing cancellation, *Proc. of the IEEE*

Algorithm, *Proc. Of the IEEE Workshop on Design and Diagnostics of Electronic Circuits and Systems*, Praha, Czech Technical University Publishing House, pp. 228–233.

The chapter is intended to provide the reader with powerful and flexible tools based on MATLAB and its open-source analogs for the processing and analyzing the results obtained by means of highly specialized software.

Particularly, in the chapter we give a brief overview of free and open-source as well as shareware software for computation of the photonic crystals characteristics. We concentrate attention mostly to their advantages and drawbacks and give short description of the data files they give results in.

The next parts of the chapter are dedicated to processing of the results obtained within the specific software.

Firstly, we consider the most basic case of plane results represented by a functional dependences. In this part, we are talking about the plane data interpolation, approximation and representation.

Data interpolation may be used as a tool intended to avoid high resolution computation (when scanning, say, through the wavelength or frequency) and, therefore, to reduce computation time which may achieve weeks in the modern numerical problems.

Data approximation allows highly efficient data analysis by means of finding the parameters of the target function. For this reason, the numerically obtained data is approximated (within the certain accuracy) by the analytical functions including variable coefficients which are usually posses physical meaning. As a result of the approximation, the unknown coefficients are found.

Correct and obvious data representation is very important when writing the scientific paper. Therefore, in the end of current part we give an idea of figure formatting.

The next part of the chapter shows reader how to process and represent three-dimensional and multi-dimensional data. Basically, the operations described for plane data is extrapolated to multidimensional data arrays.

In many cases, the results of the complex numerical research is represented in the form of multiple files with similar format. In the next part of the chapter we show how to use MATLAB for merging the data into a single multidimensional array and for representing it in convenient form.

for Photonics Applications 3

Results Processing in MATLAB for Photonics Applications 121

Within all the variety of the scientific software available, the developers are allowed to choose between low price and friendly interface. The friendly interface provides the user with a lot of default settings and saves a lot of time which may be spent to the research problems. On the other hand, being more complicated for user, free software requires deep knowledge of physical and mathematical aspects of the investigated area (which is considered by authors as an advantage of free software). Particularly, the scientist intended to investigate the optical properties of the periodic structures has to be well-educated in the fields of solid-state physics, electromagnetism, numerical methods, etc. Having these knowledge, the user may define the problem more precisely and carry out the research at higher level with deeper understanding

Therefore, this part of our review is dedicated to the software which is available online and

The MIT Photonic Bands (MPB) *MPB* (2011) software is one of the most powerful tools available freeware which is intended to compute the eigenstates of the photonic crystal of different configurations by means of PWE method. It is written by developers of the PhC theory whose investment into the physics of PhC devices was critical. Therefore, the MPB

However, since MBP runs script files with detailed description of the structure and computation conditions, using of the MBP requires from user high skills, particularly, in solid-state physics and turns computation of even the simplest case into a serious programming task. Several examples provided with the installation help to understand the principles of the structure setup and they may be used as a default structures which is easy to

The program outputs results into console which allows to read them but not analyze. Therefore, file output should be used. Particularly, the field distribution is being output into HDF5 format which is supported by visualization tools like h5utils. However, even for simple analysis, approximation and decoration the user has to use additional software such

Meep (or MEEP) *MEEP* (2011) is a free finite-difference time-domain (FDTD) simulation

It can be used to successfully compute electromagnetic field distribution inside the PhC-devices made either of linear or nonlinear materials (possessing Kerr of Pockels

In contrast to free software, the shareware ones have, in most cases, highly advanced user interface which helps novice to study the new soft and, on the other hand, saves professionals'

As for the PhC modeling, the most advanced computation tools available on the markets are RSoft (namely, BandSolve and FullWave) and Comsol Multiphysics (namely, RF module).

time keeping them from the necessity of learning side skills such as programming.

software package developed at MIT to model electromagnetic systems.

**2.1 Free and open-source software**

of the matter of the problem.

**2.1.1 MPB**

as MATLAB.

**2.1.2 MEEP**

nonlinearities).

**2.2 Shareware software**

may be easily tested and learned by reader.

modify and use for other structures.

solves the problems in a most suitable way for the scientists.

The last part of the chapter is dedicated to the animations creation out of multiple instant data shots. Such form of the data representation is useful when representing time-dependent results as well as multidimensional data arrays.

During all the chapter, we provide reader with working MATLAB codes as well as figures illustrating the programs results. We believe this information together with learning of MATLAB fundamentals *MATLAB manual* (2011) will help young scientists, master and PhD students especially in the area of opto-electronics and photonics to analyze and represent their computation results in the most effective way.

#### **2. Brief review of the Photonic crystals modeling software**

Photonic crystals are optical media possessing the periodic modulation of the refractive index. Due to this, the light behavior inside such structures is similar to the one of the electrons in the atomic structures (this gives them the name "Photonic crystals"). Main property the PhCs possesses is the photonic band gap (PBG) which is an optical analog of the electronic band gap in semiconductor materials.

Depending on the geometrical peculiarities, the PhCs are subdivided to several categories. Particularly, there are wide classes of 1D, 2D and 3D PhCs classified by number of dimensions where the variations of the refractive index appear. These classes, in turn, are divided according to the lattice type, existence or absence of the PBG, etc.

The physical principles of the PhCs are well-known and are considered in numerous articles and books Joannopoulos et al. (1995), Lourtioz et al. (2005), Sakoda (2001). There exist numerous methods for modeling the PhC characteristics, which are well described Sakoda (2001), Taflove & Hagness (2005) and are implemented in a numerous software products. This chapter is considered by us as a natural addition to the book Sukhoivanov & Guryev (2009) which describes computation methods of the PhCs characteristics. It is intended to extend and improve the efficiency of the results representation.

Particularly, this section is dedicated to brief introduction into freeware and shareware software for PhC and PhC-based devices basic modeling.

Since the PhC-based devices are investigated intensively today, the crucial moment is the right choice of the computational methods which are implemented in a number of software.

Particularly, the development of the PhC devices require detailed analysis of the PhC eigenstates (or resonant frequencies) as a dependence of the direction of the light propagation. This kind of tasks may be done by means of different methods like analytical one (however, for 1D PhC only), finite differences time domain (FDTD) method, finite elements method (FEM), plane wave expansion method (PWE). Each method has its advantages and drawbacks, they all require different computation time and accuracy and are implemented in a different kinds of software.

The other important problem for PhC investigation is the computation of the field distribution inside the defined structure. Again, in the simplest cases (particularly, in case of 1D PhC with passive optical materials) the field distribution can be found analytically. However, in the most cases the numerical methods are required such as FDTD, FEM, PWE, etc.

Most of the software (especially, free ones) are intended to implement all the aspects of a single methods. However, professional shareware soft implements several different methods allowing the user to investigate the device from different points of view.

#### **2.1 Free and open-source software**

Within all the variety of the scientific software available, the developers are allowed to choose between low price and friendly interface. The friendly interface provides the user with a lot of default settings and saves a lot of time which may be spent to the research problems. On the other hand, being more complicated for user, free software requires deep knowledge of physical and mathematical aspects of the investigated area (which is considered by authors as an advantage of free software). Particularly, the scientist intended to investigate the optical properties of the periodic structures has to be well-educated in the fields of solid-state physics, electromagnetism, numerical methods, etc. Having these knowledge, the user may define the problem more precisely and carry out the research at higher level with deeper understanding of the matter of the problem.

Therefore, this part of our review is dedicated to the software which is available online and may be easily tested and learned by reader.

#### **2.1.1 MPB**

2 MATLAB/Book4

The last part of the chapter is dedicated to the animations creation out of multiple instant data shots. Such form of the data representation is useful when representing time-dependent

During all the chapter, we provide reader with working MATLAB codes as well as figures illustrating the programs results. We believe this information together with learning of MATLAB fundamentals *MATLAB manual* (2011) will help young scientists, master and PhD students especially in the area of opto-electronics and photonics to analyze and represent their

Photonic crystals are optical media possessing the periodic modulation of the refractive index. Due to this, the light behavior inside such structures is similar to the one of the electrons in the atomic structures (this gives them the name "Photonic crystals"). Main property the PhCs possesses is the photonic band gap (PBG) which is an optical analog of the electronic band

Depending on the geometrical peculiarities, the PhCs are subdivided to several categories. Particularly, there are wide classes of 1D, 2D and 3D PhCs classified by number of dimensions where the variations of the refractive index appear. These classes, in turn, are divided

The physical principles of the PhCs are well-known and are considered in numerous articles and books Joannopoulos et al. (1995), Lourtioz et al. (2005), Sakoda (2001). There exist numerous methods for modeling the PhC characteristics, which are well described Sakoda (2001), Taflove & Hagness (2005) and are implemented in a numerous software products. This chapter is considered by us as a natural addition to the book Sukhoivanov & Guryev (2009) which describes computation methods of the PhCs characteristics. It is intended to extend

Particularly, this section is dedicated to brief introduction into freeware and shareware

Since the PhC-based devices are investigated intensively today, the crucial moment is the right choice of the computational methods which are implemented in a number of software. Particularly, the development of the PhC devices require detailed analysis of the PhC eigenstates (or resonant frequencies) as a dependence of the direction of the light propagation. This kind of tasks may be done by means of different methods like analytical one (however, for 1D PhC only), finite differences time domain (FDTD) method, finite elements method (FEM), plane wave expansion method (PWE). Each method has its advantages and drawbacks, they all require different computation time and accuracy and are implemented in a different kinds

The other important problem for PhC investigation is the computation of the field distribution inside the defined structure. Again, in the simplest cases (particularly, in case of 1D PhC with passive optical materials) the field distribution can be found analytically. However, in the

Most of the software (especially, free ones) are intended to implement all the aspects of a single methods. However, professional shareware soft implements several different methods

most cases the numerical methods are required such as FDTD, FEM, PWE, etc.

allowing the user to investigate the device from different points of view.

results as well as multidimensional data arrays.

computation results in the most effective way.

gap in semiconductor materials.

of software.

**2. Brief review of the Photonic crystals modeling software**

according to the lattice type, existence or absence of the PBG, etc.

and improve the efficiency of the results representation.

software for PhC and PhC-based devices basic modeling.

The MIT Photonic Bands (MPB) *MPB* (2011) software is one of the most powerful tools available freeware which is intended to compute the eigenstates of the photonic crystal of different configurations by means of PWE method. It is written by developers of the PhC theory whose investment into the physics of PhC devices was critical. Therefore, the MPB solves the problems in a most suitable way for the scientists.

However, since MBP runs script files with detailed description of the structure and computation conditions, using of the MBP requires from user high skills, particularly, in solid-state physics and turns computation of even the simplest case into a serious programming task. Several examples provided with the installation help to understand the principles of the structure setup and they may be used as a default structures which is easy to modify and use for other structures.

The program outputs results into console which allows to read them but not analyze. Therefore, file output should be used. Particularly, the field distribution is being output into HDF5 format which is supported by visualization tools like h5utils. However, even for simple analysis, approximation and decoration the user has to use additional software such as MATLAB.

#### **2.1.2 MEEP**

Meep (or MEEP) *MEEP* (2011) is a free finite-difference time-domain (FDTD) simulation software package developed at MIT to model electromagnetic systems.

It can be used to successfully compute electromagnetic field distribution inside the PhC-devices made either of linear or nonlinear materials (possessing Kerr of Pockels nonlinearities).

#### **2.2 Shareware software**

In contrast to free software, the shareware ones have, in most cases, highly advanced user interface which helps novice to study the new soft and, on the other hand, saves professionals' time keeping them from the necessity of learning side skills such as programming.

As for the PhC modeling, the most advanced computation tools available on the markets are RSoft (namely, BandSolve and FullWave) and Comsol Multiphysics (namely, RF module).

for Photonics Applications 5

Results Processing in MATLAB for Photonics Applications 123

As a result of computation of any PhC characteristic such as band structure, field distribution, temporal response, etc., the data is alway represented in some special format which can be

With simulation software you may obtain several different kinds of data files, namely, the plane data (two-coordinate dependence), 3D data (tree-coordinate dependence), random data (where data points are linked to the mesh nodes), structured data (the data is represented in the form of known structures), multiple files data (usually, the parametric data where each file

The data represented in a plane form usually represents the tabulated function of one variable and is really the simplest data format. In many cases, the computational programs output the

In the one wants to plot or process this kind of data, then it can be easily read by MATLAB

However, in some cases, the plane data file has the text header before the data itself (see figure 2 which demonstrates the data where the first line is a text header). This causes the problem since the data cannot be read with load command from the MATLAB. One of the solution is to use the Import Data dialog (Menu File/Import Data. . . ) which performs smart data analysis

On the other hand, it is more suitable to integrate the data import procedure into the script in order to avoid manual operations and, to be able to perform multiple import operations when processing multiple data files. For this reason, we have to use standard data reading operations to analyze file. The drawback of this method is necessity to know exactly the file structure, namely, how many strings takes the text header. Here is the simple example of

corresponds to the solution with a single value of the parameter).

data in a form of two rows as presented in figure 1

Fig. 1. The example of the plane data file

varName=load(fileName, ' ASCII')

After this, the variable varName contains an array Nx2.

and separates the numerical data from the text one.

opening the file generated by RSoft and having the text header.

script using a single command:

**3. Data files formats**

**3.1 Reading plane data**

read and processed with MATLAB.

#### **2.2.1 RSoft**

RSoft *RSoft* (2011) being developed by the RSoft-group for many years and it have come a long way since its earlier versions to its modern view.

Though its interface has been changed from version to version, it uses the same methods for computation of the PhC characteristics. Namely, BanSolve module allows one to compute the band structures of the PhC as well as its modes' field distribution by means of PWE method while FullWave implements FDTD method applying it to the field distribution computation inside the arbitrary structures and, in special case, the band structure.

The RSoft interface provides user with many useful models generated with default parameters. Moreover, wise linking to the variables allows flexible modifications of the default structures. This allows you to create the simple structure in a several mouse clicks and more advanced structures in several minutes. Among the drawbacks of the interface it should be mentioned too simple 3D structure definition which does not allow visual representation of the structure. To do this, additional software such as MayaVi should be used.

Except the possibility of modifying the structure with variables, RSoft's dialog boxes provide dozens of visual controls allowing definition of an arbitrary structure, computation conditions and the output information.

It should also be mentioned the graphic output of the RSoft. Although it allows analyze the computation results, their quality is unsatisfactory and one have to use additional software to make results representation suitable for publications. Moreover, in case of FullWave module, continuous graphic output essentially slows down the computation process.

In general, the RSoft is a powerful PhC modeling tool which allows solution of wide range of problems. However, results output is performed into the data files only and requires additional postprocessing, particularly, with Matlab.

#### **2.2.2 Comsol Multiphysics**

The Comsol Multiphysics *COMSOL* (2011) (previously known as Femlab) allows solving wide variety of problems defined by the partial differential equations in a number of fields of Physics by means of FEM. Particularly, in the area of PhC devices it allows to solve the Helmholtz equation with certain kinds of boundary conditions, thus, giving as a result the electromagnetic field distribution inside the device.

The structure definition as well as setting the properties of the structure, boundary conditions and node conditions are highly visualized and it makes no difficulties in creating computation model.

As for results representation, Comsol Multiphysics provides a lot of different visualization modes in form of 2D, 3D, contour, etc., thus, making for user unnecessary to deal with visualization of quite complicated data structures.

However, field distribution and basic data analysis is not always enough for serious scientific research. In case the PhC device has complex structure and requires detailed analysis of the field distribution (like confinement factor of the waveguide), the output files data should be investigated separately.

Therefore, independent on the convenience of the software itself, the problem of the output data processing is always arising and the best way to carry this out is by means of MATLAB.

#### **3. Data files formats**

4 MATLAB/Book4

RSoft *RSoft* (2011) being developed by the RSoft-group for many years and it have come a

Though its interface has been changed from version to version, it uses the same methods for computation of the PhC characteristics. Namely, BanSolve module allows one to compute the band structures of the PhC as well as its modes' field distribution by means of PWE method while FullWave implements FDTD method applying it to the field distribution computation

The RSoft interface provides user with many useful models generated with default parameters. Moreover, wise linking to the variables allows flexible modifications of the default structures. This allows you to create the simple structure in a several mouse clicks and more advanced structures in several minutes. Among the drawbacks of the interface it should be mentioned too simple 3D structure definition which does not allow visual representation

Except the possibility of modifying the structure with variables, RSoft's dialog boxes provide dozens of visual controls allowing definition of an arbitrary structure, computation conditions

It should also be mentioned the graphic output of the RSoft. Although it allows analyze the computation results, their quality is unsatisfactory and one have to use additional software to make results representation suitable for publications. Moreover, in case of FullWave module,

In general, the RSoft is a powerful PhC modeling tool which allows solution of wide range of problems. However, results output is performed into the data files only and requires

The Comsol Multiphysics *COMSOL* (2011) (previously known as Femlab) allows solving wide variety of problems defined by the partial differential equations in a number of fields of Physics by means of FEM. Particularly, in the area of PhC devices it allows to solve the Helmholtz equation with certain kinds of boundary conditions, thus, giving as a result the

The structure definition as well as setting the properties of the structure, boundary conditions and node conditions are highly visualized and it makes no difficulties in creating computation

As for results representation, Comsol Multiphysics provides a lot of different visualization modes in form of 2D, 3D, contour, etc., thus, making for user unnecessary to deal with

However, field distribution and basic data analysis is not always enough for serious scientific research. In case the PhC device has complex structure and requires detailed analysis of the field distribution (like confinement factor of the waveguide), the output files data should be

Therefore, independent on the convenience of the software itself, the problem of the output data processing is always arising and the best way to carry this out is by means of MATLAB.

long way since its earlier versions to its modern view.

additional postprocessing, particularly, with Matlab.

electromagnetic field distribution inside the device.

visualization of quite complicated data structures.

inside the arbitrary structures and, in special case, the band structure.

of the structure. To do this, additional software such as MayaVi should be used.

continuous graphic output essentially slows down the computation process.

**2.2.1 RSoft**

and the output information.

**2.2.2 Comsol Multiphysics**

investigated separately.

model.

As a result of computation of any PhC characteristic such as band structure, field distribution, temporal response, etc., the data is alway represented in some special format which can be read and processed with MATLAB.

With simulation software you may obtain several different kinds of data files, namely, the plane data (two-coordinate dependence), 3D data (tree-coordinate dependence), random data (where data points are linked to the mesh nodes), structured data (the data is represented in the form of known structures), multiple files data (usually, the parametric data where each file corresponds to the solution with a single value of the parameter).

#### **3.1 Reading plane data**

The data represented in a plane form usually represents the tabulated function of one variable and is really the simplest data format. In many cases, the computational programs output the data in a form of two rows as presented in figure 1


Fig. 1. The example of the plane data file

In the one wants to plot or process this kind of data, then it can be easily read by MATLAB script using a single command:

```
varName=load(fileName, ' ASCII')
```
After this, the variable varName contains an array Nx2.

However, in some cases, the plane data file has the text header before the data itself (see figure 2 which demonstrates the data where the first line is a text header). This causes the problem since the data cannot be read with load command from the MATLAB. One of the solution is to use the Import Data dialog (Menu File/Import Data. . . ) which performs smart data analysis and separates the numerical data from the text one.

On the other hand, it is more suitable to integrate the data import procedure into the script in order to avoid manual operations and, to be able to perform multiple import operations when processing multiple data files. For this reason, we have to use standard data reading operations to analyze file. The drawback of this method is necessity to know exactly the file structure, namely, how many strings takes the text header. Here is the simple example of opening the file generated by RSoft and having the text header.

for Photonics Applications 7

Results Processing in MATLAB for Photonics Applications 125

This code opens the data file as a regular file and then uses the function scanf() to read the data of the certain format. Particularly, for the file presented in the figure 2 it first reads 9 strings separated by the whitespaces and then reads each line containing 8 numbers into an

In case of plane data, each column has its meaning and have to be treated by the program accordingly. For instance, in case of the RSoft file generated by the scanning the full photonic band gaps (PBG) over some parameter, the first row indicates the value of the parameter, the second one is for the number of the PBG starting from the one with the lowest frequency, the third is for the central frequency of the PBG and so on. Detailed information on the specific

3D data is represented in the file in form of 2D array. In this case, each dimension of the array corresponds to the space dimension and, therefore, each number placed in the position (X,Y) in the file is treated as an altitude of the point with coordinates (X\*scale,Y\*scale). The value of

Usually, reading of the 3D data does not differ from reading the plane data. In case of files without headers, the data is read by a simple load() function of the MATLAB. The scale of the structure should be known in this case to be able to assign physical coordinates to the

However, in case of the header presence, we have to use the method described in the subsection 3.1 for plane data. In this case the header may contain useful information such

as linking to the physical dimensions, position of the structure, etc (see figure 3).

array 1x8. Then in the cycle it copies each row into the resulting 2D array.

files formats is usually contained in the manual for the specific software.

the variable 'scale' depends on the physical dimensions of the structure.

**end**

*%closing the data file*

**3.2 Reading 3D data**

positions in the array.

Fig. 3. The 3D data file with header

fclose(file);


Fig. 2. The example of the data file containing text header

*%Code for reading the data file containing text header*

```
%opening the file
file=fopen('scan_gaps_TE.scn');
```

```
%skipping the file header (this part should be modified according
%to the specific file format)
fscanf(file,'%s',9);
```
*%Reading the first line of the data. %variable 'read' contains the data from the read string, %variable 'num' contains the number of variables actually read. %the length of the array to read to should be selected according %to the dimensions of the data in the file* [read, num]=fscanf(file,'%f',[1,8]);

*%Resetting the counter to write into the data array* counter=1;

```
%main reading cycle. Executes while the data read from the file
%contains actual data (num>0)
```
**while**(num)

*%copying the data from temporary variable into the data %variable* data(counter,:)=read;

*%increasing the counter value* counter=counter+1;

*%reading the next line from the file* [read, num]=fscanf(file,'%f',[1,8]);

#### **end**

6 MATLAB/Book4

Fig. 2. The example of the data file containing text header

*%skipping the file header (this part should be modified according*

*%variable 'num' contains the number of variables actually read. %the length of the array to read to should be selected according*

*%main reading cycle. Executes while the data read from the file*

*%copying the data from temporary variable into the data*

*%Code for reading the data file containing text header*

*%variable 'read' contains the data from the read string,*

*%Resetting the counter to write into the data array*

file=fopen('scan\_gaps\_TE.scn');

*%opening the file*

counter=1;

**while**(num)

*%variable*

*%to the specific file format)* fscanf(file,'%s',9);

*%Reading the first line of the data.*

*%contains actual data (num*>*0)*

data(counter,:)=read;

counter=counter+1;

*%increasing the counter value*

*%reading the next line from the file* [read, num]=fscanf(file,'%f',[1,8]);

*%to the dimensions of the data in the file* [read, num]=fscanf(file,'%f',[1,8]);

*%closing the data file* fclose(file);

This code opens the data file as a regular file and then uses the function scanf() to read the data of the certain format. Particularly, for the file presented in the figure 2 it first reads 9 strings separated by the whitespaces and then reads each line containing 8 numbers into an array 1x8. Then in the cycle it copies each row into the resulting 2D array.

In case of plane data, each column has its meaning and have to be treated by the program accordingly. For instance, in case of the RSoft file generated by the scanning the full photonic band gaps (PBG) over some parameter, the first row indicates the value of the parameter, the second one is for the number of the PBG starting from the one with the lowest frequency, the third is for the central frequency of the PBG and so on. Detailed information on the specific files formats is usually contained in the manual for the specific software.

#### **3.2 Reading 3D data**

3D data is represented in the file in form of 2D array. In this case, each dimension of the array corresponds to the space dimension and, therefore, each number placed in the position (X,Y) in the file is treated as an altitude of the point with coordinates (X\*scale,Y\*scale). The value of the variable 'scale' depends on the physical dimensions of the structure.

Usually, reading of the 3D data does not differ from reading the plane data. In case of files without headers, the data is read by a simple load() function of the MATLAB. The scale of the structure should be known in this case to be able to assign physical coordinates to the positions in the array.

However, in case of the header presence, we have to use the method described in the subsection 3.1 for plane data. In this case the header may contain useful information such as linking to the physical dimensions, position of the structure, etc (see figure 3).


Fig. 3. The 3D data file with header

for Photonics Applications 9

Results Processing in MATLAB for Photonics Applications 127

Particularly, in MATLAB there is a set of functions providing correct work with the HDF5. The MATLAB command HDF5Read(...) reads the contents of the file into the structure from where it can be easily extracted by standard structure addressing command (in case of MATLAB it is the dot operator). Here is simple MATLAB program which reads results computed by MEEP

However, as it have been described in the previous section, the data in the files may be represented in the forms of an arbitrary structures which requires a lot of work to dig in and

In many cases, when you perform the parametric computations, the results have form of multiple files. Each file corresponds to a single value of the parameter and may have one

The main problem in this case is to gather the data from all the files into a single data structure. This problem will be discussed in the following section. Here we will focus on how to open multiple files. The following piece of code shows how to open the files which have similar

from HDF5 file and analyzes parameters of the results.

*%field distribution computed with MEEP*

*%required to read the actual data*

clear all

*%structure*

**for** i=1:(size(res))(1)

drawres(:,:)=res(i,:,:);

*%Drawing the frame* image(drawres∗1000); colormap('gray');

*%Drawing immediately*

drawnow;

**end**

*%Program allows reading from the the .h5 file containing the results of the*

*%The function hdf5info allows obtatining the information from the .h5 file*

*%Then we are reading the information using the query string from the hinfo*

*%In this particular case, the computation results contain several frames %with field distribution at different moments of time. The cycle reads them*

hinfo=hdf5info('/home/igor/temp/wg\_bend ez.h5');

res=hdf5read(hinfo.GroupHierarchy.Datasets(1));

*%from the data structure and then draws.*

to write corresponding unique decoders.

names with a counter attached to them:

**3.5 Reading multiple files data**

of the formats described earlier.

*%Copying the frame into a separate 2D array*

In this case, instead of just skipping the header block, we can extract useful information from it. The following code demonstrates how to read and analyze the header info when drawing the field distribution computed by RSoft:

*%Opening the data file* fff=fopen('field.fld');

*%Skipping two lines containing picture formatting info %(namely, two strings separated by the* <*CR*>*)* fscanf(fff,'%s',2);

*%Reading contents of the numerical part of the header into %variable 'field\_info'* field\_info=fscanf(fff,'%i %f %f %f OUTPUT\_REAL\_3D\n %i %f %f\n',7);

*%using the variable 'field\_info' to read properly formatted %data on the field distribution* field=fscanf(fff, '%f', [field\_info(5), field\_info(1)]);

Here, after opening the file, we are skipping two strings which mean nothing for us. Then we read all the numerical values into the array and can use this variable to determine the dimension of the array, the physical size of the computation region, etc.

#### **3.3 Reading the random data**

Using the header to define the parameters of the figure is not always enough. Particularly, it works for the uniformly distributed data points such as in FDTD method. However, if you use FEM, for example, Comsol Multiphysics, you will deal with non-uniform triangular mesh. In this case, the data is organized into structured. Each unit contains X and Y coordinate (in case of 2D model) or X, Y and Z coordinate (in case of 3D physical model) and the value of the field intensity corresponding to this point (see figure 4).


Fig. 4. The file containing random data linked to the specific points of the mesh

Reading of the date in this case is performed just as in case of the plane data, however, one have to fill two or three coordinate arrays X, Y and possibly Z, and the array of corresponding values of the field intensity.

It is also possible to convert this data into a regular mesh data which will be described in following section.

#### **3.4 Reading structured data**

In case of MPB and MEEP software, the field distribution is being output in so-called HDF5 format which is supported by many visualization tools as well as programming languages. Particularly, in MATLAB there is a set of functions providing correct work with the HDF5. The MATLAB command HDF5Read(...) reads the contents of the file into the structure from where it can be easily extracted by standard structure addressing command (in case of MATLAB it is the dot operator). Here is simple MATLAB program which reads results computed by MEEP from HDF5 file and analyzes parameters of the results.

*%Program allows reading from the the .h5 file containing the results of the %field distribution computed with MEEP*

clear all

8 MATLAB/Book4

In this case, instead of just skipping the header block, we can extract useful information from it. The following code demonstrates how to read and analyze the header info when drawing

the field distribution computed by RSoft:

*%Skipping two lines containing picture formatting info %(namely, two strings separated by the* <*CR*>*)*

*%Reading contents of the numerical part of the header into*

*%using the variable 'field\_info' to read properly formatted*

field=fscanf(fff, '%f', [field\_info(5), field\_info(1)]);

intensity corresponding to this point (see figure 4).

field\_info=fscanf(fff,'%i %f %f %f OUTPUT\_REAL\_3D\n %i %f %f\n',7);

dimension of the array, the physical size of the computation region, etc.

Fig. 4. The file containing random data linked to the specific points of the mesh

Here, after opening the file, we are skipping two strings which mean nothing for us. Then we read all the numerical values into the array and can use this variable to determine the

Using the header to define the parameters of the figure is not always enough. Particularly, it works for the uniformly distributed data points such as in FDTD method. However, if you use FEM, for example, Comsol Multiphysics, you will deal with non-uniform triangular mesh. In this case, the data is organized into structured. Each unit contains X and Y coordinate (in case of 2D model) or X, Y and Z coordinate (in case of 3D physical model) and the value of the field

Reading of the date in this case is performed just as in case of the plane data, however, one have to fill two or three coordinate arrays X, Y and possibly Z, and the array of corresponding

It is also possible to convert this data into a regular mesh data which will be described in

In case of MPB and MEEP software, the field distribution is being output in so-called HDF5 format which is supported by many visualization tools as well as programming languages.

*%Opening the data file* fff=fopen('field.fld');

fscanf(fff,'%s',2);

*%variable 'field\_info'*

*%data on the field distribution*

**3.3 Reading the random data**

values of the field intensity.

**3.4 Reading structured data**

following section.

*%The function hdf5info allows obtatining the information from the .h5 file %required to read the actual data* hinfo=hdf5info('/home/igor/temp/wg\_bend ez.h5');

*%Then we are reading the information using the query string from the hinfo %structure* res=hdf5read(hinfo.GroupHierarchy.Datasets(1));

*%In this particular case, the computation results contain several frames %with field distribution at different moments of time. The cycle reads them %from the data structure and then draws.*

```
for i=1:(size(res))(1)
```

```
%Copying the frame into a separate 2D array
drawres(:,:)=res(i,:,:);
```
*%Drawing the frame* image(drawres∗1000); colormap('gray');

*%Drawing immediately* drawnow;

**end**

However, as it have been described in the previous section, the data in the files may be represented in the forms of an arbitrary structures which requires a lot of work to dig in and to write corresponding unique decoders.

#### **3.5 Reading multiple files data**

In many cases, when you perform the parametric computations, the results have form of multiple files. Each file corresponds to a single value of the parameter and may have one of the formats described earlier.

The main problem in this case is to gather the data from all the files into a single data structure. This problem will be discussed in the following section. Here we will focus on how to open multiple files. The following piece of code shows how to open the files which have similar names with a counter attached to them:

for Photonics Applications 11

Results Processing in MATLAB for Photonics Applications 129

The interpolation is the most useful operation on the data obtained during the experimental research on numerical computation. In the mathematical subfield of numerical analysis, interpolation is a method of constructing new data points within the range of a discrete set of known data points. There are a lot of different ways to interpolate the data like linear interpolation, polynomial interpolation, spline interpolation, etc. We will not consider mathematical side of the interpolation process concentrating on the implementation of this

For the interpolation of the plane data, MATLAB uses the function interp1(). As an input arguments, the function has initial X and Y values of the data as well as new values of X at which you want to obtain the new values. It is also possible to select the interpolation method. As an output, this function returns the new values of Y found by interpolation process. Here is the simple example of the data interpolation. In this example, we have the dependence of the photonic band gap width of the 1D PhC on the layer width. The data points are calculated for the limited number of values of the witdh. However, due to interpolation,

**4.1 Data interpolation**

operation in MATLAB.

*%layers witdh*

x=data(:,1)';

*%in the datafile* y=data(:,5)';

*%for interpolation*

hold on;

*%Creating new figurefigure; %Drawing at the same figure*

y1=interp1(x,y,x1,'nearest');

plot(x,y, 'bo', 'LineWidth',2);

plot(x1,y1, 'g', 'LineWidth',2);

*%Doing the same with linear interpolation*

we have smooth curve which looks much nicer.

*%The program demonstrates different interpolation techniques %performed by MATLAB. The interpolation is applied*

*%Loading the information on the PBG width from the file*

*%Getting the values of Y for the second band (third column)*

*%Creating the mesh 10 times more dense which is required*

*%Using the function interp1, interpolating the data with nearest values*

data=load('pbg\_width.dat',' ASCII');

*%Getting the values of X from the data array*

x1=interp1(1:length(x),x,1:0.1:length(x));

*%Plotting the initial data with blue round markers*

*%Plotting the interpolated data with grin solid line*

*%to the dependence of the photonic band gap width of 1D PhC on the*

*%Scanning over all the numbers of the files* **for** num=0:90

name=strcat(path,'metadata',num2str(num,'%02i'),'.dat')

clear temporal\_response; temporal\_response=load(name, ' ASCII');

*%Here performing all the operations we need* ... ...

**end**

We presume to have file names in the following form: "metadata00.dat", "metadata01.dat", "metadata02.dat",. . . ,"metadata90.dat". In this code, the outer loop scans over all the numbers from 0 to 90. The function strcat() allows concatenating several strings into one and generate the name. Using the function num2str() we convert the numerical type into a string.

The problem, may appear with the filenames since for the numbers lower than 10 we have zero-filler in the name. We can fix the problem in two different ways. The first one is just to compare the variable 'num' with 10 and if it's lower than 10, just to add '0' to the string manually. However, this is a solution only in particular cases. Suppose, you have numbers from 0 to 100000. In this case, you will have to make 5 different comparisons each for one digit (*num* < 10, *num* < 100, *num* < 1000, *num* < 10000, *num* < 100000). This is already a sort of annoying and far from elegance.

The second possible solution is to use formatted conversion of the number. In this case we just have to use format string "%02i" to fill the spaces before the number with zeros. Here the symbol of "%" tells MATLAB about the beginning of the format string, following "0" makes MATLAB to fill the spaces with zeros, "2" tells exactly how many digits it is suppose to be in the number and "i" indicates integer data type. Detailed information about the string formatting can be found in the MATLAB documentation in the section of the function fprintf().

#### **4. Processing plane data**

The plane numerical data is represented by two arrays, one of which contains values along one coordinate axis while another one corresponds to another axis. This simplest form of the data representation is actually the most clear and understandable when you want to present your results.

However, in some cases direct plotting of the data is not really representable. Particularly, when you make an expensive experimental research like temporal response measurement of the nonlinear PhC using ultra-short pulses, or when you perform the computation which is so time- and resource-consuming that you are unable to obtain sufficient number of data points like, for instance, computation of the photonic band gap maps. In these cases, before to create the data output, you have to convert it, improve it or make some other handling. This section is dedicated to the manipulations on the data which allow to improve its representation.

#### **4.1 Data interpolation**

10 MATLAB/Book4

We presume to have file names in the following form: "metadata00.dat", "metadata01.dat", "metadata02.dat",. . . ,"metadata90.dat". In this code, the outer loop scans over all the numbers from 0 to 90. The function strcat() allows concatenating several strings into one and generate the name. Using the function num2str() we convert the numerical type into a

The problem, may appear with the filenames since for the numbers lower than 10 we have zero-filler in the name. We can fix the problem in two different ways. The first one is just to compare the variable 'num' with 10 and if it's lower than 10, just to add '0' to the string manually. However, this is a solution only in particular cases. Suppose, you have numbers from 0 to 100000. In this case, you will have to make 5 different comparisons each for one digit (*num* < 10, *num* < 100, *num* < 1000, *num* < 10000, *num* < 100000). This is already a sort of

The second possible solution is to use formatted conversion of the number. In this case we just have to use format string "%02i" to fill the spaces before the number with zeros. Here the symbol of "%" tells MATLAB about the beginning of the format string, following "0" makes MATLAB to fill the spaces with zeros, "2" tells exactly how many digits it is suppose to be in the number and "i" indicates integer data type. Detailed information about the string formatting can be found in the MATLAB documentation in the section of the function fprintf().

The plane numerical data is represented by two arrays, one of which contains values along one coordinate axis while another one corresponds to another axis. This simplest form of the data representation is actually the most clear and understandable when you want to present

However, in some cases direct plotting of the data is not really representable. Particularly, when you make an expensive experimental research like temporal response measurement of the nonlinear PhC using ultra-short pulses, or when you perform the computation which is so time- and resource-consuming that you are unable to obtain sufficient number of data points like, for instance, computation of the photonic band gap maps. In these cases, before to create the data output, you have to convert it, improve it or make some other handling. This section is dedicated to the manipulations on the data which allow to improve its representation.

name=strcat(path,'metadata',num2str(num,'%02i'),'.dat')

*%Scanning over all the numbers of the files*

temporal\_response=load(name, ' ASCII');

*%Here performing all the operations we need*

clear temporal\_response;

annoying and far from elegance.

**4. Processing plane data**

your results.

**for** num=0:90

... ...

**end**

string.

The interpolation is the most useful operation on the data obtained during the experimental research on numerical computation. In the mathematical subfield of numerical analysis, interpolation is a method of constructing new data points within the range of a discrete set of known data points. There are a lot of different ways to interpolate the data like linear interpolation, polynomial interpolation, spline interpolation, etc. We will not consider mathematical side of the interpolation process concentrating on the implementation of this operation in MATLAB.

For the interpolation of the plane data, MATLAB uses the function interp1(). As an input arguments, the function has initial X and Y values of the data as well as new values of X at which you want to obtain the new values. It is also possible to select the interpolation method. As an output, this function returns the new values of Y found by interpolation process.

Here is the simple example of the data interpolation. In this example, we have the dependence of the photonic band gap width of the 1D PhC on the layer width. The data points are calculated for the limited number of values of the witdh. However, due to interpolation, we have smooth curve which looks much nicer.

*%The program demonstrates different interpolation techniques %performed by MATLAB. The interpolation is applied %to the dependence of the photonic band gap width of 1D PhC on the %layers witdh*

*%Loading the information on the PBG width from the file* data=load('pbg\_width.dat',' ASCII');

*%Getting the values of X from the data array* x=data(:,1)';

*%Getting the values of Y for the second band (third column) %in the datafile* y=data(:,5)';

*%Creating the mesh 10 times more dense which is required %for interpolation* x1=interp1(1:length(x),x,1:0.1:length(x));

*%Creating new figurefigure; %Drawing at the same figure* hold on; *%Using the function interp1, interpolating the data with nearest values* y1=interp1(x,y,x1,'nearest'); *%Plotting the initial data with blue round markers* plot(x,y, 'bo', 'LineWidth',2); *%Plotting the interpolated data with grin solid line* plot(x1,y1, 'g', 'LineWidth',2);

*%Doing the same with linear interpolation*

for Photonics Applications 13

Results Processing in MATLAB for Photonics Applications 131

In the program, the interpolation is performed in four different ways provided by MATLAB (actually, there are also "pchip" and "v5cubic" types. However, pchip is similar to cubic and v5cubic is just an old version of cubic). From the results, it is seen that for the selected curve

Data approximation is an operation which allows to describe (approximate) the tabulated data

The approximation is always used when the general view of the analytical dependence can be predicted, however, the coefficients still remain unknown. This is usually happens in the experimental researches when it is necessary to determine the parameters of the analytical dependence. Particularly, when the experiment is intent to determine the effective refractive index of the PhC by measuring the radiation deviation at different wavelength. After the experiment, the data is approximated by parametric analytical function where the parameter is the effective refractive index. The best approximation result gives the value of

In MATLAB there is a separate toolbox providing the approximation of the tabulated data, namely, the Curve Fitting Toolbox. And again, using of the toolbox is possible in a single case only. However, using the toolbox requires definite time and when you need to perform

In this section we describe how to approximate the data in MATLAB without using the toolbox. The following program generates randomly five points and then approximates them

the best interpolation technique is a spline interpolation.

**4.2 Data approximation**

by an analytical function.

the mentioned parameter.

*%the error should be absent.*

*%Creating five X values*

y=rand(1,5)∗10;

start=[1 1 1 1 1];

*%points (X, Y)*

start=estimated\_coefs;

**for** i=1:5

*%over*

x=1:5;

**function** approximation\_plane

*%Randomly generating 5 Y values*

*%Starting five minimization interations*

multiple operations, this turns into a mess.

with polynomial function of the fourth order.

*%Creating the initial values for approximation process*

*%This function approximates five random points with 4 th order*

*%polynomial function. Since 4 th order polynom allows exact approximation,*

*%Using the function fminsearch() to optimize the 'coefs' a the selected*

estimated\_coefs = (fminsearch(@(coefs)fitfun(coefs,x,y),start));

*%changing the starting values to the optimized ones and then starting*

figure; hold on y1=interp1(x,y,x1,'linear'); plot(x,y, 'bo', 'LineWidth',2); plot(x1,y1, 'g', 'LineWidth',2);

*%Doing the same with spline interpolation* figure; hold on; y1=interp1(x,y,x1,'spline'); plot(x,y, 'bo', 'LineWidth',2); plot(x1,y1, 'g', 'LineWidth',2);

```
%Doing the same with cubic interpolation
figure;
hold on;
y1=interp1(x,y,x1,'cubic');
plot(x,y, 'bo', 'LineWidth',2);
plot(x1,y1, 'g', 'LineWidth',2);
```
As a result of the program we obtain four figures each for one interpolation type.

Fig. 5. Examples of basic interpolation. a) Nearest, b) linear, c)spline and d) cubic

In the program, the interpolation is performed in four different ways provided by MATLAB (actually, there are also "pchip" and "v5cubic" types. However, pchip is similar to cubic and v5cubic is just an old version of cubic). From the results, it is seen that for the selected curve the best interpolation technique is a spline interpolation.

#### **4.2 Data approximation**

12 MATLAB/Book4

As a result of the program we obtain four figures each for one interpolation type.

x 10−6

x 10−6

Fig. 5. Examples of basic interpolation. a) Nearest, b) linear, c)spline and d) cubic

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

(d)

(b)

x 10−6

x 10−6

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

(c)

(a)

figure; hold on

figure; hold on;

figure; hold on;

y1=interp1(x,y,x1,'linear'); plot(x,y, 'bo', 'LineWidth',2); plot(x1,y1, 'g', 'LineWidth',2);

y1=interp1(x,y,x1,'spline'); plot(x,y, 'bo', 'LineWidth',2); plot(x1,y1, 'g', 'LineWidth',2);

y1=interp1(x,y,x1,'cubic'); plot(x,y, 'bo', 'LineWidth',2); plot(x1,y1, 'g', 'LineWidth',2);

*%Doing the same with spline interpolation*

*%Doing the same with cubic interpolation*

Data approximation is an operation which allows to describe (approximate) the tabulated data by an analytical function.

The approximation is always used when the general view of the analytical dependence can be predicted, however, the coefficients still remain unknown. This is usually happens in the experimental researches when it is necessary to determine the parameters of the analytical dependence. Particularly, when the experiment is intent to determine the effective refractive index of the PhC by measuring the radiation deviation at different wavelength. After the experiment, the data is approximated by parametric analytical function where the parameter is the effective refractive index. The best approximation result gives the value of the mentioned parameter.

In MATLAB there is a separate toolbox providing the approximation of the tabulated data, namely, the Curve Fitting Toolbox. And again, using of the toolbox is possible in a single case only. However, using the toolbox requires definite time and when you need to perform multiple operations, this turns into a mess.

In this section we describe how to approximate the data in MATLAB without using the toolbox. The following program generates randomly five points and then approximates them with polynomial function of the fourth order.

*%This function approximates five random points with 4 th order %polynomial function. Since 4 th order polynom allows exact approximation, %the error should be absent.*

**function** approximation\_plane *%Creating five X values* x=1:5;

*%Randomly generating 5 Y values* y=rand(1,5)∗10;

*%Creating the initial values for approximation process* start=[1 1 1 1 1];

*%Starting five minimization interations*

**for** i=1:5

*%Using the function fminsearch() to optimize the 'coefs' a the selected %points (X, Y)*

estimated\_coefs = (fminsearch(@(coefs)fitfun(coefs,x,y),start));

*%changing the starting values to the optimized ones and then starting %over* start=estimated\_coefs;

for Photonics Applications 15

Results Processing in MATLAB for Photonics Applications 133

+(76.4662)⋅x<sup>2</sup>

+(−130.3002)⋅x+(1.3917

3

+(−17.599)⋅x

<sup>0</sup> <sup>1</sup> <sup>2</sup> <sup>3</sup> <sup>4</sup> <sup>5</sup> <sup>6</sup> −10

As an example, we give the code which approximates the dependence of central frequency of the photonic band gap of 1D PhC on the layer width, with hyperbolic function and estimates the approximation error. The data is obtained by means of plane wave expansion method Sukhoivanov & Guryev (2009) with low accuracy. With approximation it is possible

Fig. 6. Result of approximation of the data by the polynomial function

to determine the analytic function describing such a dependence.

*%This function approximates the dependence of the central %wavelength of the photonic band gap of 1D PhC on the layer*

**function** approximation\_plane\_err\_estimation *%Loading information about the central frequency of the*

data=load('pbg\_center.dat',' ASCII');

*%Getting the values of Y for the second band (third column)*

*%Now there's only four values since we use 3 rd order polynom*

*%Using the function fminsearch() to optimize the 'coefs' a the selected*

estimated\_coefs = (fminsearch(@(coefs)fitfun(coefs,x,y),start));

*%Creating the initial values for approximation process*

*%Getting the values of X from the data array*

Consider the following code:

*%width with hyperbolic function. %It also performs the error estimation.*

*%bandgap of 1D PhC from the file*

x=data(:,1)';

*%in the datafile* y=data(:,3)';

start=[1e14 0];

*%points (X, Y)*

**end**

```
%Building the smooth analytical function with coefficients obtained during
%optimization process
x1=0:0.1:6;
y1=analyt_fn(estimated_coefs,x1);
```

```
%Creating the figure
figure;
hold on;
%Plotting the initial points
plot(x,y,'bo','LineWidth',2);
%and the smooth analytical function
plot(x1,y1,'g','LineWidth',2);
```

```
%In the title writing the final polynom
title(strcat('y=',num2str(estimated_coefs(1)),'\cdotxˆ4+(',...
  num2str(estimated_coefs(2)),')\cdotxˆ3+(',...
  num2str(estimated_coefs(3)),')\cdotxˆ2+(',...
  num2str(estimated_coefs(4)),')\cdotx+(',...
  num2str(estimated_coefs(1)),')'),'FontSize',16);
```

```
%Function for optimization
```

```
function err=fitfun(coefs,x,y)
  %The function returns the difference between the analytical function
  %with coefficients 'coefs' at points X and the actual values of Y
  err=norm(analyt_fn(coefs,x) y);
```
**end**

```
%Analytical parametric function to be optimized
```

```
function analyt_fn=analyt_fn(coefs,x)
  analyt_fn=coefs(1)∗x.ˆ4+coefs(2)∗x.ˆ3+coefs(3)∗x.ˆ2+coefs(4)∗x+coefs(5);
end
```
**end**

The result of the function work is presented in figure 6. It can bee seen that approximation is quite accurate.

For the optimization procedure, MATLAB uses simplex search method which is implemented in the form of function fminsearch(). As an input parameters it has a pointer to a function to be minimized and the initial values of the parameter. In our case, the function to minimize returns the difference between actual points and the analytical function which, in turn, is defined in the function analyt\_fn().

It should be paid attention to the fact that in current example the polynomial function of the fourth power approximates five data points. In this case, there is only one strict solution. However, in most cases, the approximation is made with certain accuracy which should be properly estimated. The best way to estimate the accuracy is to measure minimum, maximum and avarage deviation.

14 MATLAB/Book4

*%Building the smooth analytical function with coefficients obtained during*

title(strcat('y=',num2str(estimated\_coefs(1)),'\cdotxˆ4+(',...

*%The function returns the difference between the analytical function %with coefficients 'coefs' at points X and the actual values of Y*

analyt\_fn=coefs(1)∗x.ˆ4+coefs(2)∗x.ˆ3+coefs(3)∗x.ˆ2+coefs(4)∗x+coefs(5);

The result of the function work is presented in figure 6. It can bee seen that approximation is

For the optimization procedure, MATLAB uses simplex search method which is implemented in the form of function fminsearch(). As an input parameters it has a pointer to a function to be minimized and the initial values of the parameter. In our case, the function to minimize returns the difference between actual points and the analytical function which, in turn, is

It should be paid attention to the fact that in current example the polynomial function of the fourth power approximates five data points. In this case, there is only one strict solution. However, in most cases, the approximation is made with certain accuracy which should be properly estimated. The best way to estimate the accuracy is to measure minimum, maximum

num2str(estimated\_coefs(2)),')\cdotxˆ3+(',... num2str(estimated\_coefs(3)),')\cdotxˆ2+(',... num2str(estimated\_coefs(4)),')\cdotx+(',... num2str(estimated\_coefs(1)),')'),'FontSize',16);

**end**

x1=0:0.1:6;

figure; hold on;

**end**

**end end**

quite accurate.

*%optimization process*

*%Creating the figure*

*%Plotting the initial points* plot(x,y,'bo','LineWidth',2); *%and the smooth analytical function* plot(x1,y1,'g','LineWidth',2);

*%Function for optimization* **function** err=fitfun(coefs,x,y)

y1=analyt\_fn(estimated\_coefs,x1);

*%In the title writing the final polynom*

err=norm(analyt\_fn(coefs,x) y);

defined in the function analyt\_fn().

and avarage deviation.

*%Analytical parametric function to be optimized* **function** analyt\_fn=analyt\_fn(coefs,x)

Fig. 6. Result of approximation of the data by the polynomial function

As an example, we give the code which approximates the dependence of central frequency of the photonic band gap of 1D PhC on the layer width, with hyperbolic function and estimates the approximation error. The data is obtained by means of plane wave expansion method Sukhoivanov & Guryev (2009) with low accuracy. With approximation it is possible to determine the analytic function describing such a dependence. Consider the following code:

*%This function approximates the dependence of the central %wavelength of the photonic band gap of 1D PhC on the layer %width with hyperbolic function. %It also performs the error estimation.*

**function** approximation\_plane\_err\_estimation *%Loading information about the central frequency of the %bandgap of 1D PhC from the file* data=load('pbg\_center.dat',' ASCII');

*%Getting the values of X from the data array* x=data(:,1)';

*%Getting the values of Y for the second band (third column) %in the datafile* y=data(:,3)';

*%Creating the initial values for approximation process %Now there's only four values since we use 3 rd order polynom* start=[1e14 0];

*%Using the function fminsearch() to optimize the 'coefs' a the selected %points (X, Y)* estimated\_coefs = (fminsearch(@(coefs)fitfun(coefs,x,y),start));

for Photonics Applications 17

Results Processing in MATLAB for Photonics Applications 135

0 0.2 0.4 0.6 0.8 1

Now the result is represented in the form shown in the figure 7. We currently used the hyperbolic function which gave good agreement with numerical results. In this case, not all the points (in general case, neither of them) are fall at the analytical curve. Thus, estimating the error, we are able to confirm the validity of the approximation. In case of large error value, the function is unable to approximate the data at any conditions. This may happen when you have selected wrong kind of approximating function or, in case of too complex function, have

In the scientific publication the data representation plays main role when it comes to correct understanding of the computation or experimental results. In two previous subsections we have described the way for data interpolation and approximation which helps to perform the data analysis as well as smooth the data. However, the data presentation is also important. Some basic techniques have already been used earlier without explanation such as manipulation of the figure's and plot's parameters. In this subsection we will try to describe

First thing which should be known about MATLAB, it is the object-oriented environment. This means that after creating some object such as plot, text, line, etc., they can be modified

As always, the most simple way to manipulate the objects parameters is to use visual environment provided by MATLAB. For this reason, it is enough just to double-click at some object while editing mode is active (for this, you have to press the mouse pointer icon at the

y=9.36e+08/(x+4.93e−07)

Minimum error =6.10e+10 Maximum error =4.70e+13 Avarage error =1.45e+13

x 10−6

*%The function returns the difference between the analytical function %with coefficients 'coefs' at points X and the actual values of Y*

> 0.6 0.8 1 1.2 1.4 1.6 1.8 <sup>2</sup> x 1015

Fig. 7. Displaying the approximation accuracy information

selected wrong values of the initial parameters.

**4.3 Data representation**

and explain the most useful of them.

and their properties can be easily changed.

err=norm(analyt\_fn(coefs,x) y);

*%Analytical parametric function to be optimized* **function** analyt\_fn=analyt\_fn(coefs,x) analyt\_fn=coefs(1)./(x+coefs(2));

**end**

**end end**

*%Building the smooth analytical function with coefficients obtained during %optimization process* x1=interp1(1:length(data(:,1)),data(:,1),1:0.1:length(data(:,1))); y1=analyt\_fn(estimated\_coefs,x1);

*%using the function 'find' searching for indices of the refined X values %which are equal to the initial ones* [row est\_ind]=find(meshgrid(x1,x) meshgrid(x,x1)' == 0);

*%Searching for minimum deviation of the analytical function from the initial %data (using est\_index to compare corresponding points)* minerr=min(abs(y1(est\_ind) y));

```
%In the same manner, searching for maximum deviation
maxerr=max(abs(y1(est_ind) y));
```

```
%Calculating avarage error
avgerr=sum(abs(y1(est_ind) y))/length(y)
```

```
%Creating the figure
figure;
hold on;
```
*%plotting the solution* plot(x,y,'bo','LineWidth',2); *%and the smooth analytical function* plot(x1,y1,'g','LineWidth',2);

```
%Writing the error values at the figure
text(min(x1),min(y1)+(max(y1) min(y1))/10,...
  strcat('Minimum error = ', num2str(minerr,'%.2e')),'FontSize', 14);
text(min(x1),min(y1)+(max(y1) min(y1))/10∗2,...
  strcat('Maximum error = ', num2str(maxerr,'%.2e')),'FontSize', 14);
text(min(x1),min(y1)+(max(y1) min(y1))/10∗3,...
  strcat('Avarage error = ', num2str(avgerr,'%.2e')),'FontSize', 14);
```
*%Creating the title to the figure*

title(strcat('y=',num2str(estimated\_coefs(1),'%.2e'),'/(x+',... num2str(estimated\_coefs(2),'%.2e'),')'),'FontSize',16);

*%Function for optimization* **function** err=fitfun(coefs,x,y) 16 MATLAB/Book4

*%Building the smooth analytical function with coefficients obtained during*

x1=interp1(1:length(data(:,1)),data(:,1),1:0.1:length(data(:,1)));

*%using the function 'find' searching for indices of the refined X values*

*%Searching for minimum deviation of the analytical function from the initial*

[row est\_ind]=find(meshgrid(x1,x) meshgrid(x,x1)' == 0);

*%data (using est\_index to compare corresponding points)*

*%In the same manner, searching for maximum deviation*

*%optimization process*

y1=analyt\_fn(estimated\_coefs,x1);

*%which are equal to the initial ones*

minerr=min(abs(y1(est\_ind) y));

maxerr=max(abs(y1(est\_ind) y));

avgerr=sum(abs(y1(est\_ind) y))/length(y)

*%Calculating avarage error*

*%Creating the figure*

*%plotting the solution*

plot(x,y,'bo','LineWidth',2); *%and the smooth analytical function* plot(x1,y1,'g','LineWidth',2);

*%Writing the error values at the figure*

*%Creating the title to the figure*

*%Function for optimization* **function** err=fitfun(coefs,x,y)

text(min(x1),min(y1)+(max(y1) min(y1))/10,...

text(min(x1),min(y1)+(max(y1) min(y1))/10∗2,...

text(min(x1),min(y1)+(max(y1) min(y1))/10∗3,...

strcat('Minimum error = ', num2str(minerr,'%.2e')),'FontSize', 14);

strcat('Maximum error = ', num2str(maxerr,'%.2e')),'FontSize', 14);

strcat('Avarage error = ', num2str(avgerr,'%.2e')),'FontSize', 14);

title(strcat('y=',num2str(estimated\_coefs(1),'%.2e'),'/(x+',... num2str(estimated\_coefs(2),'%.2e'),')'),'FontSize',16);

figure; hold on;

*%The function returns the difference between the analytical function %with coefficients 'coefs' at points X and the actual values of Y* err=norm(analyt\_fn(coefs,x) y); **end**

```
%Analytical parametric function to be optimized
function analyt_fn=analyt_fn(coefs,x)
  analyt_fn=coefs(1)./(x+coefs(2));
end
end
```
Fig. 7. Displaying the approximation accuracy information

Now the result is represented in the form shown in the figure 7. We currently used the hyperbolic function which gave good agreement with numerical results. In this case, not all the points (in general case, neither of them) are fall at the analytical curve. Thus, estimating the error, we are able to confirm the validity of the approximation. In case of large error value, the function is unable to approximate the data at any conditions. This may happen when you have selected wrong kind of approximating function or, in case of too complex function, have selected wrong values of the initial parameters.

#### **4.3 Data representation**

In the scientific publication the data representation plays main role when it comes to correct understanding of the computation or experimental results. In two previous subsections we have described the way for data interpolation and approximation which helps to perform the data analysis as well as smooth the data. However, the data presentation is also important. Some basic techniques have already been used earlier without explanation such as manipulation of the figure's and plot's parameters. In this subsection we will try to describe and explain the most useful of them.

First thing which should be known about MATLAB, it is the object-oriented environment. This means that after creating some object such as plot, text, line, etc., they can be modified and their properties can be easily changed.

As always, the most simple way to manipulate the objects parameters is to use visual environment provided by MATLAB. For this reason, it is enough just to double-click at some object while editing mode is active (for this, you have to press the mouse pointer icon at the

for Photonics Applications 19

Results Processing in MATLAB for Photonics Applications 137

Here we show more advanced plotting from the approximation example. Pay attention that after plotting the data (X-values and Y-values), we place a lot of textual and numerical parameters. Parameters are usually defined by a pair name-value. The value can be a string,

Particularly, the first plot function will draw the data with the LineWidth=2 (initially its 0.5), without the lines connecting the data points (LineStyle=none), but with circular markers

The second plot draws the approximated data with LineWidth=2, with dashed line (LineStyle='–') and with color defined by three RGB values (each value varies from 0 to 1).

As we have mentioned previously, there is another way to change the parameters of the objects using the object identifier. Note that each function which creates an objects, also returns an object identifier which can be stored into a variable and used lately to manipulate the object. Particularly, the properties of the object can be easily changed using the function set() which has as a parameters the object identifier and the set of object's properties. Consider this simple example which makes the same as previous one but modifying the properties after the object

set(p1, 'LineWidth', 2, 'LineStyle', 'none', 'Marker', 'o', 'Color', 'k');

set(p2,'LineWidth',2, 'LineStyle', ' ', 'Color', [0 0.3 1]);

set(t1, 'Position',[min(x1)+0.1,min(y1)+(max(y1) min(y1))/10,0],...

plot(x,y, 'LineWidth', 2, 'LineStyle', 'none', 'Marker', 'o', 'Color', 'k');

plot(x1,y1,'LineWidth',2, 'LineStyle', ' ', 'Color', [0 0.3 1]);

strcat('Minimum error = ', num2str(minerr)),'FontSize', 14);

strcat('Maximum error = ', num2str(maxerr)),'FontSize', 14);

strcat('Avarage error = ', num2str(avgerr)),'FontSize', 14);

a number or an array depending on the parameter's requirements.

text(min(x1)+0.1,min(y1)+(max(y1) min(y1))/10,...

text(min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗2,...

text(min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗3,...

(Marker=o) drawn with black color (Color=k).

The text marks have only the parameters Fontsize=14.

*%Plotting the initial data*

*%Adding the text marks*

is created:

... ... ...

*%Approximation process*

*%Plotting the initial data*

*%Adding the text marks*

*%Plotting the approximated data*

p1=plot(x,y);

p2=plot(x1,y1);

t1=text();

*%Plotting the approximated data*

Fig. 8. Visual editing of the object's properties

tool panel). After this, you will see the editing panel which allows editing some parameters of the selected object (see figure 8(a)).

Deeper understanding of the parameters of the objects can be achieved when opening additional parameters panel or properties inspector (see figure 8(b)). Here you can see and edit all the parameters available for the object.

However, editing parameters inside your script is the most fast and productive way to change the figure appearance (imagine that you have to redraw the figure over and over and modify dozens of parameters manually).

There are two ways to set up the parameters of the objects. The first one is to set them on the stage of their creation. And the second one is modify them after the object is created with default properties.

In the first method, the parameters are transfered as arguments to the function which creates an objects like in following code:

*%Approximation process*

...

...

...

*%Plotting the initial data* plot(x,y, 'LineWidth', 2, 'LineStyle', 'none', 'Marker', 'o', 'Color', 'k'); *%Plotting the approximated data* plot(x1,y1,'LineWidth',2, 'LineStyle', ' ', 'Color', [0 0.3 1]);

```
%Adding the text marks
text(min(x1)+0.1,min(y1)+(max(y1) min(y1))/10,...
 strcat('Minimum error = ', num2str(minerr)),'FontSize', 14);
text(min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗2,...
 strcat('Maximum error = ', num2str(maxerr)),'FontSize', 14);
text(min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗3,...
 strcat('Avarage error = ', num2str(avgerr)),'FontSize', 14);
```
Here we show more advanced plotting from the approximation example. Pay attention that after plotting the data (X-values and Y-values), we place a lot of textual and numerical parameters. Parameters are usually defined by a pair name-value. The value can be a string, a number or an array depending on the parameter's requirements.

Particularly, the first plot function will draw the data with the LineWidth=2 (initially its 0.5), without the lines connecting the data points (LineStyle=none), but with circular markers (Marker=o) drawn with black color (Color=k).

The second plot draws the approximated data with LineWidth=2, with dashed line (LineStyle='–') and with color defined by three RGB values (each value varies from 0 to 1). The text marks have only the parameters Fontsize=14.

As we have mentioned previously, there is another way to change the parameters of the objects using the object identifier. Note that each function which creates an objects, also returns an object identifier which can be stored into a variable and used lately to manipulate the object. Particularly, the properties of the object can be easily changed using the function set() which has as a parameters the object identifier and the set of object's properties. Consider this simple example which makes the same as previous one but modifying the properties after the object is created:

*%Approximation process*

... ...

18 MATLAB/Book4

(a) (b)

tool panel). After this, you will see the editing panel which allows editing some parameters

Deeper understanding of the parameters of the objects can be achieved when opening additional parameters panel or properties inspector (see figure 8(b)). Here you can see and

However, editing parameters inside your script is the most fast and productive way to change the figure appearance (imagine that you have to redraw the figure over and over and modify

There are two ways to set up the parameters of the objects. The first one is to set them on the stage of their creation. And the second one is modify them after the object is created with

In the first method, the parameters are transfered as arguments to the function which creates

Fig. 8. Visual editing of the object's properties

edit all the parameters available for the object.

of the selected object (see figure 8(a)).

dozens of parameters manually).

an objects like in following code:

default properties.

... ... ...

*%Approximation process*

...

```
%Plotting the initial data
p1=plot(x,y);
set(p1, 'LineWidth', 2, 'LineStyle', 'none', 'Marker', 'o', 'Color', 'k');
```
*%Plotting the approximated data* p2=plot(x1,y1); set(p2,'LineWidth',2, 'LineStyle', ' ', 'Color', [0 0.3 1]);

```
%Adding the text marks
t1=text();
set(t1, 'Position',[min(x1)+0.1,min(y1)+(max(y1) min(y1))/10,0],...
```
for Photonics Applications 21

Results Processing in MATLAB for Photonics Applications 139

*%The variable 'name' contains filename without the extension*

*%Into the variable 'field\_info' reading useful information*

field=fscanf(fff, '%f', [field\_info(5), field\_info(1)]);

field\_info(3)∗1e 6,... *%to maximum X*

field\_info(6)∗1e 6:... *%from minimum Y*

field\_info(7)∗1e 6,... *%to maximum Y*

field, 'LineStyle', 'none');

*%Making MATLAB to draw at the same plot*

*%Opening the file with the structure boundaries*

fff=fopen(strcat(name,'.pdo'));

*%Setting the color map* colormap('copper');

hold on;

*%Reading the data points into 2D array. The dimensions of the*

*%Again, the information about absolute coordinates is taken from the*

contourf(field\_info(2)∗1e 6:... *%from minimum X*

*%with the step computed from the number of points over X* ((field\_info(3)∗1e 6) field\_info(2)∗1e 6)/(field\_info(1) 1):

*%with the step computed from the number of points over Y*

*%plotting the field data points without contour lines*

((field\_info(7)∗1e 6) field\_info(6)∗1e 6)/(field\_info(5) 1):...

field\_info=fscanf(fff,'%i %f %f %f OUTPUT\_REAL\_3D\n %i %f %f\n',7);

name='field\_high\_state'

*%which is unuseful for us* fscanf(fff,'%s',2);

figure;

*%variable 'field\_info'*

*%Opening the file containing the field info* fff=fopen(strcat(name,'.fld'));

*%Skipping the inforation about the plotting*

*%array are taken from the variable 'field\_info'*

*%Plotting the field with filled contour plot*

```
'String',strcat('Minimum error = ', num2str(minerr)),...
    'FontSize', 14);
t2=text();
set(t2, 'Position', [min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗2,0],...
    'String',strcat('Maximum error = ', num2str(maxerr)),...
    'FontSize', 14);
t3=text();
set(t3, 'Position', [min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗3,0],...
    'String',strcat('Maximum error = ', num2str(avgerr)),...
    'FontSize', 14);
```
This code does exactly the same as previous one, however with significant difference. When the function creating the object is called, the object identifier is stored into a corresponding variable. After this, the function set() is called to modify the objects properties. Pay attention, in case of text labels, the object is created with default properties. All the properties including the text strings are set later using the identifiers.

This code seems to be unnecessarily more awkward then the previous one but actually, it has one great advantage. After you stored the variable identifier, you are able to manipulate the object in any way you want. Particularly, you can easily create an animation by varying the parameters of the object (such as camera view) in a loop and saving each frame (for information about animations see the last section).

#### **5. Simple 3D data processing**

Operating with 3D data is more complicated than with plane one. Particularly, the most simple data to be work with is data organized into a regular 2D mesh such as results of FDTD computation. However, in some cases it is important to perform the data interpolation.

#### **5.1 Data on rectangular mesh**

The most frequently, the data organized into a rectangular mesh may be found in the field distribution output of the software implementing the FDTD method. However, each of the programs creates its own data format and, therefore, should be treated individually. Here as an example, we will give the programs reading the data from RSoft and MEEP.

#### **5.1.1 Processing the data from RSoft**

When the RSoft outputs the field distribution, it usually created two different files. The first one contains the field distribution and the second one contains the information about the structure boundaries. Therefore, it is necessary to read both files for correct representation the computed data.

Consider the following code for reading the RSoft field output:

*%The program draws the field distribution computed by RSoft %Over the field distribution it draws the structure which is %read from the file generated by RSoft*

clear all

*%The variable 'name' contains filename without the extension* name='field\_high\_state'

*%Opening the file containing the field info* fff=fopen(strcat(name,'.fld'));

20 MATLAB/Book4

This code does exactly the same as previous one, however with significant difference. When the function creating the object is called, the object identifier is stored into a corresponding variable. After this, the function set() is called to modify the objects properties. Pay attention, in case of text labels, the object is created with default properties. All the properties including

This code seems to be unnecessarily more awkward then the previous one but actually, it has one great advantage. After you stored the variable identifier, you are able to manipulate the object in any way you want. Particularly, you can easily create an animation by varying the parameters of the object (such as camera view) in a loop and saving each frame (for

Operating with 3D data is more complicated than with plane one. Particularly, the most simple data to be work with is data organized into a regular 2D mesh such as results of FDTD computation. However, in some cases it is important to perform the data interpolation.

The most frequently, the data organized into a rectangular mesh may be found in the field distribution output of the software implementing the FDTD method. However, each of the programs creates its own data format and, therefore, should be treated individually. Here as

When the RSoft outputs the field distribution, it usually created two different files. The first one contains the field distribution and the second one contains the information about the structure boundaries. Therefore, it is necessary to read both files for correct representation the

an example, we will give the programs reading the data from RSoft and MEEP.

Consider the following code for reading the RSoft field output:

*%The program draws the field distribution computed by RSoft %Over the field distribution it draws the structure which is*

'String',strcat('Minimum error = ', num2str(minerr)),...

set(t2, 'Position', [min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗2,0],... 'String',strcat('Maximum error = ', num2str(maxerr)),...

set(t3, 'Position', [min(x1)+0.1,min(y1)+(max(y1) min(y1))/10∗3,0],... 'String',strcat('Maximum error = ', num2str(avgerr)),...

'FontSize', 14);

'FontSize', 14);

'FontSize', 14);

the text strings are set later using the identifiers.

information about animations see the last section).

**5. Simple 3D data processing**

**5.1 Data on rectangular mesh**

computed data.

clear all

**5.1.1 Processing the data from RSoft**

*%read from the file generated by RSoft*

t2=text();

t3=text();

*%Skipping the inforation about the plotting %which is unuseful for us* fscanf(fff,'%s',2);

*%Into the variable 'field\_info' reading useful information* field\_info=fscanf(fff,'%i %f %f %f OUTPUT\_REAL\_3D\n %i %f %f\n',7);

*%Reading the data points into 2D array. The dimensions of the %array are taken from the variable 'field\_info'* field=fscanf(fff, '%f', [field\_info(5), field\_info(1)]);

*%Plotting the field with filled contour plot* figure;

*%Again, the information about absolute coordinates is taken from the %variable 'field\_info'* contourf(field\_info(2)∗1e 6:... *%from minimum X*

*%with the step computed from the number of points over X* ((field\_info(3)∗1e 6) field\_info(2)∗1e 6)/(field\_info(1) 1):

field\_info(3)∗1e 6,... *%to maximum X*

field\_info(6)∗1e 6:... *%from minimum Y*

*%with the step computed from the number of points over Y* ((field\_info(7)∗1e 6) field\_info(6)∗1e 6)/(field\_info(5) 1):...

field\_info(7)∗1e 6,... *%to maximum Y*

*%plotting the field data points without contour lines* field, 'LineStyle', 'none');

*%Setting the color map* colormap('copper');

*%Making MATLAB to draw at the same plot* hold on;

*%Opening the file with the structure boundaries* fff=fopen(strcat(name,'.pdo'));

for Photonics Applications 23

Results Processing in MATLAB for Photonics Applications 141

In several cases, such as computation of the electromagnetic field distribution inside the PhC by means of FEM, the results appear within the nodes of irregular mesh. In this case it is impossible to perform the interpolation or other regular operation by means of standard MATLAB functions. In a professional software such as Comsol Multiphysics, this problem is overcome and now it is possible to export the results linked to a regular or user-defined

However, for example, in case of PDE Toolbox available in MATLAB, the results are exported in a form of random data where the specific value correspond to a node of triangular mesh.

In many cases, it is not necessary to output the data but we need to obtain some information from it. For instance, when one investigates the transmission characteristics of the microstructure, it is necessary to analyze the field distribution before the structure and after the structure. Namely, one knows exactly the regions of space where the field intensity should be found. In case of triangular mesh, the simplest way to do this is to find all the points within each region, the intensity in each point and then divide it by the number of points in the region. Consider the following code which finds average intensity within the region in

**6. Random data processing**

Dealing with this kind of data becomes a total mess.

space using the results exported from Comsol Multiphysics:

*%Function for computation of the transmittance of the nonlinear %photonic crystal from the field distribution computed in Comsol*

*%Loading the file with random data containing the field distribution*

field=load('results\_tri\_mesh.txt','ASCII');

*%Determining maximum and minimum values of coordinates %Considering y coordinate only since we are interested in*

*%Defining the limits for optical intensity calculation. Namely,*

*%Counter of nodes in the bottom region (corresponding to reflection)*

**6.1 Processing data on triangular mesh**

*%inside the investigated structure*

*%the whole width of the structure*

ymax=max(field(:,2)); ymin=min(field(:,2));

*%from the bottom* y00=ymin;

y01=ymin+3;

y10=ymax 3; *%to the top* y11=ymax;

*%to bottom+3 microns*

*%from the top 3 microns*

mesh.

clear all

```
%Reading header with useful information about filling
[isfill count]=fscanf(fff,'fill %i\n');
[color count] =fscanf(fff,'color %i\n');
num_points=fscanf(fff,'polygon %i\n',1);
while(num_points)
```

```
str=(fscanf(fff,'%f ',[2 num_points]))'∗1e 6;
fscanf(fff, 'end polygon\n',1);
patch(str(:,1),str(:,2),zeros(num_points,1),'FaceColor',[1 1 1],...
        'EdgeColor',[1 1 1],'LineWidth',2);
num_points=fscanf(fff,'polygon %i\n',1)
```
**end**

fclose(fff);

In this case, the first file "field\_high\_state.fld" is a simple 3D data with header containing information linking to the field distribution to real coordinate system, some auxiliary information about the number of data points, etc. Before writing the program reading and displaying such kind of data, it is useful to consult with the software manual and, which is very important, to open the file with text editor and check the header and data formats manually.

The second file "field\_high\_state.pdo" which contains information about the optical structure, contains structured data and should be read accordingly. Namely, each of the boundaries of the structure is stored in the form of polygon. The information about the polygon size is contained in the polygon's header. Before working with this file, it is also useful to take a look with a simple text editor.

#### **5.2 3D data visualization**

Unlike 2D data containing only two independent coordinates, the 3D data may be visualized in a variety of different ways. Depending on convenience, the following visualizing functions are most frequently used:


All the function actually take the same parameters, namely, X, Y arrays containing coordinates and Z 2D-array containing actual data.

Such kinds of visualization possess different functionality and speed which condition their usage. Particularly, if you need to express the intensity levels, the "contour()" function is the right choice. However, this function is too slow and trying to make animation with it requires a lot of time. The best choice for the animation is the function "image()".

#### **6. Random data processing**

22 MATLAB/Book4

patch(str(:,1),str(:,2),zeros(num\_points,1),'FaceColor',[1 1 1],...

In this case, the first file "field\_high\_state.fld" is a simple 3D data with header containing information linking to the field distribution to real coordinate system, some auxiliary information about the number of data points, etc. Before writing the program reading and displaying such kind of data, it is useful to consult with the software manual and, which is very important, to open the file with text editor and check the header and data formats

The second file "field\_high\_state.pdo" which contains information about the optical structure, contains structured data and should be read accordingly. Namely, each of the boundaries of the structure is stored in the form of polygon. The information about the polygon size is contained in the polygon's header. Before working with this file, it is also useful to take a look

Unlike 2D data containing only two independent coordinates, the 3D data may be visualized in a variety of different ways. Depending on convenience, the following visualizing functions

All the function actually take the same parameters, namely, X, Y arrays containing coordinates

Such kinds of visualization possess different functionality and speed which condition their usage. Particularly, if you need to express the intensity levels, the "contour()" function is the right choice. However, this function is too slow and trying to make animation with it requires

a lot of time. The best choice for the animation is the function "image()".

'EdgeColor',[1 1 1],'LineWidth',2);

*%Reading header with useful information about filling*

str=(fscanf(fff,'%f ',[2 num\_points]))'∗1e 6;

num\_points=fscanf(fff,'polygon %i\n',1)

fscanf(fff, 'end polygon\n',1);

[isfill count]=fscanf(fff,'fill %i\n'); [color count] =fscanf(fff,'color %i\n'); num\_points=fscanf(fff,'polygon %i\n',1);

**while**(num\_points)

**end** fclose(fff);

manually.

• surf() • mesh() • contour() • contourf() • image()

with a simple text editor.

**5.2 3D data visualization**

are most frequently used:

and Z 2D-array containing actual data.

In several cases, such as computation of the electromagnetic field distribution inside the PhC by means of FEM, the results appear within the nodes of irregular mesh. In this case it is impossible to perform the interpolation or other regular operation by means of standard MATLAB functions. In a professional software such as Comsol Multiphysics, this problem is overcome and now it is possible to export the results linked to a regular or user-defined mesh.

However, for example, in case of PDE Toolbox available in MATLAB, the results are exported in a form of random data where the specific value correspond to a node of triangular mesh. Dealing with this kind of data becomes a total mess.

#### **6.1 Processing data on triangular mesh**

In many cases, it is not necessary to output the data but we need to obtain some information from it. For instance, when one investigates the transmission characteristics of the microstructure, it is necessary to analyze the field distribution before the structure and after the structure. Namely, one knows exactly the regions of space where the field intensity should be found. In case of triangular mesh, the simplest way to do this is to find all the points within each region, the intensity in each point and then divide it by the number of points in the region. Consider the following code which finds average intensity within the region in space using the results exported from Comsol Multiphysics:

*%Function for computation of the transmittance of the nonlinear %photonic crystal from the field distribution computed in Comsol*

clear all

*%Loading the file with random data containing the field distribution %inside the investigated structure* field=load('results\_tri\_mesh.txt','ASCII');

*%Determining maximum and minimum values of coordinates %Considering y coordinate only since we are interested in %the whole width of the structure* ymax=max(field(:,2)); ymin=min(field(:,2));

*%Defining the limits for optical intensity calculation. Namely, %from the bottom* y00=ymin; *%to bottom+3 microns* y01=ymin+3; *%from the top 3 microns* y10=ymax 3; *%to the top* y11=ymax;

*%Counter of nodes in the bottom region (corresponding to reflection)*

for Photonics Applications 25

Results Processing in MATLAB for Photonics Applications 143

xlabel('Radiation intensity', 'FontSize', 14, 'FontWeight', 'bold');

The program, in general, is quite simple. The triangular mesh is written into a file as a number of point coordinates with corresponding field intensity. It is just necessary to find all the points

To improve the data quality and perform regional analysis of the field distribution, it is more convenient to use the regular data array instead of random data. Therefore, in this section, we provide reader with useful program able to convert the triangular mesh into a regular mesh

In PDE Toolbox, the data is always exported as one-dimensional or two-dimensional array (depending on the specific task). In case of two-dimensional array, each column corresponds to a single solution. It is also can be exported the mesh information. Particularly, we are interested in nodes coordinates and also particular interest represents the information about nodes which form the triangles. By default, this information is contained in the variables named "u" (the solution), "p" (coordinates of the nodes) and "t" (information about triangles formation) respectively. The information about the triangles is contained in the array "t" in a

After exporting this information, we can now start converting the mesh using the following

ylabel('Transmittance', 'FontSize', 14, 'FontWeight', 'bold');

falling at the region required and sum all the intensities within the region.

*%points*

**end**

figure;

grid on;

code:

clear data

reflect=reflect\_total/count\_reflect; transmit=transmit\_total/count\_transmit;

intensities(count\_intensity)=intensity; count\_intensity=count\_intensity+1;

*%Plotting the result and formatting the plot*

set(gca, 'FontSize', 12)

**6.2 Converting into a regular mesh**

by performing simple linear 3D interpolation.

*%The program converts the data on triangular mesh*

form of indices of the coordinate points stored in the array "p".

*%obtained in PDE Tool of the MATLAB into a data on square mesh %The data from PDE Tool, namely, array of data points 'u', %Data of the triangles of the mesh 't' and the array of*

*%coordinates 'p' should be previously exported to MATLAB workspace*

plot(intensities,reflectance, 'LineWidth', 2);

*%Assuming the total field intensity as a sum of reflectance %and transmittance, calculating the reflectance of the structure* reflectance(count\_intensity)=reflect/(reflect+transmit);

count\_reflect=1;

*%Counter of nodes in the top region (corresponding to transmission)* count\_transmit=1;

```
%Determinig the length of the field structure
len=length(field);
```

```
%Counter for different intensity values
count_intensity=1;
```
*%The cycle for the intensity*

```
for intensity=0:0.01:0.5
  %The variable to sum the field intensity in all the points of
  %the bottom region
  reflect_total=0;
```

```
%The variable to sum the field intensity in all the points of
%the top region
transmit_total=0;
```

```
%The cycle for all the points in the mesh
for iii=1:len
```

```
%If the point is inside the bottom region,
if (field(iii,2)>y00)&&(field(iii,2)<y01)
 %and if it is not NaN value
 if(abs(field(iii,3))>=0)
    %Adding its absolute value to the total reflected field
    reflect_total=reflect_total+(abs(field(iii,2+count_intensity)));
    count_reflect=count_reflect+1;
 end
```

```
end
```

```
%If the point is inside the bottom region,
 if (field(iii,2)>y10)&&(field(iii,2)<y11)
   %and if it is not NaN value
   if(abs(field(iii,3))>=0)
      %Adding its absolute value to the total transmitted field
      transmit_total=transmit_total+(abs(field(iii,2+count_intensity)));
      count_transmit=count_transmit+1;
   end
 end
end
```
*%Normalizing the field intensities by the number of*

*%points* reflect=reflect\_total/count\_reflect; transmit=transmit\_total/count\_transmit;

*%Assuming the total field intensity as a sum of reflectance %and transmittance, calculating the reflectance of the structure* reflectance(count\_intensity)=reflect/(reflect+transmit); intensities(count\_intensity)=intensity; count\_intensity=count\_intensity+1;

#### **end**

24 MATLAB/Book4

*%Counter of nodes in the top region (corresponding to transmission)*

*%The variable to sum the field intensity in all the points of*

*%The variable to sum the field intensity in all the points of*

*%Adding its absolute value to the total reflected field*

*%Adding its absolute value to the total transmitted field*

transmit\_total=transmit\_total+(abs(field(iii,2+count\_intensity)));

reflect\_total=reflect\_total+(abs(field(iii,2+count\_intensity)));

*%Determinig the length of the field structure*

*%The cycle for all the points in the mesh*

*%and if it is not NaN value* **if**(abs(field(iii,3))>=0)

*%If the point is inside the bottom region,* **if** (field(iii,2)>y00)&&(field(iii,2)<y01)

count\_reflect=count\_reflect+1;

*%If the point is inside the bottom region,* **if** (field(iii,2)>y10)&&(field(iii,2)<y11)

count\_transmit=count\_transmit+1;

*%Normalizing the field intensities by the number of*

*%and if it is not NaN value* **if**(abs(field(iii,3))>=0)

*%Counter for different intensity values*

count\_reflect=1;

count\_transmit=1;

len=length(field);

count\_intensity=1;

*%The cycle for the intensity* **for** intensity=0:0.01:0.5

*%the bottom region* reflect\_total=0;

*%the top region* transmit\_total=0;

**for** iii=1:len

**end end**

**end end end**

```
%Plotting the result and formatting the plot
figure;
plot(intensities,reflectance, 'LineWidth', 2);
xlabel('Radiation intensity', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Transmittance', 'FontSize', 14, 'FontWeight', 'bold');
set(gca, 'FontSize', 12)
grid on;
```
The program, in general, is quite simple. The triangular mesh is written into a file as a number of point coordinates with corresponding field intensity. It is just necessary to find all the points falling at the region required and sum all the intensities within the region.

#### **6.2 Converting into a regular mesh**

To improve the data quality and perform regional analysis of the field distribution, it is more convenient to use the regular data array instead of random data. Therefore, in this section, we provide reader with useful program able to convert the triangular mesh into a regular mesh by performing simple linear 3D interpolation.

In PDE Toolbox, the data is always exported as one-dimensional or two-dimensional array (depending on the specific task). In case of two-dimensional array, each column corresponds to a single solution. It is also can be exported the mesh information. Particularly, we are interested in nodes coordinates and also particular interest represents the information about nodes which form the triangles. By default, this information is contained in the variables named "u" (the solution), "p" (coordinates of the nodes) and "t" (information about triangles formation) respectively. The information about the triangles is contained in the array "t" in a form of indices of the coordinate points stored in the array "p".

After exporting this information, we can now start converting the mesh using the following code:

*%The program converts the data on triangular mesh %obtained in PDE Tool of the MATLAB into a data on square mesh %The data from PDE Tool, namely, array of data points 'u', %Data of the triangles of the mesh 't' and the array of %coordinates 'p' should be previously exported to MATLAB workspace*

clear data

for Photonics Applications 27

Results Processing in MATLAB for Photonics Applications 145

*%Then checking if any of the triangles contains the point of the*

*%For this reason, calculating the area of the triangle* area\_mesh=abs(tri\_area([p(1,t(1,triangles(tc)))... p(2,t(1,triangles(tc)))], ... [p(1,t(2,triangles(tc)))... p(2,t(2,triangles(tc)))],... [p(1,t(3,triangles(tc)))... p(2,t(3,triangles(tc)))]));

*%and the sum of areas of three triangles formed by the %point of the square mesh and each two points of the*

area\_sum=abs(tri\_area([xnew(xc) ynew(yc)],... [p(1,t(2,triangles(tc)))... p(2,t(2,triangles(tc)))],... [p(1,t(3,triangles(tc)))... p(2,t(3,triangles(tc)))]))+... abs(tri\_area([p(1,t(1,triangles(tc)))... p(2,t(1,triangles(tc)))],... [xnew(xc) ynew(yc)],... [p(1,t(3,triangles(tc)))... p(2,t(3,triangles(tc)))]))+... abs(tri\_area([p(1,t(1,triangles(tc)))... p(2,t(1,triangles(tc)))],... [p(1,t(2,triangles(tc)))... p(2,t(2,triangles(tc)))],... [xnew(xc) ynew(yc)])); *%The points is inside the triangular element if the sum of %the areas of three triangles is less or equivalent to the*

num\_t= 1;

*%rectangular mesh*

**for** tc=1:length(triangles)

*%triangular mesh element*

*%one of the triangluar element*

num\_t=triangles(tc); *%and finising the cycle*

break; **end end**

continue;

*%the compoutation area* **if**(num\_t== 1)

datanew(yc,xc)=NaN;

**if**(abs(area\_mesh area\_sum)<=0) *%writing the number of the triangle*

*%if the trianle not found then the requested point is outside*

clear datanew

*%For convenience, we will use the 2D array to store the coordinates %and corresponding data values %First writing the coordinates of the points* data=p';

*%and then, the values computed by PDE Toolbox in these points* data(:,3)=u(:,3)';

*%Finding minimum and maximum values of the coorditates* xmin=min(data(:,1)); xmax=max(data(:,1));

ymin=min(data(:,2)); ymax=max(data(:,2));

*%Defining the number of square mesh nodes in both dimensions %Pay attention that the best way is to sellect the number of %points in correspondence with the size on the structure in %each dimension* xnpoints=50; ynpoints=50;

```
%Generating coordinated of the rectangular mesh
xnew=xmin:(xmax xmin)/xnpoints:xmax;
ynew=ymin:(ymax ymin)/ynpoints:ymax;
```

```
%Now we have to find the data value in each of the points of the mesh
for xc=1:length(xnew)
```

```
for yc=1:length(ynew)
```

```
%To determine the point of the triangular mesh which is the closest
        %one to the current point of rectangular mesh, we compute the
        %Cartesian distances between them.
len=(sqrt((data(:,1) xnew(xc)).ˆ2+(data(:,2) ynew(yc)).ˆ2));
```

```
%Using the fucntion find(), we are searching for the minimum distance.
i1=find(len==min(len), 1, 'first');
if(len(i1)==0)
  datanew(yc,xc)=data(i1,3);
  continue;
```

```
end
```
*%Now we are searching for all the triangles containing the point found*

```
[row, triangles]=find(t==i1);
```

```
num_t= 1;
```
26 MATLAB/Book4

*%For convenience, we will use the 2D array to store the coordinates*

*%and then, the values computed by PDE Toolbox in these points*

*%Finding minimum and maximum values of the coorditates*

*%Defining the number of square mesh nodes in both dimensions %Pay attention that the best way is to sellect the number of %points in correspondence with the size on the structure in*

*%Now we have to find the data value in each of the points of the mesh*

*%Cartesian distances between them.*

i1=find(len==min(len), 1, 'first');

datanew(yc,xc)=data(i1,3);

[row, triangles]=find(t==i1);

len=(sqrt((data(:,1) xnew(xc)).ˆ2+(data(:,2) ynew(yc)).ˆ2));

*%To determine the point of the triangular mesh which is the closest %one to the current point of rectangular mesh, we compute the*

*%Using the fucntion find(), we are searching for the minimum distance.*

*%Now we are searching for all the triangles containing the point found*

*%Generating coordinated of the rectangular mesh* xnew=xmin:(xmax xmin)/xnpoints:xmax; ynew=ymin:(ymax ymin)/ynpoints:ymax;

clear datanew

data(:,3)=u(:,3)';

xmin=min(data(:,1)); xmax=max(data(:,1));

ymin=min(data(:,2)); ymax=max(data(:,2));

*%each dimension* xnpoints=50; ynpoints=50;

**for** xc=1:length(xnew) **for** yc=1:length(ynew)

**if**(len(i1)==0)

continue;

**end**

data=p';

*%and corresponding data values*

*%First writing the coordinates of the points*

```
%Then checking if any of the triangles contains the point of the
%rectangular mesh
  for tc=1:length(triangles)
     %For this reason, calculating the area of the triangle
     area_mesh=abs(tri_area([p(1,t(1,triangles(tc)))...
                    p(2,t(1,triangles(tc)))], ...
                   [p(1,t(2,triangles(tc)))...
                    p(2,t(2,triangles(tc)))],...
                   [p(1,t(3,triangles(tc)))...
                    p(2,t(3,triangles(tc)))]));
     %and the sum of areas of three triangles formed by the
     %point of the square mesh and each two points of the
     %triangular mesh element
     area_sum=abs(tri_area([xnew(xc) ynew(yc)],...
                   [p(1,t(2,triangles(tc)))...
                    p(2,t(2,triangles(tc)))],...
                   [p(1,t(3,triangles(tc)))...
                    p(2,t(3,triangles(tc)))]))+...
          abs(tri_area([p(1,t(1,triangles(tc)))...
                    p(2,t(1,triangles(tc)))],...
                   [xnew(xc) ynew(yc)],...
                   [p(1,t(3,triangles(tc)))...
                   p(2,t(3,triangles(tc)))]))+...
          abs(tri_area([p(1,t(1,triangles(tc)))...
                   p(2,t(1,triangles(tc)))],...
                   [p(1,t(2,triangles(tc)))...
                   p(2,t(2,triangles(tc)))],...
                   [xnew(xc) ynew(yc)]));
 %The points is inside the triangular element if the sum of
 %the areas of three triangles is less or equivalent to the
 %one of the triangluar element
     if(abs(area_mesh area_sum)
                                     <=0)
       %writing the number of the triangle
       num_t=triangles(tc);
       %and finising the cycle
       break;
     end
  end
  %if the trianle not found then the requested point is outside
  %the compoutation area
  if(num_t== 1)
     datanew(yc,xc)=NaN;
     continue;
```
for Photonics Applications 29

Results Processing in MATLAB for Photonics Applications 147

(a)

(b)

it finds the closest point, it uses the array "t" to find all the triangles containing this point. Then it checks which triangle contains current point of the rectangular mesh. For this reason it uses the method of the areas summation. If the one is not found, the point is outside the computation area and NaN is written into a resulting array. Otherwise, it builds the equation of the plane using the normal vector and one point. Then, using this equation, it calculates the

The conversion of from the triangular mesh into a rectangular one has one serious drawback. Since, triangular mesh may be highly non-uniform, it is possible to lost some points

Here, two points mark the node of the rectangular mesh and the closest point of the triangular mesh. The grayed triangles are the ones containing the closest point. Pay attention that the node of the rectangular mesh does not belong to any of the triangles. In this case, the program puts NaN into this node, thus creating the hole in the resulting field distribution (see

independently on the meshes resolution. Consider situation presented in figure 9(a)

Fig. 9. Field distribution of the leaky mode of the PhC fiber containing the defect of

triangular to rectangular mesh conversion

value of the data at the required point.

figure 9(b)).

#### **end**

*%otherwise, reading the coordinate points and %the field values within these points* x1=data(t(1, num\_t),1); y1=data(t(1, num\_t),2); z1=data(t(1, num\_t),3);

x2=data(t(2, num\_t),1); y2=data(t(2, num\_t),2); z2=data(t(2, num\_t),3);

x3=data(t(3, num\_t),1); y3=data(t(3, num\_t),2); z3=data(t(3, num\_t),3);

```
%Building the equation of the plane using vector
%cross product
n=cross([x1 x2 y1 y2 z1 z2], [x3 x2 y3 y2 z3 z2]);
A=n(1);
B=n(2);
C=n(3);
D= n(1)∗x1 n(2)∗y1 n(3)∗z1;
%and from the plane equation finding new value of the
%field intensity within the point of square mesh
datanew(yc,xc)=( D A∗xnew(xc) B∗ynew(yc))/C;
```
#### **end end**

```
%Refining the mesh using the interp2() function
[xrefined yrefined]=meshgrid(xmin:(xmax xmin)/ynpoints/10:xmax,...
                ymin:(ymax ymin)/ynpoints/10:ymax);
datanew_interp=interp2(xnew, ynew, datanew, xrefined, yrefined);
```
*%And plotting the resulting field distribution* surfc(datanew\_interp,'LineStyle', 'none');figure(gcf)

**function** tri\_area=tri\_area(p1,p2,p3) tri\_area=0.5∗(p1(1)∗p2(2) + p1(2)∗p3(1) + p3(2)∗p2(1) p3(1)∗p2(2) ... p3(2)∗p1(1) p2(1)∗p1(2));

#### **end**

Here in this program we are using the PDE Toolbox data structures to demonstrate the mesh conversion technique. In case of using other computation software, the results output may differ from the one shown here.

In this case, the program generates the regular square mesh. Then it scans over all the points in the mesh and for each point searches the closest one of the initial triangular mesh. After 28 MATLAB/Book4

**end**

*%otherwise, reading the coordinate points and*

*%Building the equation of the plane using vector*

n=cross([x1 x2 y1 y2 z1 z2], [x3 x2 y3 y2 z3 z2]);

*%and from the plane equation finding new value of the %field intensity within the point of square mesh*

datanew(yc,xc)=( D A∗xnew(xc) B∗ynew(yc))/C;

[xrefined yrefined]=meshgrid(xmin:(xmax xmin)/ynpoints/10:xmax,... ymin:(ymax ymin)/ynpoints/10:ymax);

tri\_area=0.5∗(p1(1)∗p2(2) + p1(2)∗p3(1) + p3(2)∗p2(1) p3(1)∗p2(2) ... p3(2)∗p1(1) p2(1)∗p1(2));

Here in this program we are using the PDE Toolbox data structures to demonstrate the mesh conversion technique. In case of using other computation software, the results output may

In this case, the program generates the regular square mesh. Then it scans over all the points in the mesh and for each point searches the closest one of the initial triangular mesh. After

datanew\_interp=interp2(xnew, ynew, datanew, xrefined, yrefined);

surfc(datanew\_interp,'LineStyle', 'none');figure(gcf)

*%the field values within these points*

x1=data(t(1, num\_t),1); y1=data(t(1, num\_t),2); z1=data(t(1, num\_t),3);

x2=data(t(2, num\_t),1); y2=data(t(2, num\_t),2); z2=data(t(2, num\_t),3);

x3=data(t(3, num\_t),1); y3=data(t(3, num\_t),2); z3=data(t(3, num\_t),3);

D= n(1)∗x1 n(2)∗y1 n(3)∗z1;

*%Refining the mesh using the interp2() function*

*%And plotting the resulting field distribution*

**function** tri\_area=tri\_area(p1,p2,p3)

differ from the one shown here.

*%cross product*

A=n(1); B=n(2); C=n(3);

**end end**

**end**

(b)

Fig. 9. Field distribution of the leaky mode of the PhC fiber containing the defect of triangular to rectangular mesh conversion

it finds the closest point, it uses the array "t" to find all the triangles containing this point. Then it checks which triangle contains current point of the rectangular mesh. For this reason it uses the method of the areas summation. If the one is not found, the point is outside the computation area and NaN is written into a resulting array. Otherwise, it builds the equation of the plane using the normal vector and one point. Then, using this equation, it calculates the value of the data at the required point.

The conversion of from the triangular mesh into a rectangular one has one serious drawback. Since, triangular mesh may be highly non-uniform, it is possible to lost some points independently on the meshes resolution. Consider situation presented in figure 9(a)

Here, two points mark the node of the rectangular mesh and the closest point of the triangular mesh. The grayed triangles are the ones containing the closest point. Pay attention that the node of the rectangular mesh does not belong to any of the triangles. In this case, the program puts NaN into this node, thus creating the hole in the resulting field distribution (see figure 9(b)).

for Photonics Applications 31

Results Processing in MATLAB for Photonics Applications 149

*%The first raw of the data corresponds to the X values and is copied*

*%the array data\_y is formed automatically from number of file*

*%is being copied into corresponding 3D array*

*%Making MATLAB to add new plot to the existing figure*

*%The actual 3D data is contained within the raws from 2 to 9 and*

*%The cycle for the number of levels now nl scans over another dimension*

*%Copying the data for a single band into a temporal 2D array*

surf(data\_x, data\_y, draw\_data,'MeshStyle','row');

set(gca,'XGrid', 'on','YGrid', 'on','ZGrid', 'on', ...

separately each surface corresponding to each PhC band.

xlabel('Wave vector', 'FontSize', 14, 'FontWeight', 'bold'); ylabel('Radius r/a', 'FontSize', 14, 'FontWeight', 'bold');

zlabel('Frequency, \omegaa/2\pic', 'FontSize', 14, 'FontWeight', 'bold');

The function simply merges the multiple files into a single data structure and then plots

Animations may become important when you are creating the presentation. In this case, dynamic representation of the results gives listener deep understanding of the problem.

'CameraPosition', [7.9 2.89 4.6]);

*%Closing the file* fclose(f);

*%Into array data\_x* data\_x(1,:)=data(1,:);

data\_y(nc+1)=nc∗0.01;

*%Creating the figure*

**end**

figure;

hold on;

**for** nl=1:8

**end**

ylim([0 0.5]); colormap('gray');

data\_z(nc+1,:,:)=data(2:9,:);

*%of the array (namely, to a band number).*

*%Creating the labels and decorating the plot*

draw\_data(:,:)=data\_z(:,nl,:);

*%Plotting the surface*

**8. Building animations**

This drawback, while being almost unessential for analysis, makes the method useless for data representation. In this case, the bes solution is to use something like neural networks interpolation which are presented in the following section.

#### **7. Multiple files data processing**

In most cases, you have the results in a form of a multiple files. For instance, when you computing the characteristics for the structure which parameters are varied during the computation process. In this case, each file corresponds to the specific value of the parameter.

#### **7.1 Gathering the multiple files with MATLAB**

The main problem in this case is to gather the files in a single data structure. This particular task depends on the kind of data. For example, in case of plane data you may want to plot the surface. However, when each file contains 3D data, the only possibility to visualize it is to make slices or animate the data. Particularly, when time-dependent field distribution is computed, it is wise to animate the results for better presentation.

In the first section we have considered how to open several files in a row. Here is the MATLAB program which gathers the data of the PhC band structures computed by RSoft into a single 3D array and plots the surfaces each one of which corresponds to a single PhC band. At this, the band structures of the PhC having different elements radius are suppose to be saved is a separate files with names like "bs\_separate00\_eigs\_TE.dat".

*%The program opens the files with the results of the band structure %computation as a function of the refractive index %obtained by RSoft and merges them into a single data structure. %Then each band is displayed as a single surface*

clear all;

*%Parts of the filename* name='bs\_single'; name\_fin='\_eigs\_TE.dat';

*%Cycle for the numbers of files* **for** nc=0:50

*%Forming the name* full\_name=strcat(name, num2str(nc,'%02i'), name\_fin);

*%Opening the file* f=fopen(full\_name,'r');

*%Skipping the header* skip=fscanf(f,'%s\n',5);

*%Reading the data from file into a 2D array* data=fscanf(f,'%f ',[9 16]);

30 MATLAB/Book4

This drawback, while being almost unessential for analysis, makes the method useless for data representation. In this case, the bes solution is to use something like neural networks

In most cases, you have the results in a form of a multiple files. For instance, when you computing the characteristics for the structure which parameters are varied during the computation process. In this case, each file corresponds to the specific value of the parameter.

The main problem in this case is to gather the files in a single data structure. This particular task depends on the kind of data. For example, in case of plane data you may want to plot the surface. However, when each file contains 3D data, the only possibility to visualize it is to make slices or animate the data. Particularly, when time-dependent field distribution is

In the first section we have considered how to open several files in a row. Here is the MATLAB program which gathers the data of the PhC band structures computed by RSoft into a single 3D array and plots the surfaces each one of which corresponds to a single PhC band. At this, the band structures of the PhC having different elements radius are suppose to be saved is a

interpolation which are presented in the following section.

computed, it is wise to animate the results for better presentation.

separate files with names like "bs\_separate00\_eigs\_TE.dat".

*%The program opens the files with the results of the band structure*

*%obtained by RSoft and merges them into a single data structure.*

full\_name=strcat(name, num2str(nc,'%02i'), name\_fin);

**7. Multiple files data processing**

**7.1 Gathering the multiple files with MATLAB**

*%computation as a function of the refractive index*

*%Then each band is displayed as a single surface*

clear all;

**for** nc=0:50

*%Parts of the filename* name='bs\_single';

name\_fin='\_eigs\_TE.dat';

*%Cycle for the numbers of files*

f=fopen(full\_name,'r');

data=fscanf(f,'%f ',[9 16]);

*%Reading the data from file into a 2D array*

*%Forming the name*

*%Opening the file*

*%Skipping the header* skip=fscanf(f,'%s\n',5);

```
%Closing the file
  fclose(f);
  %The first raw of the data corresponds to the X values and is copied
  %Into array data_x
  data_x(1,:)=data(1,:);
  %the array data_y is formed automatically from number of file
  data_y(nc+1)=nc∗0.01;
  %The actual 3D data is contained within the raws from 2 to 9 and
  %is being copied into corresponding 3D array
  data_z(nc+1,:,:)=data(2:9,:);
end
%Creating the figure
figure;
%Making MATLAB to add new plot to the existing figure
hold on;
%The cycle for the number of levels now nl scans over another dimension
%of the array (namely, to a band number).
for nl=1:8
  %Copying the data for a single band into a temporal 2D array
  draw_data(:,:)=data_z(:,nl,:);
  %Plotting the surface
  surf(data_x, data_y, draw_data,'MeshStyle','row');
end
%Creating the labels and decorating the plot
xlabel('Wave vector', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Radius r/a', 'FontSize', 14, 'FontWeight', 'bold');
zlabel('Frequency, \omegaa/2\pic', 'FontSize', 14, 'FontWeight', 'bold');
set(gca,'XGrid', 'on','YGrid', 'on','ZGrid', 'on', ...
                     'CameraPosition', [7.9 2.89 4.6]);
ylim([0 0.5]);
colormap('gray');
```
The function simply merges the multiple files into a single data structure and then plots separately each surface corresponding to each PhC band.

#### **8. Building animations**

Animations may become important when you are creating the presentation. In this case, dynamic representation of the results gives listener deep understanding of the problem.

for Photonics Applications 33

Results Processing in MATLAB for Photonics Applications 151

ylabel('Frequency, \omegaa/2\pic', 'FontSize', 14, 'FontWeight', 'bold');

Here in this program, three different kinds of animations in MATLAB are represented. Namely, using the command "drawnow", the animation is performed in the first time. Then, after it have been stored into an array "ani", the animation is played by the command "movie" which makes nothing more but displaying frame by frame the array. After this, the animation is stored into a ".avi" file which can be played separately in a video player or can be inserted

In this chapter, we have shared the experience of using the MATLAB for the photonic crystals modeling. Computation techniques considered in this chapter represent the examples aplicable to different fields of physics and electronics. The codes demonstrating the PhC computations, have been selected by authors following their experience and represent

Particularly, there have been given the brief review of several well-knows software projects for PhC devices investigation, presented the techniques for working with different useful data formats for different kinds of the data representation; presented several working codes for reading and analyzing the data from RSoft, COMSOL multiphysics, MPB and MEEP;

particular problems being solved during their investigation in the PhC field.

data\_y=data(2:9,:);

*%Drawing the bands*

plot(data\_x, data\_y(nl,:),'k','LineWidth',2);

xlabel('Wave vector', 'FontSize', 14, 'FontWeight', 'bold');

*%Setting up the parameters of the plot*

*%Saving current figure as a frame* ani(nc+1)=getframe(gcf);

*%Playing back all the array of frames*

*%Saving the animation into a video file* movie2avi(ani,'bs\_on\_radius');

set(gca, 'XGrid', 'on','YGrid', 'on');

*%Making MATLAB to display the plot immidiately*

**for** nl=1:8

**end**

hold on;

ylim([0 1.5]);

drawnow;

hold off;

movie(ani);

into a presentation.

**9. Conclusion**

**end**

Fig. 10. The result of merging the data from multiple files

MATLAB provides users with methods for creating basic animations capturing frame by frame from a graph area. To record an animation, you first have to perform the graphic output of the results. To output the data in the real time, the function "drawnow" is required. Then, using the function "getframe" each frame is saved into an array which then can be animated in MATLAB or stored into a file to insert into the presentation.

As an example, we will consider previous case with the band structure computation and will animate the band structure in time instead of representing it in a form of surfaces. Consider the following program:

*%The program reads the series of files from the specific location %and builds the animation which is then written into a file*

```
%Clearing the workspace
clear all;
%The first and the last constant parts of the name
name='bs_single';
name_fin='_eigs_TE.dat';
figure;
%The cycle through all the file numbers
for nc=0:50
  %Forming the filename
  full_name=strcat(name, num2str(nc,'%02i'), name_fin);
  %Reading the current file data
  f=fopen(full_name,'r');
  skip=fscanf(f,'%s\n',5);
  data=fscanf(f,'%f ',[9 16]);
  fclose(f);
  %Extracting the X data
```

```
data_x(1,:)=data(1,:);
```
*%Extracting the Y datafor 8 bands*

```
data_y=data(2:9,:);
```
32 MATLAB/Book4

**Radius r/a Wave vector**

MATLAB provides users with methods for creating basic animations capturing frame by frame from a graph area. To record an animation, you first have to perform the graphic output of the results. To output the data in the real time, the function "drawnow" is required. Then, using the function "getframe" each frame is saved into an array which then can be animated

As an example, we will consider previous case with the band structure computation and will animate the band structure in time instead of representing it in a form of surfaces. Consider

<sup>1</sup> <sup>0</sup> 0.1 0.2 0.3 0.4 0.5

0

Fig. 10. The result of merging the data from multiple files

in MATLAB or stored into a file to insert into the presentation.

*%The program reads the series of files from the specific location %and builds the animation which is then written into a file*

full\_name=strcat(name, num2str(nc,'%02i'), name\_fin);

*%The first and the last constant parts of the name*

the following program:

*%Clearing the workspace*

name='bs\_single';

name\_fin='\_eigs\_TE.dat';

*%Forming the filename*

*%Extracting the X data* data\_x(1,:)=data(1,:);

*%Extracting the Y datafor 8 bands*

*%The cycle through all the file numbers*

*%Reading the current file data* f=fopen(full\_name,'r'); skip=fscanf(f,'%s\n',5); data=fscanf(f,'%f ',[9 16]);

clear all;

figure;

**for** nc=0:50

fclose(f);

0

0.5

1

**Frequency,** ω**a/2**π**c**

1.5

2

0.5

```
%Drawing the bands
for nl=1:8
  plot(data_x, data_y(nl,:),'k','LineWidth',2);
  hold on;
end
```

```
%Setting up the parameters of the plot
ylim([0 1.5]);
xlabel('Wave vector', 'FontSize', 14, 'FontWeight', 'bold');
ylabel('Frequency, \omegaa/2\pic', 'FontSize', 14, 'FontWeight', 'bold');
set(gca, 'XGrid', 'on','YGrid', 'on');
```

```
%Making MATLAB to display the plot immidiately
drawnow;
```

```
%Saving current figure as a frame
ani(nc+1)=getframe(gcf);
hold off;
```
#### **end**

```
%Playing back all the array of frames
movie(ani);
```

```
%Saving the animation into a video file
movie2avi(ani,'bs_on_radius');
```
Here in this program, three different kinds of animations in MATLAB are represented. Namely, using the command "drawnow", the animation is performed in the first time. Then, after it have been stored into an array "ani", the animation is played by the command "movie" which makes nothing more but displaying frame by frame the array. After this, the animation is stored into a ".avi" file which can be played separately in a video player or can be inserted into a presentation.

#### **9. Conclusion**

In this chapter, we have shared the experience of using the MATLAB for the photonic crystals modeling. Computation techniques considered in this chapter represent the examples aplicable to different fields of physics and electronics. The codes demonstrating the PhC computations, have been selected by authors following their experience and represent particular problems being solved during their investigation in the PhC field.

Particularly, there have been given the brief review of several well-knows software projects for PhC devices investigation, presented the techniques for working with different useful data formats for different kinds of the data representation; presented several working codes for reading and analyzing the data from RSoft, COMSOL multiphysics, MPB and MEEP;

**Part 3** 

**Power Systems Applications** 

processing data represented in different forms such as plane data, 3D (rectangular and triangular mesh) data, the data represented by several files.

The examples of working codes given in the chapter allow wide range of investigators and, in the first place, PhD and master students, to start quickly the solution of their own original problems.

#### **10. References**

*COMSOL* (2011).

URL: *http://www.comsol.com*

Joannopoulos, J., Meade, R. & Winn, J. (1995). *Photonic Crystals: Molding the Flow of Light*, Princeton University Press, Princeton.

Lourtioz, J.-M., Benisty, H. & Berger, V. (2005). *Photonic Crystals, Towards Nanoscale Photonic Devices*, Springer-Verlag.

*MATLAB manual* (2011).

URL: *http://www.mathworks.com/help/techdoc/*

*MEEP* (2011).

URL: *http://ab-initio.mit.edu/wiki/index.php/Meep\_manual*

*MPB* (2011).

URL: *http://ab-initio.mit.edu/wiki/index.php/MPB\_manual*

*RSoft* (2011).

URL: *http://www.rsoftdesign.com/*

Sakoda, K. (2001). *Optical Properties of Photonic Crystals*, Springer-Verlag.

Sukhoivanov, I. & Guryev, I. (2009). *Photonic crystals – physics and practical modeling*, Optical sciences, 1 edn, Springer, Berlin.

Taflove, A. & Hagness, S. (2005). *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3 edn, Artech House, Norwood, MA.
