where *εi*, *μi* and *di* are the permittivity, permeability and thickness of the *i th* layer. The total wavevector in the *i*th layer is given by *ki* = √*εiμiω*/*c* where *ω* is the frequency of the electromagnetic field and *c* is the speed of light in vacuum. Also a normalized in plane wavevector is defined as *u* = *kx*/*k*1, where *kx* is the in plane wavevector component which is independent of the layer. The value of *ε* for the metal will in general be complex and frequency dependent, and is most easily obtained by fitting a high order polynomial to experimentally measured data (e.g. Drachev (2008) for silver). Alternatively for not to thin metal layers one may to a good approximation use a Drude-Sommerfeld model (e.g. Kittel (1995)) where the parameters can later be systematically varied to account for e.g. changes in temperature Elsayad and Heinze (2010b).

To determine the dispersion relation and the cut-off energy of a given structure one may also make use of the natural matrix calculation approaches of MATLAB. To do so one firstly writes a general expression for the parallel electric fields in each layer. For a layer-*i* extending from *z* = *zmin* to *z* = *zmax* this would be *A*(1) *i ekzi*(*z*+*zmin*) + *A*(2) *i e*−*kzi*(*z*+*zmax*), where *kzi* is the out of plane wavevector in the layer and *A*(1) *i* & *A*(2) *i* are at this point arbitrary complex coefficients. One may then write the corresponding magnetic induction fields for each layer using Maxwell's equations (i.e. by taking the normal derivative), so that for each layer one has two equations with two unknowns. This is done for each layer *i* = 1, 2...*N*. Finally one constructs a 2*N* × 2*N* with the coefficients for each layer constituting the entries in each row. The determinant of this matrix will then give the dispersion relation *ω*(*kx*) - see e.g. Dionne et al. (2008). Since *kx* is complex one can not simply assign a solver such as fzero to the task. It is it turns out often necessary to evaluate this by brute force - i.e. letting MATLAB evaluate the determinant for a given complex transverse wavevector and checking how close it is to zero. In general it is important to provide a good starting guess of Re(*kx*) and Im(*kx*) for each frequency. For the cases of interest, where the frequencies of interest are close to the cut-off frequency, a good starting guess for the former is the light-line Re(*kx*) = *nω*/*c* of the medium with refractive index *n* in which the mode becomes unbound, and then to progressively search further into the evanescent region, i.e. Re(*kx*) = *nω*/*c* + *δ* with *δ* ∼ 0.01*nω*/*c*. For the imaginary part Im(*kx*) - which diverges at the cut-off - it is most efficient to employ an algorithm which calculates the gradient of the determinant between two nearby values of Im(*kx*) [i.e. taking the difference in the determinant value and dividing by the difference in Im(*kx*)], and subsequently search for larger/smaller Im(*kx*) values if the determinant is smaller/larger than zero. The stability of the solution should be checked to assure Lim*δ*→0{*ω*(*kx* + *δ*) − *ω*(*kx* − *δ*) → 0}. The above mentioned approach typically requires that Im(*kx*) is well behaved and monotonous in the vicinity of the cut-off, which is fortunately often the case. Since the equation has to be satisfied for both Re(*kx*) and Im(*kx*) optimizations of both of these have to be performed simultaneously. For completeness an analysis can also be performed for transverse electric (TE) as well as transverse magnetic (TM) polarizations, although for the cases of the very thin non-magnetic plasmonic structures only the latter are usually relevant. For the case of highly asymmetric structures of few layers, coupling between the interfaces may be negligible and a single plasmonic mode at one interface can to a good approximation be assumed solely responsible for the dispersion of the cut-off mode Burke and Stegeman (1986). In this case the cut-off energy can for a non-magnetic structure be estimated by *Ec* = *hc*¯ *εjεk*[*εi*(*εj* + *εk*)]−1, where *εi* is the permittivity of the medium where the cut-off occurs, and *εj* and *εk* are the permittivities on either side of the interface at which the mode is originally localized.

The principle of MeSuMo requires that the cut-off of the structure falls within the emission spectrum of the emitter. Since many common fluorophores and dyes have a fairly broad emission spectrum this criteria can usually be met quite easily and is generally quite robust. To calculate the change in intensity at a given frequency one firstly defines the change in the decay rate of a fluorophore as a function of distance and frequency. This may (at not too small distances) be modeled quite well in the dipole approximation using the classical Chance-Prock-Silbey (CPS) model Chance (Prock and Silbey). Under these conditions the decay rate in the +*z*ˆ and −*z*ˆ direction (away from and towards the substrate respectively) for perpendicular (⊥) and parallel (||) electric dipole like emitters can be obtained by:

$$\begin{aligned} \hat{\Gamma}\_{+}^{\perp} &= q - \frac{3q}{4} \text{Im} \int\_{0}^{1} (\mathcal{Z}\_{1} + \mathcal{Z}\_{2}) \, du \\ \hat{\Gamma}\_{-}^{\perp} &= \frac{3q}{4} \text{Im} \left( \int\_{0}^{1} \mathcal{Z}\_{1} du - \int\_{1}^{\infty} \mathcal{Z}\_{2} du \right) \\ \hat{\Gamma}\_{+}^{\parallel} &= q - \frac{3q}{8} \text{Im} \int\_{0}^{1} (\mathcal{Z}\_{3} - \mathcal{Z}\_{4}) du \end{aligned} \tag{6}$$

$$\hat{\mathbf{f}}\_{-}^{||} = \frac{3q}{8} \text{Im} \left( \int\_{0}^{1} \mathcal{Z}\_{3} du + \int\_{1}^{\infty} \mathcal{Z}\_{4} du \right) \tag{7}$$

in which the integrands are

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

complete equations (e.g. as a function) into MATLAB. For example, for the case of a four

*jkl* exp[2*i*(*εjμj*(*εiμi*)−1 − *u*2)

*jkl* exp[2*i*(*εjμj*(*εiμi*)−1 − *u*2)

*kl* exp[2*i*(*εkμk*(*εiμi*)−1 − *u*2)

*kl* exp[2*i*(*εkμk*(*εiμi*)−1 − *u*2)

*ij* = *μj*(1 − *u*2)1/2 − *μi*[*εjμj*(*εiμi*)−1 − *u*2]

*ij* = *εj*(1 − *u*2)1/2 − *εi*[*εjμj*(*εiμi*)−1 − *u*2]

total wavevector in the *i*th layer is given by *ki* = √*εiμiω*/*c* where *ω* is the frequency of the electromagnetic field and *c* is the speed of light in vacuum. Also a normalized in plane wavevector is defined as *u* = *kx*/*k*1, where *kx* is the in plane wavevector component which is independent of the layer. The value of *ε* for the metal will in general be complex and frequency dependent, and is most easily obtained by fitting a high order polynomial to experimentally measured data (e.g. Drachev (2008) for silver). Alternatively for not to thin metal layers one may to a good approximation use a Drude-Sommerfeld model (e.g. Kittel (1995)) where the parameters can later be systematically varied to account for e.g. changes in temperature

To determine the dispersion relation and the cut-off energy of a given structure one may also make use of the natural matrix calculation approaches of MATLAB. To do so one firstly writes a general expression for the parallel electric fields in each layer. For a layer-*i* extending

coefficients. One may then write the corresponding magnetic induction fields for each layer using Maxwell's equations (i.e. by taking the normal derivative), so that for each layer one has two equations with two unknowns. This is done for each layer *i* = 1, 2...*N*. Finally one constructs a 2*N* × 2*N* with the coefficients for each layer constituting the entries in each row. The determinant of this matrix will then give the dispersion relation *ω*(*kx*) - see e.g. Dionne et al. (2008). Since *kx* is complex one can not simply assign a solver such as fzero to the task. It is it turns out often necessary to evaluate this by brute force - i.e. letting MATLAB evaluate the determinant for a given complex transverse wavevector and checking how close it is to zero. In general it is important to provide a good starting guess of Re(*kx*) and Im(*kx*) for each frequency. For the cases of interest, where the frequencies of interest are close to the cut-off frequency, a good starting guess for the former is the light-line Re(*kx*) = *nω*/*c* of the medium with refractive index *n* in which the mode becomes unbound,

*i ekzi*(*z*+*zmin*) + *A*(2)

*i* & *A*(2)

1 2 *kidj*]

> 1 2 *kidj*]

1 2 *kidk*]

1 2 *kidk*]

1/2

1/2

*μj*(1 − *u*2)1/2 + *μi*[*εjμj*(*εiμi*)−1 − *u*2]1/2 (4)

*εj*(1 − *u*2)1/2 + *εi*[*εjμj*(*εiμi*)−1 − *u*2]1/2 (5)

(2)

(3)

*th* layer. The

*i e*−*kzi*(*z*+*zmax*), where *kzi* is the

*i* are at this point arbitrary complex

layered structure (layers indexed by *i*, *j*, *k*, *l*) this would be:

*P*(*S*) *ij* + *r*

1 + *r P*(*S*) *ij r P*(*S*)

*P*(*S*) *jk* + *r*

1 + *r P*(*S*) *jk r P*(*S*)

*P*(*S*)

*P*(*S*)

where *εi*, *μi* and *di* are the permittivity, permeability and thickness of the *i*

*r P*(*S*) *ijkl* = *r*

*r P*(*S*) *jkl* = *r*

Elsayad and Heinze (2010b).

*rS*

*rP*

from *z* = *zmin* to *z* = *zmax* this would be *A*(1)

out of plane wavevector in the layer and *A*(1)

with

and

$$\mathcal{Z}\_1 = 2u^3 \frac{1 - |r^P|^2}{(u^2 - 1)^{\frac{1}{2}}} \qquad \mathcal{Z}\_2 = u^3 \frac{r^P}{(u^2 - 1)^{\frac{1}{2}}} \exp[2(u^2 - 1)^{\frac{1}{2}}k\_1 z]$$

$$\mathcal{Z}\_3 = u \frac{(1 - |r^S|^2) + (1 - u^2)(1 - |r^P|^2)}{(u^2 - 1)^{\frac{1}{2}}} \qquad \mathcal{Z}\_4 = 2u \frac{r^S + (1 - u^2)r^P}{(u^2 - 1)^{\frac{1}{2}}} \exp[2(u^2 - 1)^{\frac{1}{2}}k\_1 z] \tag{8}$$

proceed in several ways. It is of no practical relevance in our applications if one requires information on the position and orientation of each and every contributing fluorophore. We thus can proceed only via an effective medium approximation. We will assume a density *ρ* of fluorophores which may only vary in the axial (*z*) direction. In the case for *ρ*(0 < *z* < *zmax*) =

369 Using MATLAB to Achieve Nanoscale Optical Sectioning

Alternatively one may represent the contribution from a collection of fluorophores by assigning a negative imaginary response function to the top layer (layer-1) in which they are immersed - i.e. set Im(*ε*1)< 0 for non-magnetic structures. Whilst this may be computationally simpler for constant densities, for the case of a non-trivial *z* dependence the calculations

For imaginary (1 − *u*2)1/2*k*1 the coupling fields are propagating and equation (12) can be reduced to the result for e.g. the field in a large semi-transparent 1D box filled with emitters. For real (1 − *u*2)1/2*k*1 the field will on the other hand decay rapidly with increasing *z* as one would expect. The contribution to *I*(*ωab*) from equation (12) will in itself however only be

There will be an additional second order contribution to the excitation field and hence the emission intensity at *ωab* from excited surface modes of energy *ωax* where *x* �= *b*. This contribution is only significant because the SPP modes have very short decay times - typically Γ ∼ O(10-100fs), such that their uncertainty broadening is relatively large - O(1 − 10eV). For *ωax* ∼ *ωab*, or more specifically |*ωax* − *ωab*|∼O(meV), they will contribute significantly to **E***R*(*ωab*, *z*). It follows that when the energy *ωab* corresponds to mainly real values of (1 − *u*2)1/2*k*1 (evanescent fields) whereas *ωax* to mainly imaginary values of (1 − *u*2)1/2*k*1 (propagating fields) this contribution becomes large near the surface as fluorophores far away from the substrate can effectively pump those nearby. To the lowest order the contribution

> Γ*aβ*(*u*, *z*) Γ*aβ*

and Γ*ab*(*u*, *z*) and **E***r*(*ω*, *u*, *z*) are the corresponding values of the total decay rate and reflected field as a function of transverse wavevector [i.e. the integrands of equations (2) and (2) &/or (3)]. We note that equation (13) accounts for the energy overlap between a state of energy *ωb* and a continuum of states whose energy is given by the dummy variable *ωβ*. The *ωα* integration is performed over a range of frequencies in the vicinity of the cut-off (typically corresponding to a wavelength range of 30-80nm), over which the energy range coincides with the emission spectrum of the fluorophore (i.e. are allowed transitions for the particular fluorophore). The weighting (Franck Condon) factors *fαβ* are empirically determined from

Γ*aβ*(*u*, *z*)

(1 − *u*2)1/2*k*1

Γ*ab*(*u*, *z*) Γ*a*(*z*) **E**0

*Laβ***E***r*(*ωaβ*, *u*, *z*)*Lβαδ*(*ωα* − *ωax*) (13)

(*ωa* − *ωβ*)2 + [Γ*aβ*(*u*, *z*)/2]2 (14)

*r* (*ωab*, *u*, *z*/2). (12)

−(1−*u*2)1/2*k*1*zmax* ) *ρ*

in the Vicinity of Metamaterial Substrates by Simulating Emitter-Substrate Interactions

become significantly more elaborate. We thus proceed with the former method.

*ρ* and *ρ*(*z* > *zmax*) = 0 one obtains by integrating over *z* ={0, ∞}:

**E***r*(*ωab*, *u*, *z*)=(1 − *e*

small even for large *ρ* (> 0.01nm−1).

can be written as:

where

**E**0

*r* (*ωax*, *u*, *z*) =

*dωα*

*Laβ* = 1 2*π*

*dωβ faβ*

where *q* is the quantum yield of the emitter, the superscripts *P*(*S*) denote the transverse electric(magnetic) cases, and the *z* and *ω* dependence of Γˆ ⊥(||) +(−) is implied. The integrals can in most cases quite easily be evaluated numerically (e.g. using quad). It may be the case that the structures have resonances (as can be determined by plotting the integrands as a function of *u*), in which case one may want to re-write the integrals to assure one integrates in a complex plane around them and account for the pole(s) in the conventional fashion. Integration can in almost all cases safely be cut-off at *u* ∼ 1000, with the contribution beyond *u* ≈ 100 becoming negligibly small for *z* > 10nm. In what follows we will sometimes make use of the notation Γ*ij* = 1 3 [(Γ⊥ + + Γ⊥ −) + 2(Γ|| + + Γ|| −)]*ωij* to signify the isotropic average total decay rate from state *i* to *j* separated in energy by *ωij* = *ωi* − *ωj*. The energy dependence of Γ*ij* is understood to be contained in *rS*,*P* = *rS*,*P*(*ωij*), *u* = *u*(*ωij*) and *k*1 = *k*1(*ωij*).

Besides the modification in the decay rate which has been described above, the measured emission intensity of an emitter will also depend on it's excitation rate. In particular the measurable emission intensity for the transition from the excited state *a* to a lower state *b*, which corresponds to a wavelength *λab* = 2*πc*/*ωab* will be given by:

$$I(\omega\_{ab}) \propto |\hat{\mu}\_a \cdot \mathbf{E}\_{ex}(\omega\_{ab})|^2 f\_{ab} \Gamma^R\_{ab} \Gamma^{-1}\_a \tag{9}$$

where **E***ex*(*ωab*) is the excitation field at the frequency *ωab*, and *μ*ˆ *a* is the dipole moment of the excited state *a*. Γ*R ab* and Γ*a* are the radiative decay rate for the transition *a* → *b* and the total decay rate of the state *a*. *fab* is the Franck Condon coefficient between *a* and *b*, which is related to the rotational and vibrational configurations of the two energy states relative to each other, and can be thought of as the intrinsic favourability of the transition. For a realistic emitter there will in general be many other significant "down transitions" originating from the same excited state. For the decay of the state *a* to any other state *x*, one may thus write:

$$I(\omega\_{ax}) = \frac{f\_{ax}}{f\_{ab}} \frac{\Gamma\_{ax}^{R}}{\Gamma\_{ab}^{R}} I(\omega\_{ab}) \tag{10}$$

Subsequent decay of all states *x* (and state *b*) to a common ground state is assumed to proceed non-radiatively and at a rate several orders of magnitude larger than the radiative decay rate. For the case of spontaneous emission we consider it can for all intense and purposes be assumed to be instantaneous. If we neglect direct resonant energy transfer between emitters (i.e. assume sufficient dilution), then the total excitation field with a frequency *ωab* at the location of the fluorophore will be the sum of the external incident field **E**0(*ωab*), the reflection of this field from the structure *re*−2(1−*u*2)1/2*k*1*z***E**0(*ωab*) where *r* is the reflection coefficient, and the field generated from excitations in the material by all the fluorophores **E***r*(*ωab*, *z*), i.e.

$$\mathbf{E}\_{\rm ex}(\omega\_{ab\prime}\boldsymbol{\mu}, \boldsymbol{z}) = [1 + re^{-2(1-\boldsymbol{u}^{2})^{1/2}k\_{1}z}]\mathbf{E}^{0}(\omega\_{ab}) + \mathbf{E}\_{r}(\omega\_{ab\prime}\boldsymbol{\mu}, \boldsymbol{z})\tag{11}$$

The polarization for the reflected field can in most cases be taken as that of the incident excitation field. For the case of a single fluorophore **E***r*(*ωab*, *z*) is given by the integrands of equations (6) and/or (7), which are identical to the reflected field for a given transverse wavevector component *u*, and will henceforth be written as **E**0 *r* (*ωab*, *u*, *z*). To calculate the reflected field from all of the fluorophores at various distances from the substrate one can proceed in several ways. It is of no practical relevance in our applications if one requires information on the position and orientation of each and every contributing fluorophore. We thus can proceed only via an effective medium approximation. We will assume a density *ρ* of fluorophores which may only vary in the axial (*z*) direction. In the case for *ρ*(0 < *z* < *zmax*) = *ρ* and *ρ*(*z* > *zmax*) = 0 one obtains by integrating over *z* ={0, ∞}:

$$\mathbf{E}\_{r}(\omega\_{ab\nu}\boldsymbol{u},\boldsymbol{z}) = (1 - e^{-(1-\boldsymbol{u}^{2})^{1/2}k\_{1}z\_{\max}}) \frac{\rho}{(1-\boldsymbol{u}^{2})^{1/2}k\_{1}} \frac{\Gamma\_{ab}(\boldsymbol{u},\boldsymbol{z})}{\Gamma\_{a}(\boldsymbol{z})} \mathbf{E}\_{r}^{0}(\omega\_{ab\nu}\boldsymbol{u},\boldsymbol{z}/2). \tag{12}$$

Alternatively one may represent the contribution from a collection of fluorophores by assigning a negative imaginary response function to the top layer (layer-1) in which they are immersed - i.e. set Im(*ε*1)< 0 for non-magnetic structures. Whilst this may be computationally simpler for constant densities, for the case of a non-trivial *z* dependence the calculations become significantly more elaborate. We thus proceed with the former method.

For imaginary (1 − *u*2)1/2*k*1 the coupling fields are propagating and equation (12) can be reduced to the result for e.g. the field in a large semi-transparent 1D box filled with emitters. For real (1 − *u*2)1/2*k*1 the field will on the other hand decay rapidly with increasing *z* as one would expect. The contribution to *I*(*ωab*) from equation (12) will in itself however only be small even for large *ρ* (> 0.01nm−1).

There will be an additional second order contribution to the excitation field and hence the emission intensity at *ωab* from excited surface modes of energy *ωax* where *x* �= *b*. This contribution is only significant because the SPP modes have very short decay times - typically Γ ∼ O(10-100fs), such that their uncertainty broadening is relatively large - O(1 − 10eV). For *ωax* ∼ *ωab*, or more specifically |*ωax* − *ωab*|∼O(meV), they will contribute significantly to **E***R*(*ωab*, *z*). It follows that when the energy *ωab* corresponds to mainly real values of (1 − *u*2)1/2*k*1 (evanescent fields) whereas *ωax* to mainly imaginary values of (1 − *u*2)1/2*k*1 (propagating fields) this contribution becomes large near the surface as fluorophores far away from the substrate can effectively pump those nearby. To the lowest order the contribution can be written as:

$$\mathbf{E}\_r^0(\omega\_{a\mathbf{x}\prime}\boldsymbol{\mu},z) = \int d\omega\_a \int d\omega\_\beta f\_{a\beta} \frac{\Gamma\_{a\beta}(\boldsymbol{\mu},z)}{\Gamma\_{a\beta}} \mathcal{L}\_{a\beta} \mathbf{E}\_r(\omega\_{a\beta\prime}\boldsymbol{\mu},z) \mathcal{L}\_{\beta a} \delta(\omega\_a - \omega\_{a\mathbf{x}}) \tag{13}$$

where

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

where *q* is the quantum yield of the emitter, the superscripts *P*(*S*) denote the transverse

most cases quite easily be evaluated numerically (e.g. using quad). It may be the case that the structures have resonances (as can be determined by plotting the integrands as a function of *u*), in which case one may want to re-write the integrals to assure one integrates in a complex plane around them and account for the pole(s) in the conventional fashion. Integration can in almost all cases safely be cut-off at *u* ∼ 1000, with the contribution beyond *u* ≈ 100 becoming negligibly small for *z* > 10nm. In what follows we will sometimes make use of the notation

*i* to *j* separated in energy by *ωij* = *ωi* − *ωj*. The energy dependence of Γ*ij* is understood to be

Besides the modification in the decay rate which has been described above, the measured emission intensity of an emitter will also depend on it's excitation rate. In particular the measurable emission intensity for the transition from the excited state *a* to a lower state *b*,

where **E***ex*(*ωab*) is the excitation field at the frequency *ωab*, and *μ*ˆ *a* is the dipole moment of the

decay rate of the state *a*. *fab* is the Franck Condon coefficient between *a* and *b*, which is related to the rotational and vibrational configurations of the two energy states relative to each other, and can be thought of as the intrinsic favourability of the transition. For a realistic emitter there will in general be many other significant "down transitions" originating from the same

*fab*

Subsequent decay of all states *x* (and state *b*) to a common ground state is assumed to proceed non-radiatively and at a rate several orders of magnitude larger than the radiative decay rate. For the case of spontaneous emission we consider it can for all intense and purposes be assumed to be instantaneous. If we neglect direct resonant energy transfer between emitters (i.e. assume sufficient dilution), then the total excitation field with a frequency *ωab* at the location of the fluorophore will be the sum of the external incident field **E**0(*ωab*), the reflection of this field from the structure *re*−2(1−*u*2)1/2*k*1*z***E**0(*ωab*) where *r* is the reflection coefficient, and the field generated from excitations in the material by all the fluorophores **E***r*(*ωab*, *z*), i.e.

The polarization for the reflected field can in most cases be taken as that of the incident excitation field. For the case of a single fluorophore **E***r*(*ωab*, *z*) is given by the integrands of equations (6) and/or (7), which are identical to the reflected field for a given transverse

reflected field from all of the fluorophores at various distances from the substrate one can

wavevector component *u*, and will henceforth be written as **E**0

Γ*R ax* Γ*R ab* +(−) is implied. The integrals can in

*a* (9)

*I*(*ωab*) (10)

*r* (*ωab*, *u*, *z*). To calculate the

−)]*ωij* to signify the isotropic average total decay rate from state

2 *fab*Γ*R*

*ab* and Γ*a* are the radiative decay rate for the transition *a* → *b* and the total

**E***ex*(*ωab*, *u*, *z*)=[1 + *re*−2(1−*u*2)1/2*k*1*z*]**E**0(*ωab*) + **E***r*(*ωab*, *u*, *z*) (11)

*ab*Γ−1

electric(magnetic) cases, and the *z* and *ω* dependence of Γˆ ⊥(||)

Γ*ij* = 1

3 [(Γ⊥

excited state *a*. Γ*R*

+ + Γ⊥

−) + 2(Γ||

+ + Γ||

contained in *rS*,*P* = *rS*,*P*(*ωij*), *u* = *u*(*ωij*) and *k*1 = *k*1(*ωij*).

which corresponds to a wavelength *λab* = 2*πc*/*ωab* will be given by:

*I*(*ωab*) ∝ |*μ*ˆ *a* · **E***ex*(*ωab*)|

excited state. For the decay of the state *a* to any other state *x*, one may thus write:

*I*(*ωax*) = *fax*

$$L\_{a\S} = \frac{1}{2\pi} \frac{\Gamma\_{a\S}(\mu, z)}{(\omega\_a - \omega\_\beta)^2 + [\Gamma\_{a\S}(\mu, z)/2]^2} \tag{14}$$

and Γ*ab*(*u*, *z*) and **E***r*(*ω*, *u*, *z*) are the corresponding values of the total decay rate and reflected field as a function of transverse wavevector [i.e. the integrands of equations (2) and (2) &/or (3)]. We note that equation (13) accounts for the energy overlap between a state of energy *ωb* and a continuum of states whose energy is given by the dummy variable *ωβ*. The *ωα* integration is performed over a range of frequencies in the vicinity of the cut-off (typically corresponding to a wavelength range of 30-80nm), over which the energy range coincides with the emission spectrum of the fluorophore (i.e. are allowed transitions for the particular fluorophore). The weighting (Franck Condon) factors *fαβ* are empirically determined from

**Listing 1. Numerical integration**

2 % -- material permittivities 3 sio2 = 2.13; glass = 2.46; 4 water = 1.78; si3n4 = 4.15;

7 kx = 0 : 0.000001 : 0.02\*pi;

10 % -- thicknesses for each layer

16 lower = 0; % limits for integrals

20 function [I1] = integrand\_I1(u)

28 function r\_1234 = r\_1234\_PS(u,polariz) 29 global d; global e; global k; global mu;

39 num = r\_12\_P + r\_234\_P.\*exp\_term; 40 den = 1 + r\_12\_P .\* r\_234\_P.\*exp\_term;

43 num = r\_12\_S + r\_234\_S.\*exp\_term; 44 den = 1 + r\_12\_S .\* r\_234\_S.\*exp\_term;

49 function r\_234 = r\_234\_PS(u,polariz) 50 global d; global e; global k; global mu;

15 % -- integrate Integrand\_1

23 A=1- abs(r\_P).^2; 24 B = sqrt(u.^2 - 1); 25 I1 = 2.\*u.^3.\*A./B;

 r\_12\_P = r\_ij\_PS(u,1,2,"P"); r\_12\_S = r\_ij\_PS(u,1,2,"S"); r\_234\_P = r\_234\_PS(u,"P"); r\_234\_S = r\_234\_PS(u,"S");

6 w\_vec = 0 : 0.01 : 4;

11 mu = [1 1 1 1]; 12 d = [1e6 15 15 1e6];

17 upper = 1;

5 % -- ranges for frequency & wavevectors

9 % -- permittivities & permeabilities &

13 e = [ sio2 permitivity\_silver(w,d(2)) si3n4 glass];

18 b\_up\_a = quad(@(u) integrand\_I1(u),lower,upper) );

31 exp\_term = exp(2i.\*sqrt(e(2).\*mu(2).\*(e(1).\*mu(1)).^(-1)-u.^2).\*k(1).\*d(2));

371 Using MATLAB to Achieve Nanoscale Optical Sectioning

in the Vicinity of Metamaterial Substrates by Simulating Emitter-Substrate Interactions

21 [r\_S,r\_P] = r\_1234\_vec = r\_1234\_PS(u,'P');

1 function main()

8

14

19

22

30

36

37 switch(polariz) 38 case {"P"}

42 case {"S"}

46 end 47 end 48

51

41 r\_1234 = num ./ den;

45 r\_1234 = num ./ den;

26 end 27

the emission spectrum of an isolated fluorophore scaled so that the value at *ωab* it is unity. We note that it is usually practical or necessary to convert this integration into a discrete sum with an energy spacing of Δ*ω* ∼ 1meV. Information particular to the type of fluorophore (e.g. energy level positions/separations and corresponding Franck Condon coefficients) may be implemented when constructing the summation.

Since we will not focus on single molecule studies we assume an isotropic average of dipole orientations for all cases. This can most easily be done by setting **E***iso R* = (1/3)**E**⊥ *R* + (2/3)**E**|| *R* and/or Γ*iso* = (1/3)Γ⊥ + (2/3)Γ|| as appropriate. We also note to calculate the final measurable emission intensity the integral for Γ*R* which appears in equations 9 & 10 (and is given by equations 6 & 7) should only be performed up to √*εi*·NA, where NA is the numerical aperture of the measuring objective. It is possible to now predict the change in the measurable emission intensity for arbitrary *z* and fluorophores, near a metamaterial structure. In figure 1 we outline the basic backbone of a program that calculates the matrix enhan(j,i) which quantifies the measurable fluorescence enhancement at discrete distances (index j) and wavelengths (index i). The order of computations (top to bottom) is such that several independent calculations can depending on available resources be performed in parallel. The calculation is performed over all transverse wavevectors up to at least |*u*| ∼ 10*z*−1 and summed at the end (with the corresponding weighting factors to mimic integration).

Fig. 1. Structure of program for calculating a matrix that predicts fluorescence enhancement

We now show how to implement equations 6 &7 (with the integrands defined by equation 9) in MATLAB, which are required for solving equation 10. To focus on the most relevant part of the code, we expect the reader to have some familiarity with basic MATLAB programming syntax (i.e. definition of functions and vector-based MATLAB notation). An example code for calculating the reflection coefficients [2 - 5] and subsequently performing the integration may look like:

**Listing 1. Numerical integration**

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

the emission spectrum of an isolated fluorophore scaled so that the value at *ωab* it is unity. We note that it is usually practical or necessary to convert this integration into a discrete sum with an energy spacing of Δ*ω* ∼ 1meV. Information particular to the type of fluorophore (e.g. energy level positions/separations and corresponding Franck Condon coefficients) may be

Since we will not focus on single molecule studies we assume an isotropic average of dipole

and/or Γ*iso* = (1/3)Γ⊥ + (2/3)Γ|| as appropriate. We also note to calculate the final measurable emission intensity the integral for Γ*R* which appears in equations 9 & 10 (and is given by equations 6 & 7) should only be performed up to √*εi*·NA, where NA is the numerical aperture of the measuring objective. It is possible to now predict the change in the measurable emission intensity for arbitrary *z* and fluorophores, near a metamaterial structure. In figure 1 we outline the basic backbone of a program that calculates the matrix enhan(j,i) which quantifies the measurable fluorescence enhancement at discrete distances (index j) and wavelengths (index i). The order of computations (top to bottom) is such that several independent calculations can depending on available resources be performed in parallel. The calculation is performed over all transverse wavevectors up to at least |*u*| ∼ 10*z*−1 and

summed at the end (with the corresponding weighting factors to mimic integration).

calculate cut-off energy ω(c) define Ω(i), i=1,2...10 in vicinity of ω(ab) and ω(c)

dispersion

i+1

j+1

look like:

structure (d1,epsilon1, mu1, etc.)

calculate decay rates

Fig. 1. Structure of program for calculating a matrix that predicts fluorescence enhancement We now show how to implement equations 6 &7 (with the integrands defined by equation 9) in MATLAB, which are required for solving equation 10. To focus on the most relevant part of the code, we expect the reader to have some familiarity with basic MATLAB programming syntax (i.e. definition of functions and vector-based MATLAB notation). An example code for calculating the reflection coefficients [2 - 5] and subsequently performing the integration may

z-range z(j) j=1,2..100

define reflection coeff. & decay rates

calculate reflection coefficients

calculate total enhancement for Ω(i)

total enhancement matrix as function of z(j) and Ω(i) enhan(i,j)

calculate (12)-(14)

calculate (11)

fluorophore density rho

excitation field

*R* = (1/3)**E**⊥

*R* + (2/3)**E**||

*R*

implemented when constructing the summation.

normalize & fit determine ω(ab) and ω range

& absorption spectrum metal

fluorophore emission

orientations for all cases. This can most easily be done by setting **E***iso*

```
1 function main()
2 % -- material permittivities
3 sio2 = 2.13; glass = 2.46;
4 water = 1.78; si3n4 = 4.15;
5 % -- ranges for frequency & wavevectors
6 w_vec = 0 : 0.01 : 4;
7 kx = 0 : 0.000001 : 0.02*pi;
8
9 % -- permittivities & permeabilities &
10 % -- thicknesses for each layer
11 mu = [1 1 1 1];
12 d = [1e6 15 15 1e6];
13 e = [ sio2 permitivity_silver(w,d(2)) si3n4 glass];
14
15 % -- integrate Integrand_1
16 lower = 0; % limits for integrals
17 upper = 1;
18 b_up_a = quad(@(u) integrand_I1(u),lower,upper) );
19
20 function [I1] = integrand_I1(u)
21 [r_S,r_P] = r_1234_vec = r_1234_PS(u,'P');
22
23 A=1- abs(r_P).^2;
24 B = sqrt(u.^2 - 1);
25 I1 = 2.*u.^3.*A./B;
26 end
27
28 function r_1234 = r_1234_PS(u,polariz)
29 global d; global e; global k; global mu;
30
31 exp_term = exp(2i.*sqrt(e(2).*mu(2).*(e(1).*mu(1)).^(-1)-u.^2).*k(1).*d(2));
32 r_12_P = r_ij_PS(u,1,2,"P");
33 r_12_S = r_ij_PS(u,1,2,"S");
34 r_234_P = r_234_PS(u,"P");
35 r_234_S = r_234_PS(u,"S");
36
37 switch(polariz)
38 case {"P"}
39 num = r_12_P + r_234_P.*exp_term;
40 den = 1 + r_12_P .* r_234_P.*exp_term;
41 r_1234 = num ./ den;
42 case {"S"}
43 num = r_12_S + r_234_S.*exp_term;
44 den = 1 + r_12_S .* r_234_S.*exp_term;
45 r_1234 = num ./ den;
46 end
47 end
48
49 function r_234 = r_234_PS(u,polariz)
50 global d; global e; global k; global mu;
51
```
polynomial and evaluated at the discrete wavelengths indexed by i], and saved as a 1 × 10

373 Using MATLAB to Achieve Nanoscale Optical Sectioning

in the Vicinity of Metamaterial Substrates by Simulating Emitter-Substrate Interactions

Data analysis is performed on a pixel by pixel basis. For each pixel the spectrum area is normalized, fitted with a high order polynomial, evaluated at i discrete wavelengths, and written into a 1 × 10 matrix data(i). The ratio enhanobserved(i) = data(i)/Ilambda0(i) is then calculated. A 1 × *j* matrix zn(j) is defined which will represent the fraction of the emitters in the studied area at a distance corresponding to the index j from the interface. The task is now to solve enhan\*zn-enhanobserved=0, which is in principle straight forward using a function such as e.g. linsolve. The result is best presented in a plot of zn as a function of z for each or a collection of pixels, which will show the distribution of fluorophore distances from the substrate at this particular lateral position. Depending on the sample studied and the information of interest this may be fitted with a suitable distribution function to obtain an average and spread in distances of emitters from

In many cases the initial fluorophore-concentration, metamaterial or coupling parameters may not be the best or correct parameters for the fit. This is usually evident when the results show an unexpected distribution for zn. Distances less than 10nm should be ignored as the approximations of the models (dipole approximation and local dielectric functions) are no longer valid at such small distances. The most significant source of uncertainty is the predefined and assumed constant value of the fluorophore concentration rho. To account for this one may redefine rho as being proportional to the integrated area of zn over z, and repeat the calculation of zn. Also to account for an axial variation in rho - as would be relevant for samples where the fluorophore density varies strongly with distance z - one may define rho(j) = rho\*(a + zn(j)/norm(zn)), where a is normalization constant chosen so that on summation over the distances (*i*) one recovers the total intensity. In this case all subsequent calculations have to explicitly account for the additional distance dependence in *rho (j)*, which significantly increases the computational time. This may be efficiently achieved by defining rho(j) in a smaller parameter space indexed by n that corresponds to dividing rho into components where rho(j) is constant and non-zero only over a specified region. Once zn is calculated in this manner the optimization may be repeated for the new rho until the the calculation of the distribution of zn is stable and smooth or varies on the expected

Whilst elaborate custom substrates may be designed with effective parameters which yield the desired cut-off energy and dispersion relation, it is often possible to use combinations of no more than 3-4 layers with no lateral structuring to achieve satisfactory results. Here we give an example where we use a sample consisting of a thin silver film2 (15nm) deposited on a thick quartz substrate (*ε* = 2.13), which is subsequently coated with a high permittivity dielectric film (*ε* = 4.15). The sample is immersed in a medium with a refractive index larger than that of

2 Smoothness of the silver is guaranteed by using a Germanium wetting layer and is less than 0.4nm

matrix Ilambda0(i).

the substrate within the analyzed region.

RMS as measured by Atomic Force Microscopy.

**4. Optimizing the model and code**

scale.

**5. Example**

```
52 exp_term = exp(2i.*k(1).*d(3).*sqrt(e(3).*mu(3).*(e(1).*mu(1)).^-1-u.^2));
53
54 switch(polariz)
55 case {"P"}
56 r_23_P = r_ij_PS(u,2,3,"P");
57 r_34_P = r_ij_PS(u,3,4,"P");
58 num = r_23_P + r_34_P.*exp_term;
59 den = 1 + r_23_P .* r_34_P.*exp_term;
60 r_234 = num ./ den;
61 case {"S"}
62 r_23_S = r_ij_PS(u,2,3,"S");
63 r_34_S = r_ij_PS(u,3,4,"S");
64 num = r_23_S + r_34_S.*exp_term;
65 den = 1 + r_23_S .* r_34_S.*exp_term;
66 r_234 = num ./ den;
67 end
68 end
69
70 function r_ij = r_ij_PS(u,i,j,polariz)
71 global e; global mu;
72
73 A = (1 - u.^2).^0.5;
74 B = ( e(j).*mu(j).* (e(i).*mu(i)).^-1 - u.^2).^0.5;
75
76 switch(polariz)
77 case {"S"}
78 num = mu(j).*A - mu(i).* B;
79 den = mu(j).*A + mu(i).* B;
80 r_ij = num ./ den;
81 case {"P"}
82 num = e(j).*A - e(i).* B;
83 den = e(j).*A + e(i).* B;
84 r_ij = num ./ den;
85 end
86 end
```
#### **3. Fitting the experimental data to the model**

Data can be acquired using a standard scanning epifluorescence microscope with a spectrometer and photomultiplier tubes. It is generally most practical to save images as .lsm files which can be directly read into MATLAB [In our setup each file contains a scanned image of the sample at ten different wavelengths (e.g. bins centered at *λ* = 490nm, 500nm, 510nm, etc.)]. For dynamic studies images at different time intervals can also be stored in a single .lsm file. Many spectrometers come with MATLAB drivers which allow for acquisition to also be controlled easily through MATLAB. Else for real time (video type) imaging it is necessary to save the data incrementally into a directory accessible to MATLAB. The spectral profile of the fluorophores of interest on an uncoated substrate immersed in a similar refractive index medium as the actual samples should also be available. The area under this spectrum should be normalized and extrapolated to a discrete matrix [e.g. fitted to a high order polynomial and evaluated at the discrete wavelengths indexed by i], and saved as a 1 × 10 matrix Ilambda0(i).

Data analysis is performed on a pixel by pixel basis. For each pixel the spectrum area is normalized, fitted with a high order polynomial, evaluated at i discrete wavelengths, and written into a 1 × 10 matrix data(i). The ratio enhanobserved(i) = data(i)/Ilambda0(i) is then calculated. A 1 × *j* matrix zn(j) is defined which will represent the fraction of the emitters in the studied area at a distance corresponding to the index j from the interface. The task is now to solve enhan\*zn-enhanobserved=0, which is in principle straight forward using a function such as e.g. linsolve. The result is best presented in a plot of zn as a function of z for each or a collection of pixels, which will show the distribution of fluorophore distances from the substrate at this particular lateral position. Depending on the sample studied and the information of interest this may be fitted with a suitable distribution function to obtain an average and spread in distances of emitters from the substrate within the analyzed region.

## **4. Optimizing the model and code**

In many cases the initial fluorophore-concentration, metamaterial or coupling parameters may not be the best or correct parameters for the fit. This is usually evident when the results show an unexpected distribution for zn. Distances less than 10nm should be ignored as the approximations of the models (dipole approximation and local dielectric functions) are no longer valid at such small distances. The most significant source of uncertainty is the predefined and assumed constant value of the fluorophore concentration rho. To account for this one may redefine rho as being proportional to the integrated area of zn over z, and repeat the calculation of zn. Also to account for an axial variation in rho - as would be relevant for samples where the fluorophore density varies strongly with distance z - one may define rho(j) = rho\*(a + zn(j)/norm(zn)), where a is normalization constant chosen so that on summation over the distances (*i*) one recovers the total intensity. In this case all subsequent calculations have to explicitly account for the additional distance dependence in *rho (j)*, which significantly increases the computational time. This may be efficiently achieved by defining rho(j) in a smaller parameter space indexed by n that corresponds to dividing rho into components where rho(j) is constant and non-zero only over a specified region. Once zn is calculated in this manner the optimization may be repeated for the new rho until the the calculation of the distribution of zn is stable and smooth or varies on the expected scale.

### **5. Example**

10 Will-be-set-by-IN-TECH

Data can be acquired using a standard scanning epifluorescence microscope with a spectrometer and photomultiplier tubes. It is generally most practical to save images as .lsm files which can be directly read into MATLAB [In our setup each file contains a scanned image of the sample at ten different wavelengths (e.g. bins centered at *λ* = 490nm, 500nm, 510nm, etc.)]. For dynamic studies images at different time intervals can also be stored in a single .lsm file. Many spectrometers come with MATLAB drivers which allow for acquisition to also be controlled easily through MATLAB. Else for real time (video type) imaging it is necessary to save the data incrementally into a directory accessible to MATLAB. The spectral profile of the fluorophores of interest on an uncoated substrate immersed in a similar refractive index medium as the actual samples should also be available. The area under this spectrum should be normalized and extrapolated to a discrete matrix [e.g. fitted to a high order

52 exp\_term = exp(2i.\*k(1).\*d(3).\*sqrt(e(3).\*mu(3).\*(e(1).\*mu(1)).^-1-u.^2));

53

54 switch(polariz) 55 case {"P"}

61 case {"S"}

67 end 68 end 69

72

75

56 r\_23\_P = r\_ij\_PS(u,2,3,"P"); 57 r\_34\_P = r\_ij\_PS(u,3,4,"P");

62 r\_23\_S = r\_ij\_PS(u,2,3,"S"); 63 r\_34\_S = r\_ij\_PS(u,3,4,"S");

60 r\_234 = num ./ den;

66 r\_234 = num ./ den;

71 global e; global mu;

73 A = (1 - u.^2).^0.5;

80 r\_ij = num ./ den;

84 r\_ij = num ./ den;

76 switch(polariz) 77 case {"S"}

81 case {"P"}

85 end 86 end

58 num = r\_23\_P + r\_34\_P.\*exp\_term; 59 den = 1 + r\_23\_P .\* r\_34\_P.\*exp\_term;

64 num = r\_23\_S + r\_34\_S.\*exp\_term; 65 den = 1 + r\_23\_S .\* r\_34\_S.\*exp\_term;

74 B = ( e(j).\*mu(j).\* (e(i).\*mu(i)).^-1 - u.^2).^0.5;

70 function r\_ij = r\_ij\_PS(u,i,j,polariz)

78 num = mu(j).\*A - mu(i).\* B; 79 den = mu(j).\*A + mu(i).\* B;

82 num = e(j).\*A - e(i).\* B; 83 den = e(j).\*A + e(i).\* B;

**3. Fitting the experimental data to the model**

Whilst elaborate custom substrates may be designed with effective parameters which yield the desired cut-off energy and dispersion relation, it is often possible to use combinations of no more than 3-4 layers with no lateral structuring to achieve satisfactory results. Here we give an example where we use a sample consisting of a thin silver film2 (15nm) deposited on a thick quartz substrate (*ε* = 2.13), which is subsequently coated with a high permittivity dielectric film (*ε* = 4.15). The sample is immersed in a medium with a refractive index larger than that of

2 Smoothness of the silver is guaranteed by using a Germanium wetting layer and is less than 0.4nm RMS as measured by Atomic Force Microscopy.

**6. Modifications & extensions**

adjusted thresholding values.

A natural extension to the presented nano-sectioning method would be to combine it with a different technique which allows for lateral (xy) superresolution. One such technique that is relatively straight forward to implement is FPALM which was already alluded to above. In FPALM photoactivatable fluorophores are stochastically turned on by a weak activations source so that individual molecules separated by distances larger than the diffraction limit can consecutively be localized provided one has knowledge of their point spread function and the instrument and detector response. Numerous algorithms exist for such 2D localization purposes most of which rely on fitting a normal 2D (Gaussian) intensity distribution to a region of interest surrounding each molecule. Several such routines have been made available in the form of MATLAB codes. For densely labelled samples it is often desirable to be able to fit multiple overlapping Gaussians so that the isolated single-molecule activation/excitation condition can be slightly relaxed. This is particularly useful where even a very weak activation signal is sufficient to activate a significant portion of fluorophores - which we have found to be the case for photoactivatable fluorophores near the metallic structures of interest. A MATLAB code capable of fitting "overlapping intensity distributions" may be written employing a Levenberg-Marquardt based multi Gaussian least squares optimization algorithm. The result of a multiple Gaussian fit is demonstrated in figure 3 for the distribution from two nearby molecules. We note that such algorithms are also conducive to parallel processing and thus in principle fast dynamic studies. Considering the fact that typical amount of data needed to be processed for each microscopy experiment is in the order of several hundreds of megabytes; the possibility to use acceleration power of modern GPUs to enable near real-time observation of reconstruction is welcomed. At this moment MATLAB allows us to develop our implementations for several image processing and optimization routines and possibly validate the results from our experimental parallel implementations. In following text we will focus on the use of Image processing and Optimization toolboxes to outline our framework for detection regions of interest around sparsely activated fluorophores. We expect that reader is familiar with basics of image processing and numerical optimization methods and so we only propose possible solution for stated problem. For further details about mathematical

375 Using MATLAB to Achieve Nanoscale Optical Sectioning

in the Vicinity of Metamaterial Substrates by Simulating Emitter-Substrate Interactions

background regarding image analysis we suggest Sonka et al. (2001).

Current work is underway in our laboratory on combining this with the optical nanosectioning technique outlined above. A current obstacle appears to be the difficulty in rapidly photo-bleaching of the emitters - which is an essential feature of performing fast FPALM - near the substrates. This is a well known feature of fluorophores near metallic substrates and is understood to be due to the relative suppression of the triplet decay path, which is associated with the majority of the photobleaching. FPALM based nanosectioning as mentioned above is however still possible but at the price of lowered signal to noise ratio, and with properly designed image segmentation we are able extract regions aimed for further processing. The code shown below represents the necessary steps of the image processing required for detection of relevant fluorophores. This is based on dilation functions from the mathematical morphology toolset and subsequent segmentation with experimentally

the quartz (*ε* = 2.4). The cut-off energy can be shown to occur within the emission spectrum of Alexa488 - a standard fluorescent marker with an emission spectrum peaking at a wavelength of around 528nm and extending all the way up to >550nm. To demonstrate the accuracy of the technique, fibroblast cells were grown on the laminin (thickness < 2nm, *ε* ≈ 2.15) coated substrate and a particular protein found at adhesion sites close to the substrate called *Paxillin* Schaller (2001) was stained with Alexa488. The protein Paxillin has recently been shown to be separated from the substrate matrix by a distance of around 20-40nm using a 3D version of a technique based on Fluorescence Photoactivatable Localization Microscopy (FPALM) - known as interference-PALM or iPALM Kanchanawong (2010). In our study similar cells were used and the separation from the substrate surface can be expected to be comparable. In figure 2 we show a confocal fluorescence image of the cell we subsequently analyze (excitation wavelength = 488nm, emission filter ≈500-550nm). We focus on a region marked by the red box. A spectral analysis of the signal in this region shows an enhancement at longer wavelengths as expected. Using the procedure outlined in the previous sections we are able to calculate zn for the pixels in this box over the range z = 1...130. The fluorophore concentration (rho) was initially taken to be constant up to 2.4*μ*m and the presented histogram in figure 2 is the result following four iterations of redefining rho according to the calculated zn in a n=5 parameter space (see previous section). The size of the bin width for the spectral analysis was 9.7nm. As can be seen in figure 2 where the distance distribution (zn is plotted as a function of z), the fluorophores and hence the Paxillin are indeed localized at around the expected distance from the substrate (20-40nm). We note that compared to the iPALM setup used for these studies our technique is, once optimized, very simple and straight forward to implement. The analysis is sufficiently fast to allow for fast real-time studies (which is often impossible with e.g. iPALM) and also more versatile since there is no requirement for special photoactivatable fluorophorescent markers. We are currently using it to study the dynamics of other labeled proteins found in the vicinity of cell-substrate adhesion sites including the zinc-binding phosphoprotein *Zyxin* and the microfilament associated protein *VASP*.

Fig. 2. Left: Axial distance profile of the protein *paxillin* (zn versus z) as inferred from Δ*λ* =9.7nm resolution spectrum (*λ* =480→800nm) of pixels at adhesion site. Right: image of cell. Red square shows pixels over which the spectrum was analyzed to obtain image on the left. [NIH 3T3 fibroblast with Alexa488 stained paxillin, scanning confocal microscope image at *λ* =519(±9.7nm), with 1.4 NA objective and *λ* =465nm excitation. Scale-bar = 10*μ*m].

#### **6. Modifications & extensions**

12 Will-be-set-by-IN-TECH

the quartz (*ε* = 2.4). The cut-off energy can be shown to occur within the emission spectrum of Alexa488 - a standard fluorescent marker with an emission spectrum peaking at a wavelength of around 528nm and extending all the way up to >550nm. To demonstrate the accuracy of the technique, fibroblast cells were grown on the laminin (thickness < 2nm, *ε* ≈ 2.15) coated substrate and a particular protein found at adhesion sites close to the substrate called *Paxillin* Schaller (2001) was stained with Alexa488. The protein Paxillin has recently been shown to be separated from the substrate matrix by a distance of around 20-40nm using a 3D version of a technique based on Fluorescence Photoactivatable Localization Microscopy (FPALM) - known as interference-PALM or iPALM Kanchanawong (2010). In our study similar cells were used and the separation from the substrate surface can be expected to be comparable. In figure 2 we show a confocal fluorescence image of the cell we subsequently analyze (excitation wavelength = 488nm, emission filter ≈500-550nm). We focus on a region marked by the red box. A spectral analysis of the signal in this region shows an enhancement at longer wavelengths as expected. Using the procedure outlined in the previous sections we are able to calculate zn for the pixels in this box over the range z = 1...130. The fluorophore concentration (rho) was initially taken to be constant up to 2.4*μ*m and the presented histogram in figure 2 is the result following four iterations of redefining rho according to the calculated zn in a n=5 parameter space (see previous section). The size of the bin width for the spectral analysis was 9.7nm. As can be seen in figure 2 where the distance distribution (zn is plotted as a function of z), the fluorophores and hence the Paxillin are indeed localized at around the expected distance from the substrate (20-40nm). We note that compared to the iPALM setup used for these studies our technique is, once optimized, very simple and straight forward to implement. The analysis is sufficiently fast to allow for fast real-time studies (which is often impossible with e.g. iPALM) and also more versatile since there is no requirement for special photoactivatable fluorophorescent markers. We are currently using it to study the dynamics of other labeled proteins found in the vicinity of cell-substrate adhesion sites including the zinc-binding phosphoprotein *Zyxin*

and the microfilament associated protein *VASP*.

Fig. 2. Left: Axial distance profile of the protein *paxillin* (zn versus z) as inferred from

Δ*λ* =9.7nm resolution spectrum (*λ* =480→800nm) of pixels at adhesion site. Right: image of cell. Red square shows pixels over which the spectrum was analyzed to obtain image on the left. [NIH 3T3 fibroblast with Alexa488 stained paxillin, scanning confocal microscope image at *λ* =519(±9.7nm), with 1.4 NA objective and *λ* =465nm excitation. Scale-bar = 10*μ*m].

A natural extension to the presented nano-sectioning method would be to combine it with a different technique which allows for lateral (xy) superresolution. One such technique that is relatively straight forward to implement is FPALM which was already alluded to above. In FPALM photoactivatable fluorophores are stochastically turned on by a weak activations source so that individual molecules separated by distances larger than the diffraction limit can consecutively be localized provided one has knowledge of their point spread function and the instrument and detector response. Numerous algorithms exist for such 2D localization purposes most of which rely on fitting a normal 2D (Gaussian) intensity distribution to a region of interest surrounding each molecule. Several such routines have been made available in the form of MATLAB codes. For densely labelled samples it is often desirable to be able to fit multiple overlapping Gaussians so that the isolated single-molecule activation/excitation condition can be slightly relaxed. This is particularly useful where even a very weak activation signal is sufficient to activate a significant portion of fluorophores - which we have found to be the case for photoactivatable fluorophores near the metallic structures of interest. A MATLAB code capable of fitting "overlapping intensity distributions" may be written employing a Levenberg-Marquardt based multi Gaussian least squares optimization algorithm. The result of a multiple Gaussian fit is demonstrated in figure 3 for the distribution from two nearby molecules. We note that such algorithms are also conducive to parallel processing and thus in principle fast dynamic studies. Considering the fact that typical amount of data needed to be processed for each microscopy experiment is in the order of several hundreds of megabytes; the possibility to use acceleration power of modern GPUs to enable near real-time observation of reconstruction is welcomed. At this moment MATLAB allows us to develop our implementations for several image processing and optimization routines and possibly validate the results from our experimental parallel implementations. In following text we will focus on the use of Image processing and Optimization toolboxes to outline our framework for detection regions of interest around sparsely activated fluorophores. We expect that reader is familiar with basics of image processing and numerical optimization methods and so we only propose possible solution for stated problem. For further details about mathematical background regarding image analysis we suggest Sonka et al. (2001).

Current work is underway in our laboratory on combining this with the optical nanosectioning technique outlined above. A current obstacle appears to be the difficulty in rapidly photo-bleaching of the emitters - which is an essential feature of performing fast FPALM - near the substrates. This is a well known feature of fluorophores near metallic substrates and is understood to be due to the relative suppression of the triplet decay path, which is associated with the majority of the photobleaching. FPALM based nanosectioning as mentioned above is however still possible but at the price of lowered signal to noise ratio, and with properly designed image segmentation we are able extract regions aimed for further processing. The code shown below represents the necessary steps of the image processing required for detection of relevant fluorophores. This is based on dilation functions from the mathematical morphology toolset and subsequent segmentation with experimentally adjusted thresholding values.

**Listing 4. Fluorophore segmentation based on thresholding and morphology**

in the Vicinity of Metamaterial Substrates by Simulating Emitter-Substrate Interactions

[pos] = find\_local\_maxima(img, 0.7);% -- detect fluorophore positions

for k=1:1: length(x) % -- iterate over all fluorophore positions

 % -- rescale the window with nearest neighborhood interpolation % & prepare the meshgrid of the same size for Gaussian 2d fitting

% -- optimize parameters based on image data and gaussian model

% -- calculate goodness of fit for possible discarding of

 % -- calculate distance measure of model and image data abs\_dif\_fit\_values = [abs\_dif\_fit\_values; sum(abs(D(:))]; fprintf("Absolute Difference of Fit : %.3f \n", sum(abs(D(:)) );

[F,Yfit, D] = gaussian\_fit\_fun(params\_fit, px, py, ups\_win);

[win,corner\_x,corner\_y] = getROI(img,[y(k),x(k)],win\_width,win\_height);

ups\_win = imresize(win,[upscale\_factor\*Mw,upscale\_factor\*Nw],"nearest");

 % params\_fit = lsqnonlin(@(x) gaussian\_fit\_fun(x, px, py,window),param0); params\_fit = lev\_marq\_gaussian2d(ups\_win,px,py,[A,x0,sx,y0,sy 0]', ...

[x,y] = remove\_bound(x,y,win\_W,win\_H,M,N);% -- remove spots near the boundary

377 Using MATLAB to Achieve Nanoscale Optical Sectioning

win\_W = 8; % -- moving window size

stack = lsm\_read("1.lsm"); % -- load LSM stack

% -- find relevant positions of fluorophores

img = im2double(img); % -- change to 0-1 scale

[y,x] = find(pos == 1); % -- get their coordinates

upscaled\_window = zeros(M\*upscale\_factor,N\*upscale\_factor);

% -- get corners of the window from low-res image

% -- prepare guess estimates for Gaussian fitting

 upscale\_factor = 5; abs\_dif\_fit\_values = [];

img = stack(1).data;

% -- preallocated upscaled ROI

% -- where molecule was detected

[Mw,Nw] = size(win);

A = 10; z0 = 0;

 [ny,nx] = size(ups\_win); [px,py] = meshgrid(1:nx,1:ny);

 x0 = ny / 2; y0 = nx / 2; sx = 0.2 \* upscale\_factor; sy = 0.2 \* upscale\_factor; param0 = [A, x0, sx, y0, sy, z0];

1500,0.01,1e-6,0.01);

% -- incorrect values

[M,N] = size(img);

win\_H = 8;

end

**Listing 2. Fluorophore segmentation based on thresholding and morphology**

Depending on the sample studied, one may observe an uneven background signal which systematically introduce errors into the model fitting which may prove detrimental. For these reasons we estimate the background intensity map based on uniform filtering over large region, e.g.:

**Listing 3. Fluorophore segmentation with background correction**

```
1 corrected_bg = img_dbl - imfilter(img_dbl,fspecial("average",50));
2 img_corr_bg_closed_bw = im2bw(corrected_bg,0.60);
3 img_corr_bg_closed_bw_dil = imdilate(img_corr_bg_closed_bw,strel("disk",2));
```
Fig. 3. a. sample of activated fluorophores with noisy background b. segmented fluorophores after background correction c. overlay of fluorophores positions and corrected background

The main algorithm framework for detection and precise position localization of activated fluorophores consists of two parts. Lines (1)-(9) set the size of the region of interest around the detected fluorophore candidates, *upscalef actor* determine up to a certain level precision for the fitting (with higher values at the expense of increased computational resources). The following lines are for loading of the acquired data-stack and the rescaling to the MATLAB native double format (images with scales from 0 to 1). The subsequent section is based on the code presented in listings 4 or 3 to detect local maxima and translate them to (x,y) coordinates. Here the for loop iterates over the range of all (x,y) fluorophore coordinates, and the code is built from 2 parts, extraction of the image region of interest around the position of the fluorophore and fitting to the image data model (a Gaussian function based on a nonlinear iterative least-squares fitting). For the main algorithm of Levenberg-Marquard, interested readers are referred to Marquardt (1963)

Will-be-set-by-IN-TECH

Depending on the sample studied, one may observe an uneven background signal which systematically introduce errors into the model fitting which may prove detrimental. For these reasons we estimate the background intensity map based on uniform filtering over large

corrected\_bg = img\_dbl - imfilter(img\_dbl,fspecial("average",50));

img\_corr\_bg\_closed\_bw\_dil = imdilate(img\_corr\_bg\_closed\_bw,strel("disk",2));

Fig. 3. a. sample of activated fluorophores with noisy background b. segmented fluorophores after background correction c. overlay of fluorophores positions and corrected background The main algorithm framework for detection and precise position localization of activated fluorophores consists of two parts. Lines (1)-(9) set the size of the region of interest around the detected fluorophore candidates, *upscalef actor* determine up to a certain level precision for the fitting (with higher values at the expense of increased computational resources). The following lines are for loading of the acquired data-stack and the rescaling to the MATLAB native double format (images with scales from 0 to 1). The subsequent section is based on the code presented in listings 4 or 3 to detect local maxima and translate them to (x,y) coordinates. Here the for loop iterates over the range of all (x,y) fluorophore coordinates, and the code is built from 2 parts, extraction of the image region of interest around the position of the fluorophore and fitting to the image data model (a Gaussian function based on a nonlinear iterative least-squares fitting). For the main algorithm of Levenberg-Marquard, interested

**Listing 2. Fluorophore segmentation based on thresholding and morphology**

**Listing 3. Fluorophore segmentation with background correction**

img\_corr\_bg\_closed\_bw = im2bw(corrected\_bg,0.60);

 img = lsm\_read(fileName); img\_dbl = im2double(img); se = strel("disk" , 5); img\_dil = imdilate(img\_dbl,se); img\_dil\_bw = im2bw(img\_dil,threshold);

readers are referred to Marquardt (1963)

region, e.g.:

**Listing 4. Fluorophore segmentation based on thresholding and morphology**

```
1 win_W = 8; % -- moving window size
2 win_H = 8;
3 upscale_factor = 5;
4 abs_dif_fit_values = [];
5
6 stack = lsm_read("1.lsm"); % -- load LSM stack
7 img = stack(1).data;
8 img = im2double(img); % -- change to 0-1 scale
9 [M,N] = size(img);
10
11 % -- find relevant positions of fluorophores
12 [pos] = find_local_maxima(img, 0.7);% -- detect fluorophore positions
13 [y,x] = find(pos == 1); % -- get their coordinates
14 [x,y] = remove_bound(x,y,win_W,win_H,M,N);% -- remove spots near the boundary
15
16 % -- preallocated upscaled ROI
17 upscaled_window = zeros(M*upscale_factor,N*upscale_factor);
18
19 for k=1:1: length(x) % -- iterate over all fluorophore positions
20 % -- get corners of the window from low-res image
21 % -- where molecule was detected
22 [win,corner_x,corner_y] = getROI(img,[y(k),x(k)],win_width,win_height);
23
24 % -- rescale the window with nearest neighborhood interpolation
25 % & prepare the meshgrid of the same size for Gaussian 2d fitting
26 [Mw,Nw] = size(win);
27 ups_win = imresize(win,[upscale_factor*Mw,upscale_factor*Nw],"nearest");
28 [ny,nx] = size(ups_win);
29 [px,py] = meshgrid(1:nx,1:ny);
30
31 % -- prepare guess estimates for Gaussian fitting
32 A = 10; z0 = 0;
33 x0 = ny / 2; y0 = nx / 2;
34 sx = 0.2 * upscale_factor;
35 sy = 0.2 * upscale_factor;
36 param0 = [A, x0, sx, y0, sy, z0];
37
38 % -- optimize parameters based on image data and gaussian model
39 % params_fit = lsqnonlin(@(x) gaussian_fit_fun(x, px, py,window),param0);
40 params_fit = lev_marq_gaussian2d(ups_win,px,py,[A,x0,sx,y0,sy 0]', ...
41 1500,0.01,1e-6,0.01);
42 % -- calculate goodness of fit for possible discarding of
43 % -- incorrect values
44 [F,Yfit, D] = gaussian_fit_fun(params_fit, px, py, ups_win);
45
46 % -- calculate distance measure of model and image data
47 abs_dif_fit_values = [abs_dif_fit_values; sum(abs(D(:))];
48 fprintf("Absolute Difference of Fit : %.3f \n", sum(abs(D(:)) );
49 end
```
most software and self written codes, the intuitive matrix approach of MATLAB makes it well suited for the required tasks. Certain toolboxes may prove useful,3 however no toolboxes are essential or required. The presented technique is likely to find a range of applications in the lifesciences related to studying membrane traffic or indeed any events that can be made to occur near a coated substrate and nanoscale axial dynamics are of interest. In practice a fundamental limitations of the technique is the rate and accuracy with which the spectral information over a give region can be acquired. Whilst the use of multiple photodetectors combined with fast scanning of the sample may be the most accurate and best suited for not so bright samples, one can in principle perform the spectral separation optically using dichroic mirrors & filters and an array of sensitive CCD cameras (or several with split chips). The latter would remove the need for scanning, with the imaging rate now limited primarily by photon

379 Using MATLAB to Achieve Nanoscale Optical Sectioning

in the Vicinity of Metamaterial Substrates by Simulating Emitter-Substrate Interactions

D. Axelrod (1981). Cell-substrate contacts illuminated by total internal reflection fluorescence.

P.A. Belov and Y. Hao (2006). Subwavelength imaging at optical frequencies using a

D.J. Bergman and M.I. Stockman (2003). Surface Plasmon Amplification by Stimulated

Betzig, E. et al. (1991). Breaking the Diffraction Barrier: Optical Microscopy on a Nanometric

M. Born and E. Wolf (1997). *Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light*, Cambridge University Press, Cambridge, UK. J.J. Burke and G.I. Stegeman (1986). Surface-polariton-like waves guided by thin, lossy metal

R.R. Chance, A. Prock and R. Silbey (1978). Molecular fluorescence and energy transfer near

C. Cremer and T. Cremer (1978). Considerations on a laser-scanning-microscope with high resolution and depth of field. *Microscopica Acta*, Vol.81, No.1, pp.31-44. S.A. Darmanyan and A.V. Zayats (2003). Light tunneling via resonant surface plasmon

J.A. Dionne et al. (2008). Are negative index materials achievable with surface plasmon

V.P. Drachev (2008). The Ag dielectric function in plasmonic metamaterials. *Opt. Exp.*, Vol.16,

3 In particular the *Symbolic Math* and *Curve-fitting* toolboxes for modeling the structures and fitting the spectrums, and the *Instrument Control* and *Image Acquisition* toolboxes for controlling optical

films: An analytical study. *Phys. Rev. B*, Vol.67, pp.035424-035430.

polariton states and the enhanced transmission of periodically nanostructured metal

waveguides? A case study of three plasmonic geometries. *Opt. Exp.*, Vol.16, No.23,

in the canalization regime. *Phys. Rev. B*, Vol.73, pp.113110-113113.

Nanosystems. *Phys. Rev. Lett.*, Vol.90, pp.027402-027405.

Scale. *Science*, Vol.251, No.5000, pp.1468-1470.

films. *Phys. Rev. B*, Vol.33, pp.51865201.

pp.19001-19016.

No.2, pp.1186- 1195.

components & shutters and data acquisition in MATLAB.

interfaces. *Adv. Chem. Phys.*, Vol.37, pp.165.

transmission device formed by a periodic layered metal-dielectric structure operating

Emission of Radiation: Quantum Generation of Coherent Surface Plasmons in

statistics.

**8. References**

*Journ Cell Biol*, Vol.89, pp.141-145.

Fig. 4. a. Extracted ROI from raw microscopy data b. Initial shape of Gaussian model c. Fitted model

Fig. 5. Fitting two overlapping Gaussians to data using a Levenberg-Marquardt based routine. Left: raw data. Center: Gaussian Fits. Right: Wellness of fit as a function of number of iterations.

#### **7. Concluding remarks**

In summary we have outlined the data analysis protocol for a novel fluorescent imaging technique capable of performing sub-wavelength scale optical sectioning in the vicinity of metamaterial coated substrates - MeSuMo fluorescence microscopy. An outline of the approach one would take for analyzing the response properties of particular materials and how to go about fitting data to the model is explained from a perspective that may be readily implemented into MATLAB routines. Whilst the calculation can be performed by most software and self written codes, the intuitive matrix approach of MATLAB makes it well suited for the required tasks. Certain toolboxes may prove useful,3 however no toolboxes are essential or required. The presented technique is likely to find a range of applications in the lifesciences related to studying membrane traffic or indeed any events that can be made to occur near a coated substrate and nanoscale axial dynamics are of interest. In practice a fundamental limitations of the technique is the rate and accuracy with which the spectral information over a give region can be acquired. Whilst the use of multiple photodetectors combined with fast scanning of the sample may be the most accurate and best suited for not so bright samples, one can in principle perform the spectral separation optically using dichroic mirrors & filters and an array of sensitive CCD cameras (or several with split chips). The latter would remove the need for scanning, with the imaging rate now limited primarily by photon statistics.

## **8. References**

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

Fig. 4. a. Extracted ROI from raw microscopy data b. Initial shape of Gaussian model c.

Fig. 5. Fitting two overlapping Gaussians to data using a Levenberg-Marquardt based routine. Left: raw data. Center: Gaussian Fits. Right: Wellness of fit as a function of number

In summary we have outlined the data analysis protocol for a novel fluorescent imaging technique capable of performing sub-wavelength scale optical sectioning in the vicinity of metamaterial coated substrates - MeSuMo fluorescence microscopy. An outline of the approach one would take for analyzing the response properties of particular materials and how to go about fitting data to the model is explained from a perspective that may be readily implemented into MATLAB routines. Whilst the calculation can be performed by

Fitted model

of iterations.

**7. Concluding remarks**

3 In particular the *Symbolic Math* and *Curve-fitting* toolboxes for modeling the structures and fitting the spectrums, and the *Instrument Control* and *Image Acquisition* toolboxes for controlling optical components & shutters and data acquisition in MATLAB.

**16**

*1Brazil 2Ireland*

**A Methodology and Tool to**

*1Univ Estadual Paulista, UNESP; 2University of Limerick, UL;*

**Translate MATLAB**®**/Simulink**® **Models of**

Alexandre César Rodrigues da Silva1 and Ian Andrew Grout2

Nowadays, due to the increasing complexity of circuits and systems incorporating mixedsignal components, many designers recognise the advantages of, and in some cases the necessity for, mixed-signal modelling and simulation of the circuits and systems as a whole prior to, during and post production of the circuit or system. This is particularly so within the microelectronics (semiconductor) arena. This, coupled with the emergence of a range of description languages that find use in the modelling and analysis via simulation of microelectronic circuit designs, means the ability now exists to design complex signal processing functions in both the digital and analogue domains. Modelling and simulation has increased the designer efficiency and ability to develop increasingly complex and useful electronic circuits, particularly at the integrated circuit (IC) level. Each type of circuit considered has its own particular requirements in terms of design performance (specifications), design methodologies required to realise a design, modelling and simulation

This chapter will consider design aspects relating specifically to data converters, in particular the digital-to-analogue converter (DAC). In many electronic systems, the data converter (both the analogue-to-digital [A/D] and digital-to-analogue [D/A] converter) provides an important electronic function. The developed converter design will be provided to the end-user as either a discrete integrated circuit or as a macro cell within a larger

The data converter provides the primary interface block between the analogue 'real' world and the digital signal processing circuitry typically found in many control, instrumentation and, generically, 'mechatronic' systems. The converter design performance is continually increasing to meet the changing needs of the end-user applications, in particular: an increase in resolution (number of bits); an increase in speed of operation - driven primarily by the communications applications; the operation of the circuit on lower power supply voltages and with lower power consumption - driven primarily by the need for portable electronic

**1. Introduction**

design, see Fig. 1.

circuit applications.

toolsets utilised, designer experience and skills set.

**Mixed-Signal Circuits to VHDL-AMS**

## **A Methodology and Tool to Translate MATLAB**®**/Simulink**® **Models of Mixed-Signal Circuits to VHDL-AMS**

Alexandre César Rodrigues da Silva1 and Ian Andrew Grout2 *1Univ Estadual Paulista, UNESP; 2University of Limerick, UL; 1Brazil 2Ireland*

## **1. Introduction**

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

380 Engineering Education and Research Using MATLAB

K. Elsayad and K.G. Heinze (2010). Temperature dependence of the near-field superlensing

K. Elsayad and K.G. Heinze (2010). Multifrequency parallelized near-field optical imaging with anistropic metal-dielectric stacks. *Phys. Rev. A*, Vol.81, pp.053840-053845.

P. Kanchanawong et al. (2010). Nanoscale architecture of integrin-based cell adhesions. *Nat.*,

D. Marquardt (1963). An Algorithm for Least Squares Estimation on Nonlinear Parameters.

S.A. Ramakrishna et al. (2003). Imaging the near field. *Journ. Mod. Opt.*, Vol.50, No.9,

M.D. Schaller (2001). Paxillin: a focal adhesion-associated adaptor protein. *Oncogene*, Vol.20,

M. Sonka et al. (2001). Image Processing: Analysis and Machine Vision, ISBN: 049508252X.

S.W. Hell (2007). Far-Field Optical Nanoscopy. *Science*, Vol.316, No.5828, pp.1153-1158. Z. Jacob et al. (2010). Engineering photonic density of states using metamaterials. *App. Phys.*

C.Kittel (1995). *Introduction to Solid State Physics*, Wiley, ISBN<PHONE_NUMBER>.

N.I. Zheludev et al. (2008). Lasing spaser. *Nat. Phot.*, Vol.2, pp.351-354.

77573L, San Diego, USA, Aug. 2010.

SIAMJ. J. APPL. MATH.11, pp. 431-441.

*B*, Vol.100, No.1, pp.215-218.

Vol.468, pp.580-584.

No.44 pp.6459-6472.

pp.1419-1430.

effect for single metal layers and stacked metal-dielectric films, *Proceedings of SPIE*,

Nowadays, due to the increasing complexity of circuits and systems incorporating mixedsignal components, many designers recognise the advantages of, and in some cases the necessity for, mixed-signal modelling and simulation of the circuits and systems as a whole prior to, during and post production of the circuit or system. This is particularly so within the microelectronics (semiconductor) arena. This, coupled with the emergence of a range of description languages that find use in the modelling and analysis via simulation of microelectronic circuit designs, means the ability now exists to design complex signal processing functions in both the digital and analogue domains. Modelling and simulation has increased the designer efficiency and ability to develop increasingly complex and useful electronic circuits, particularly at the integrated circuit (IC) level. Each type of circuit considered has its own particular requirements in terms of design performance (specifications), design methodologies required to realise a design, modelling and simulation toolsets utilised, designer experience and skills set.

This chapter will consider design aspects relating specifically to data converters, in particular the digital-to-analogue converter (DAC). In many electronic systems, the data converter (both the analogue-to-digital [A/D] and digital-to-analogue [D/A] converter) provides an important electronic function. The developed converter design will be provided to the end-user as either a discrete integrated circuit or as a macro cell within a larger design, see Fig. 1.

The data converter provides the primary interface block between the analogue 'real' world and the digital signal processing circuitry typically found in many control, instrumentation and, generically, 'mechatronic' systems. The converter design performance is continually increasing to meet the changing needs of the end-user applications, in particular: an increase in resolution (number of bits); an increase in speed of operation - driven primarily by the communications applications; the operation of the circuit on lower power supply voltages and with lower power consumption - driven primarily by the need for portable electronic circuit applications.

A Methodology and Tool to

Signal) language.

Translate MATLAB®/Simulink® Models of Mixed-Signal Circuits to VHDL-AMS 383

interaction between the various domains and determine that the behavioural models in the different domains can be created using the VHDL-AMS (VHDL Analogue and Mixed-

In Christen's work (Christen E. & Bakalar L., 1999), a revision of the VHDL-AMS language for analogue applications and mixed-signal is provided. The main elements of the language are studied. Doboli's work (Doboli A. & Vermuri R., 2003) and that of Pêcheaus (Pêcheaus et

Trofimov and Mosin (Travimov M. & Mosin S., 2004) note that one way of increasing the speed of undertaking project development is to use the top-down methodology. In analogue and mixed-signal circuit projects, the generation of analogue models is at a high level of abstraction. The use of models at high levels of abstraction permits a complete simulation of a mixed-signal device. However, an adequate tool must be used for modelling and

Grout and Keane (Grout I. A. & Keane K., 2000) show a prototype of one tool (software toolbox) that analyses and processes a Simulink® block diagram model producing a VHDL representation of the model. The derived VHDL model will be a combination of behavioural, RTL (Register Transfer Logic) and structural definitions mapped directly from the Simulink® model. This approach was considered to enable a user to develop and simulate a digital control algorithm using MATLAB® and once the initial algorithm

In a study undertaken the following year, Grout (Grout I. A., 2001) show the application of the developed tool to transform a process modelled into Simulink® in a described model in

Zorzi (Zorzi M. et al., 2004) describe a tool, called *I.M.A.Ge-AMS*, for writing new Spice models using the VHDL-AMS standard. The core of the *I.M.A.Ge-AMS* tool is a compiler derived from the VHDL-AMS compiler, developed for another EDA tool called *SANSA*. In 2004, a paper by Grout and O'Shea (Grout I. A. & O´Shea T., 2004) discusses the need to provide suitable provision for, and flexibility in the use of, modelling and simulation languages suitable for supporting the range of mixed-signal microelectronic circuit design, test and test development activities that are encountered in today's complex microelectronic products. In this work, emphasis is placed on language support for test development activities. The target area considered is Delta-Sigma (Δ-Σ) modulation for on-chip signal generation in the support of mixed-signal built-in self-test (MS-BIST). In the same year, Zorzi's work (Zorzi M. et al., 2004) shows the tool to integrate the VHDL-AMS with Spice. The developed tool allows the user to write new models using VHDL-AMS and automatically generates the embedded codes to be used by Spice. This approach shows the advantage of using a behavioural language with a simulator that is the best known of all

The following section will present the basic theory of digital-to-analogue converter operation and modelling with emphasis in the R-2R ladder networks. The mathematical

In this section, the concepts used in digital-to-analogue converter design modelling are reviewed. This type of circuit was selected because it is the most frequently encountered mixed-signal circuit. The data converter can be either an analogue-to-digital (A/D)

al., 2005) also deal with the manner of modelling using VHDL-AMS.

simulation, and there must be a great interaction between the project tools.

development was completed, to then convert this to synthesisable VHDL code.

VHDL. A closed-loop control system is used in this work.

analogue simulators at the present time.

model of R-2R ladder network is also presented.

**2. Digital-to-analogue converter modelling**

The focus area of this work was in the area of device modelling and simulation during the circuit design process, in particular, the conversion (translation) of model design descriptions between the different languages typically used by the design community, along with design synthesis requirements. The design and modelling of data converter algorithms was investigated using MATLAB®/Simulink® (The MathWorks Inc.) and then the conversion of the MATLAB®/Simulink® circuit model (of the data converter design) into a VHDL-AMS description was considered.

Fig. 1. Data converter designs embedded within a digital processor

A number of authors have already undertaken work in this area. In the next section, a summary of some of the prior work undertaken that have influenced the present work is presented.

## **1.1 Summary of related work**

Electronic devices become more complex each day. To deal with this increasing complexity, design synthesis methodologies, assisted by suitable software tools, have been developed in order to allow designers to initially model their circuit or system at different design abstraction levels during the various design cycles and then ultimately synthesis of some or the whole design into an electronic circuit implementation within a particular design project.

In the work of Edenfeld (Edenfeld D. et al., 2004), the future trend of semiconductor technology is presented, with the advent of systems with mixing signals (mixed-signal). He noted that even in highly abstract projects, it is possible to discover potential problems in the initial stage of development, therefore reducing the project time and decreasing the cost.

In MacMillen's work (MacMillen D. et al., 2000), a complete overview of the technologies, algorithms and methodologies that are used in the EDA tools (electronic design automation) and the economic impact of these technologies is undertaken. They consider the different steps of the project: simulation, verification, synthesis and test. Also discussed are the types of necessary toolsets to support a project environment.

The necessity of behavioural modelling within mixed technology systems, that is systems with electronic, electric and non-electric (mechanical, pneumatic, etc.) parts operating together, is described in Wilson's work (Wilson P. R. et al., 1997). They present the

The focus area of this work was in the area of device modelling and simulation during the circuit design process, in particular, the conversion (translation) of model design descriptions between the different languages typically used by the design community, along with design synthesis requirements. The design and modelling of data converter algorithms was investigated using MATLAB®/Simulink® (The MathWorks Inc.) and then the conversion of the MATLAB®/Simulink® circuit model (of the data converter design) into a VHDL-AMS

> Analogue output

A number of authors have already undertaken work in this area. In the next section, a summary of some of the prior work undertaken that have influenced the present work is

Electronic devices become more complex each day. To deal with this increasing complexity, design synthesis methodologies, assisted by suitable software tools, have been developed in order to allow designers to initially model their circuit or system at different design abstraction levels during the various design cycles and then ultimately synthesis of some or the whole design into an electronic circuit implementation within a particular design

In the work of Edenfeld (Edenfeld D. et al., 2004), the future trend of semiconductor technology is presented, with the advent of systems with mixing signals (mixed-signal). He noted that even in highly abstract projects, it is possible to discover potential problems in the initial stage of development, therefore reducing the project time and decreasing the cost. In MacMillen's work (MacMillen D. et al., 2000), a complete overview of the technologies, algorithms and methodologies that are used in the EDA tools (electronic design automation) and the economic impact of these technologies is undertaken. They consider the different steps of the project: simulation, verification, synthesis and test. Also discussed are the types

The necessity of behavioural modelling within mixed technology systems, that is systems with electronic, electric and non-electric (mechanical, pneumatic, etc.) parts operating together, is described in Wilson's work (Wilson P. R. et al., 1997). They present the

Processor: digital signal processing and data converter control

Data Data

Control

ADC DAC

Fig. 1. Data converter designs embedded within a digital processor

of necessary toolsets to support a project environment.

description was considered.

Analogue input

**1.1 Summary of related work**

presented.

project.

interaction between the various domains and determine that the behavioural models in the different domains can be created using the VHDL-AMS (VHDL Analogue and Mixed-Signal) language.

In Christen's work (Christen E. & Bakalar L., 1999), a revision of the VHDL-AMS language for analogue applications and mixed-signal is provided. The main elements of the language are studied. Doboli's work (Doboli A. & Vermuri R., 2003) and that of Pêcheaus (Pêcheaus et al., 2005) also deal with the manner of modelling using VHDL-AMS.

Trofimov and Mosin (Travimov M. & Mosin S., 2004) note that one way of increasing the speed of undertaking project development is to use the top-down methodology. In analogue and mixed-signal circuit projects, the generation of analogue models is at a high level of abstraction. The use of models at high levels of abstraction permits a complete simulation of a mixed-signal device. However, an adequate tool must be used for modelling and simulation, and there must be a great interaction between the project tools.

Grout and Keane (Grout I. A. & Keane K., 2000) show a prototype of one tool (software toolbox) that analyses and processes a Simulink® block diagram model producing a VHDL representation of the model. The derived VHDL model will be a combination of behavioural, RTL (Register Transfer Logic) and structural definitions mapped directly from the Simulink® model. This approach was considered to enable a user to develop and simulate a digital control algorithm using MATLAB® and once the initial algorithm development was completed, to then convert this to synthesisable VHDL code.

In a study undertaken the following year, Grout (Grout I. A., 2001) show the application of the developed tool to transform a process modelled into Simulink® in a described model in VHDL. A closed-loop control system is used in this work.

Zorzi (Zorzi M. et al., 2004) describe a tool, called *I.M.A.Ge-AMS*, for writing new Spice models using the VHDL-AMS standard. The core of the *I.M.A.Ge-AMS* tool is a compiler derived from the VHDL-AMS compiler, developed for another EDA tool called *SANSA*.

In 2004, a paper by Grout and O'Shea (Grout I. A. & O´Shea T., 2004) discusses the need to provide suitable provision for, and flexibility in the use of, modelling and simulation languages suitable for supporting the range of mixed-signal microelectronic circuit design, test and test development activities that are encountered in today's complex microelectronic products. In this work, emphasis is placed on language support for test development activities. The target area considered is Delta-Sigma (Δ-Σ) modulation for on-chip signal generation in the support of mixed-signal built-in self-test (MS-BIST). In the same year, Zorzi's work (Zorzi M. et al., 2004) shows the tool to integrate the VHDL-AMS with Spice. The developed tool allows the user to write new models using VHDL-AMS and automatically generates the embedded codes to be used by Spice. This approach shows the advantage of using a behavioural language with a simulator that is the best known of all analogue simulators at the present time.

The following section will present the basic theory of digital-to-analogue converter operation and modelling with emphasis in the R-2R ladder networks. The mathematical model of R-2R ladder network is also presented.

## **2. Digital-to-analogue converter modelling**

In this section, the concepts used in digital-to-analogue converter design modelling are reviewed. This type of circuit was selected because it is the most frequently encountered mixed-signal circuit. The data converter can be either an analogue-to-digital (A/D)

A Methodology and Tool to

**2.1 The R-2R ladder network**

between 0 volts and *Vref* volts.

critical to the function of the R-2R ladder.

Fig. 3. The R-2R ladder network with N bits

Translate MATLAB®/Simulink® Models of Mixed-Signal Circuits to VHDL-AMS 385

switch (*Sn*) controls a current source (*In*). When a switch is ON, its respective current source contributes its current to the output current (*Isum*), which becomes proportional to the

A popular approach used within the digital-to-analogue converter available in IC form is the R-2R current ladder. Only two resistor values are required, regardless of the number of bits used, which simplifies resistor trimming during the IC fabrication process. The principle of the R-2R ladder is that at each node the ladder is divided in two. The current divides equally because the resistor values are 2R for each leg. All the models used in this work use the R-2R current ladder network. The R-2R ladder network is described in the next section.

The R-2R ladder network provides for a simple and inexpensive way in which to perform digital-to-analogue conversion. Its popularity is due to the networks' inherent accuracy superiority (when compared with other models) and ease of manufacture. Fig. 3 presents a diagram of the basic R-2R ladder network with *N* bits. Note that the network consists of only two resistors values. One resistor has *R* value and the other one has *2R* value (twice the value of *R*) no matter how many bits make up the ladder. The particular value of *R* is not

The node labelled '*grd*' (at one end of the termination resistor) is connected to *ground*. The termination resistor assures that the Thevenin resistance of the network, as measured to ground looking toward the LSB (with all the bits grounded), is *R*. The Thevenin resistance of an R-2R ladder is always *R*, regardless of the number of bits in the ladder. From Fig. 3, it can be noted that the information is presented to the ladder as individual bits of a digital word switched between a reference voltage, labelled '*Vref*' and ground. The output voltage *Vout* is dependent on the number and location of the bits switched to *Vref* or ground. *Vout* will vary

If all inputs are connected to ground, 0 volts are produced at the output *Vout*. If all inputs are connected to *Vref*, the output voltage approaches *Vref*, and if some inputs are connected to ground and some to *Vref* then an output voltage between 0 volts and *Vref* appears at *Vout*. These inputs range from the most significant bit (MSB) to the least significant bit (LSB).

parallel input presented as the digital input (*digital data input*).

(converting signals from the analogue to the digital domain) or a digital-to-analogue (D/A) converter (converting signals from the digital to the analogue domain).

With the increasing requirements on the design (and testing) of the data converter, the development of electronic design automation (EDA) tools for supporting design and testing of A/D and D/A converters are becoming more important. In addition, with the advances obtained in the development and use of hardware description languages (HDLs), new EDA tools which integrate different environments are being created. This has permitted the migration from development of structural representation (logic synthesis) to behavioural specifications (higher levels abstraction). The new behavioural specification has removed many of the traditional constraints to design and has given the designer freedom of choice as to the best solution to any particular problem.

One of the basic ideas developed was to create the model using mathematical blocks exercising a top-down design methodology for both continuous time and discrete time in nature. With the Simulink® toolbox, these models are represented using a graphical block diagram format, in a visual form that can aid in supporting the designer's understanding of the problem.

Digital-to-analogue conversion is the process in which data in discrete values are converted to a continuous variable form (current or voltage). A characteristic feature in digital-toanalogue conversion is the relative value of bits. In a binary code, each more significant bit has twice the value of the previous bit. That is, in the representation 2n, where n = 0, 1, 2, 3, 4… etc., the n + 1th bit is in a more significant position than the bit in the nth position. An analogue voltage generated from a binary code will have components that are related by powers of 2.

In the present models for commercially available DACs, the output voltage of logic gates (output digital level) is used to activate electronic switches that apply zero voltage to the inputs of the digital-to-analogue converter for the logic 0 state and 5.0 V, or another reference voltage, for the logic 1 state. Fig. 2 shows the simplified circuit of the basic digitalto-analogue converter with 4 bits. It converts each bit of its inputs (digital data) that has a binary value of 1 to a current, and sums these currents into *Isum* at its analogue output. Each

Fig. 2. Simplified circuit diagram of a basic digital-to-analogue converter with 4 bits

switch (*Sn*) controls a current source (*In*). When a switch is ON, its respective current source contributes its current to the output current (*Isum*), which becomes proportional to the parallel input presented as the digital input (*digital data input*).

A popular approach used within the digital-to-analogue converter available in IC form is the R-2R current ladder. Only two resistor values are required, regardless of the number of bits used, which simplifies resistor trimming during the IC fabrication process. The principle of the R-2R ladder is that at each node the ladder is divided in two. The current divides equally because the resistor values are 2R for each leg. All the models used in this work use the R-2R current ladder network. The R-2R ladder network is described in the next section.

#### **2.1 The R-2R ladder network**

384 Engineering Education and Research Using MATLAB

(converting signals from the analogue to the digital domain) or a digital-to-analogue (D/A)

With the increasing requirements on the design (and testing) of the data converter, the development of electronic design automation (EDA) tools for supporting design and testing of A/D and D/A converters are becoming more important. In addition, with the advances obtained in the development and use of hardware description languages (HDLs), new EDA tools which integrate different environments are being created. This has permitted the migration from development of structural representation (logic synthesis) to behavioural specifications (higher levels abstraction). The new behavioural specification has removed many of the traditional constraints to design and has given the designer freedom of choice

One of the basic ideas developed was to create the model using mathematical blocks exercising a top-down design methodology for both continuous time and discrete time in nature. With the Simulink® toolbox, these models are represented using a graphical block diagram format, in a visual form that can aid in supporting the designer's understanding of

Digital-to-analogue conversion is the process in which data in discrete values are converted to a continuous variable form (current or voltage). A characteristic feature in digital-toanalogue conversion is the relative value of bits. In a binary code, each more significant bit has twice the value of the previous bit. That is, in the representation 2n, where n = 0, 1, 2, 3, 4… etc., the n + 1th bit is in a more significant position than the bit in the nth position. An analogue voltage generated from a binary code will have components that are related by

In the present models for commercially available DACs, the output voltage of logic gates (output digital level) is used to activate electronic switches that apply zero voltage to the inputs of the digital-to-analogue converter for the logic 0 state and 5.0 V, or another reference voltage, for the logic 1 state. Fig. 2 shows the simplified circuit of the basic digitalto-analogue converter with 4 bits. It converts each bit of its inputs (digital data) that has a binary value of 1 to a current, and sums these currents into *Isum* at its analogue output. Each

Fig. 2. Simplified circuit diagram of a basic digital-to-analogue converter with 4 bits

converter (converting signals from the digital to the analogue domain).

as to the best solution to any particular problem.

the problem.

powers of 2.

The R-2R ladder network provides for a simple and inexpensive way in which to perform digital-to-analogue conversion. Its popularity is due to the networks' inherent accuracy superiority (when compared with other models) and ease of manufacture. Fig. 3 presents a diagram of the basic R-2R ladder network with *N* bits. Note that the network consists of only two resistors values. One resistor has *R* value and the other one has *2R* value (twice the value of *R*) no matter how many bits make up the ladder. The particular value of *R* is not critical to the function of the R-2R ladder.

The node labelled '*grd*' (at one end of the termination resistor) is connected to *ground*. The termination resistor assures that the Thevenin resistance of the network, as measured to ground looking toward the LSB (with all the bits grounded), is *R*. The Thevenin resistance of an R-2R ladder is always *R*, regardless of the number of bits in the ladder. From Fig. 3, it can be noted that the information is presented to the ladder as individual bits of a digital word switched between a reference voltage, labelled '*Vref*' and ground. The output voltage *Vout* is dependent on the number and location of the bits switched to *Vref* or ground. *Vout* will vary between 0 volts and *Vref* volts.

Fig. 3. The R-2R ladder network with N bits

If all inputs are connected to ground, 0 volts are produced at the output *Vout*. If all inputs are connected to *Vref*, the output voltage approaches *Vref*, and if some inputs are connected to ground and some to *Vref* then an output voltage between 0 volts and *Vref* appears at *Vout*. These inputs range from the most significant bit (MSB) to the least significant bit (LSB).

A Methodology and Tool to

Inc. [c]).

Translate MATLAB®/Simulink® Models of Mixed-Signal Circuits to VHDL-AMS 387

MATLAB®/Simulink® or any other mathematical modelling environment. Note that this

In this work, models of commercially available digital-to-analogue converters were developed. All converters considered had a resolution of 8 bits and used the R-2R ladder approach to implement the conversion circuitry. The IC designs used were the DAC08 (Analog Devices Inc. [a]), AD7524 (Analog Devices Inc. [b]) and AD5450 (Analog Devices

Today, the designer of digital circuits and systems typically uses hardware description languages (HDLs) such as VHDL or Verilog®-HDL in order model the circuit or system at a high level of design abstraction. Although these HDLs provide language constructs for behavioural simulation, their language subsets are far too restrictive for system level design. On the other hand, the MATLAB®/Simulink® modelling and simulation environment provides a powerful high level mathematical modelling environment for systems that can be widely used for high level algorithm development. It is a high-performance language for technical computing. It has a large collection of computational algorithms ranging from elementary functions, like summation and multiplication, through to complex arithmetic

A hardware description language, such as VHDL, can readily be used to describe a digital circuit or system. However, this language has a very limited capability when dealing with complex mathematical operations. The use of an interface, however, between MATLAB®/Simulink® and certain HDLs can fill the gap between an HDL and a high level

Fig. 4. Mathematical approach for modelling the R-2R ladder using Simulink®

digital-to-analogue converters selected for this work are presented.

**3. Data converter models implemented in MATLAB®**

In the next section, models developed in MATLAB®/Simulink® that implemented the

A highly useful feature of Simulink® used in this work is the possibility of using it in the development of hybrid system models. In the DAC models developed, the use of *sum*, *multiplier*, *logic gates* and *shift register* components, along with components to generate input

**/Simulink®**

model is similar to R-2R ladder network presented in Fig. 3.

functions, such as fast Fourier transforms (FFTs).

mathematical modelling and simulation environment.

The MSB, when activated, causes the greatest change in the output voltage and the LSB, when activated, will cause the smallest change in the output voltage. The output voltage caused by connecting a particular bit to *Vref* with all other bits in grounded is *Vout= Vref / 2N*, where *N* is the bit number. Then for *N = 1*, *Vout = Vref / 2*, for *N = 2*, *Vout = Vref / 4*, etc. Table 1 shows the effect of individual bit locations to the *Nth* bit. Notice that since the bit *N =1* has the greatest effect on the output voltage, it is designated the MSB.

Table 1. Effect of individual bits on *Vout*

The R-2R ladder is a linear circuit. Then, by the principle of superposition, the output can be calculated by taking the sum of the effect of all bits connected to *Vref*. For example, if bits *1* and *3* are connected to *Vref* with all other inputs grounded, the output voltage is calculated by *Vout = (Vref /2) + (Vref /8)* which reduces to *Vout = (5.Vref) /8*.

As the R-2R ladder is a binary circuit, the effect of each successive bit approaching the LSB (less significant bit) is half that of the previous bit. If this sequence is extended to a ladder of infinite bits, the effect of the LSB on *Vout* approaches 0. Conversely, the full-scale output of the network approaches *Vref*.

#### **2.2 Mathematical model that represents the R-2R ladder**

As shown in Table 1, each bit contributes one small part of the reference voltage in the generation of the output voltage. The sum of the contribution of each individual bit then represents the analogue signal output. For example, in a DAC with 8 bits, the most significant bit contributes 0.5 of the reference voltage. The least significant bit contributes only 0.00390625 of the reference voltage. Each bit then contributes one specific gain value to the output voltage. The sum of all gains multiplied by the reference voltage then generates the final output voltage.

A simple scheme for analogue signal generation using only the mathematical approach was created for this work. As shown in Fig. 4, the mathematical approach of R-2R ladder consists of an arrangement of *gain*, *sum* and *multiplier* blocks. Each of the *gain* blocks has a weighting factor that is added when its input is supplied with a logic 1 and multiplied by reference voltage, so generating the final output voltage.

In the example presented in Fig. 4, using *Vref = 1V* when *B3 = B0 = 1* and *B2 = B1 = 0*, the output value is *0.5625* which represents the addition of each input. The advantage of the approach is that the DAC using R-2R ladder can be easily simulated using the

The MSB, when activated, causes the greatest change in the output voltage and the LSB, when activated, will cause the smallest change in the output voltage. The output voltage caused by connecting a particular bit to *Vref* with all other bits in grounded is *Vout= Vref / 2N*, where *N* is the bit number. Then for *N = 1*, *Vout = Vref / 2*, for *N = 2*, *Vout = Vref / 4*, etc. Table 1 shows the effect of individual bit locations to the *Nth* bit. Notice that since the bit

> Bit number Vout 1 (MSB) Vr/2 2 Vr/4 3 Vr/8 4 Vr/16 5 Vr/32 6 Vr/64 7 Vr/128 8 Vr/256 9 Vr/512 10 Vr/1024 n (LSB) Vr/2n

The R-2R ladder is a linear circuit. Then, by the principle of superposition, the output can be calculated by taking the sum of the effect of all bits connected to *Vref*. For example, if bits *1* and *3* are connected to *Vref* with all other inputs grounded, the output voltage is calculated

As the R-2R ladder is a binary circuit, the effect of each successive bit approaching the LSB (less significant bit) is half that of the previous bit. If this sequence is extended to a ladder of infinite bits, the effect of the LSB on *Vout* approaches 0. Conversely, the full-scale output of

As shown in Table 1, each bit contributes one small part of the reference voltage in the generation of the output voltage. The sum of the contribution of each individual bit then represents the analogue signal output. For example, in a DAC with 8 bits, the most significant bit contributes 0.5 of the reference voltage. The least significant bit contributes only 0.00390625 of the reference voltage. Each bit then contributes one specific gain value to the output voltage. The sum of all gains multiplied by the reference voltage then generates

A simple scheme for analogue signal generation using only the mathematical approach was created for this work. As shown in Fig. 4, the mathematical approach of R-2R ladder consists of an arrangement of *gain*, *sum* and *multiplier* blocks. Each of the *gain* blocks has a weighting factor that is added when its input is supplied with a logic 1 and multiplied by reference

In the example presented in Fig. 4, using *Vref = 1V* when *B3 = B0 = 1* and *B2 = B1 = 0*, the output value is *0.5625* which represents the addition of each input. The advantage of the approach is that the DAC using R-2R ladder can be easily simulated using the

*N =1* has the greatest effect on the output voltage, it is designated the MSB.

Table 1. Effect of individual bits on *Vout*

the network approaches *Vref*.

the final output voltage.

by *Vout = (Vref /2) + (Vref /8)* which reduces to *Vout = (5.Vref) /8*.

**2.2 Mathematical model that represents the R-2R ladder**

voltage, so generating the final output voltage.

MATLAB®/Simulink® or any other mathematical modelling environment. Note that this model is similar to R-2R ladder network presented in Fig. 3.

In this work, models of commercially available digital-to-analogue converters were developed. All converters considered had a resolution of 8 bits and used the R-2R ladder approach to implement the conversion circuitry. The IC designs used were the DAC08 (Analog Devices Inc. [a]), AD7524 (Analog Devices Inc. [b]) and AD5450 (Analog Devices Inc. [c]).

Today, the designer of digital circuits and systems typically uses hardware description languages (HDLs) such as VHDL or Verilog®-HDL in order model the circuit or system at a high level of design abstraction. Although these HDLs provide language constructs for behavioural simulation, their language subsets are far too restrictive for system level design. On the other hand, the MATLAB®/Simulink® modelling and simulation environment provides a powerful high level mathematical modelling environment for systems that can be widely used for high level algorithm development. It is a high-performance language for technical computing. It has a large collection of computational algorithms ranging from elementary functions, like summation and multiplication, through to complex arithmetic functions, such as fast Fourier transforms (FFTs).

A hardware description language, such as VHDL, can readily be used to describe a digital circuit or system. However, this language has a very limited capability when dealing with complex mathematical operations. The use of an interface, however, between MATLAB®/Simulink® and certain HDLs can fill the gap between an HDL and a high level mathematical modelling and simulation environment.

Fig. 4. Mathematical approach for modelling the R-2R ladder using Simulink®

In the next section, models developed in MATLAB®/Simulink® that implemented the digital-to-analogue converters selected for this work are presented.

#### **3. Data converter models implemented in MATLAB® /Simulink®**

A highly useful feature of Simulink® used in this work is the possibility of using it in the development of hybrid system models. In the DAC models developed, the use of *sum*, *multiplier*, *logic gates* and *shift register* components, along with components to generate input

A Methodology and Tool to

**3.1 DAC08 modelled in MATLAB®**

weight of each individual bit in the DAC08.

Table 2. Individual bit weighting in the DAC08 model

Fig. 6. DAC08 model implemented in Simulink®

Translate MATLAB®/Simulink® Models of Mixed-Signal Circuits to VHDL-AMS 389

The DAC08 is a multiplying D/A converter that uses the R-2R ladder in which an output current is the product of a digital number and an input reference current. Each digital input has a weight represented by the gain. These gains are added when their inputs have 1 logic and are multiplied by reference voltage, generating the output voltage. Table 2 shows the

> Bit position Weight S7 (MSB) 0.5 S6 0.25 S5 0.125 S4 0.0625 S3 0.03125 S2 0.015625 S1 0.00781250 S0 (LSB) 0.00390625

With this information it is possible to create the DAC08 converter model using the components available in the Simulink® libraries. Fig. 6 shows the DAC08 model implemented in Simulink®. It can be seen that the model is composed of the R-2R ladder network block, a

The multiplier block generates output voltage of the binary value in the digital input. The reference voltage will define the maximum output voltage. The commutation between the

As the R-2R ladder block is analogue in nature, there is a need to use the correct data type between the digital inputs and the analogue ladder network. In the R-2R ladder, the blocks

product block and a subtraction block to generate the output currents (*Io* and *Iob*).

outputs, labelled *Io* and *Iob*, is performed by the subtraction block.

**/Simulink®**

signals and undertake results analysis was required. All these basic components are available in libraries provided with Simulink®. These components are found in different libraries named *toolboxes*. A list of professional *toolboxes* currently available can be obtained from The MathWorks Inc. This list is by no means static; more toolboxes are being created every year. For example, there are toolboxes for communication systems, control systems, frequency-domain design, fuzzy logic and digital signal processing. With the components available to the user, it is possible to create models of the DACs selected in order to undertake a suitable design modelling and simulation study. Components are available to implement the analogue and digital circuit parts of the data converter. In addition, there are components to generate the input signals and to perform results analysis.

Using the available *and*, *or* and *not* logic gates, all required digital parts were created; shift *register*, *counters*, *latches* and all the logic necessary to implement the ADC7524 and ADC5450 *control* using *subsystems*.

As the model increases in size and complexity, the model visualization can be kept manageable by grouping blocks into suitably defined *subsystems*. For example, an R-2R ladder network can be generated as a subsystem and used in all the models. Fig. 5 shows the R-2R ladder created as a *subsystem*. This *subsystem* may then be used in all the models.

Fig. 5. R-2R ladder created as a *subsystem*

The *subsystem* named *R-2R* ladder is a block that has inputs labelled *S0* to *S7* (digital inputs) and an output *Vout* (analogue output). As shown in Fig. 5, the *subsystem* was created using the available *gain*, *sum*, *input* and *output* components, found in the Simulink® libraries.

#### **3.1 DAC08 modelled in MATLAB® /Simulink®**

388 Engineering Education and Research Using MATLAB

signals and undertake results analysis was required. All these basic components are available in libraries provided with Simulink®. These components are found in different libraries named *toolboxes*. A list of professional *toolboxes* currently available can be obtained from The MathWorks Inc. This list is by no means static; more toolboxes are being created every year. For example, there are toolboxes for communication systems, control systems, frequency-domain design, fuzzy logic and digital signal processing. With the components available to the user, it is possible to create models of the DACs selected in order to undertake a suitable design modelling and simulation study. Components are available to implement the analogue and digital circuit parts of the data converter. In addition, there are

Using the available *and*, *or* and *not* logic gates, all required digital parts were created; shift *register*, *counters*, *latches* and all the logic necessary to implement the ADC7524 and ADC5450

As the model increases in size and complexity, the model visualization can be kept manageable by grouping blocks into suitably defined *subsystems*. For example, an R-2R ladder network can be generated as a subsystem and used in all the models. Fig. 5 shows the R-2R ladder created as a *subsystem*. This *subsystem* may then be used in all the models.

The *subsystem* named *R-2R* ladder is a block that has inputs labelled *S0* to *S7* (digital inputs) and an output *Vout* (analogue output). As shown in Fig. 5, the *subsystem* was created using the available *gain*, *sum*, *input* and *output* components, found in the

components to generate the input signals and to perform results analysis.

*control* using *subsystems*.

Fig. 5. R-2R ladder created as a *subsystem*

Simulink® libraries.

The DAC08 is a multiplying D/A converter that uses the R-2R ladder in which an output current is the product of a digital number and an input reference current. Each digital input has a weight represented by the gain. These gains are added when their inputs have 1 logic and are multiplied by reference voltage, generating the output voltage. Table 2 shows the weight of each individual bit in the DAC08.

Table 2. Individual bit weighting in the DAC08 model

With this information it is possible to create the DAC08 converter model using the components available in the Simulink® libraries. Fig. 6 shows the DAC08 model implemented in Simulink®. It can be seen that the model is composed of the R-2R ladder network block, a product block and a subtraction block to generate the output currents (*Io* and *Iob*).

The multiplier block generates output voltage of the binary value in the digital input. The reference voltage will define the maximum output voltage. The commutation between the outputs, labelled *Io* and *Iob*, is performed by the subtraction block.

Fig. 6. DAC08 model implemented in Simulink®

As the R-2R ladder block is analogue in nature, there is a need to use the correct data type between the digital inputs and the analogue ladder network. In the R-2R ladder, the blocks

A Methodology and Tool to

B1 IS THE MORE SIGNIFICANT BIT

Fig. 8. Simulink® model used to test the DAC08

Translate MATLAB®/Simulink® Models of Mixed-Signal Circuits to VHDL-AMS 391

Fig. 9. Output analogue waveform resulting from the digital-to-analogue conversion

construction of latch with 12, 16 and 20 bits. Fig. 11 shows the *Data Latch* block.

*subsystem* named *convertor* generated to implement the data type conversion that was used outside of the R-2R ladder block. The 8-bit input for the *Data Latch* block was constructed by two *latch* blocks with 4 input bits and by a single control circuit. This modularity is for easy

labelled *Data type conversion* were used. These blocks implement the necessary conversion between different data types. Fig. 7 shows the used ladder R-2R network implemented as a subsystem. To test this model, an 8-bit sine waveform was generated as the converter signal input. To use the waveform as input of the model, the vectors *B1* (MSB) to *B8* (LSB) were generated from a MATLAB® M-file (with a file extension *.m*) and were used by Simulink® through the *From Workspace* block. The Simulink® model used to test the DAC08 is shown in Fig. 8 and the output signal, resulting from digital to analogue conversion, is shown in Fig. 9.

Fig. 7. R-2R ladder network implemented as *subsystem*

#### **3.2 AD7524 modelled in MATLAB® /Simulink®**

The AD7524 consists of an R-2R ladder and an input latch with 8 bits. The load cycle of the AD7524 is similar to the write cycle of a random access memory (RAM). In the *write mode*, when *CSb* (chip select – active low) and *WRb* (write – active low) are both LOW, the AD7524 is in the WRITE mode and the AD7524 analogue output responds to data activity at the *DB0-DB7* data bus inputs.

In all tests undertaken in the study, the *CSb* and *WRb* were set to a low level. In this mode, the AD7524 acts like a non-latched input D/A converter. Fig. 10 shows the AD7524 model implemented in Simulink®. We can see that the model is comprised of an R-2R ladder network block, a product block and a subtraction block to generate the outputs (*Out1* and *Out2*). This model also has *Data Latch* blocks with 8 bits of input. The R-2R ladder block is similar to the one used by the ADC08 converter presented in Fig. 5. In this model there is a

B1 IS THE MORE SIGNIFICANT BIT

390 Engineering Education and Research Using MATLAB

labelled *Data type conversion* were used. These blocks implement the necessary conversion between different data types. Fig. 7 shows the used ladder R-2R network implemented as a subsystem. To test this model, an 8-bit sine waveform was generated as the converter signal input. To use the waveform as input of the model, the vectors *B1* (MSB) to *B8* (LSB) were generated from a MATLAB® M-file (with a file extension *.m*) and were used by Simulink® through the *From Workspace* block. The Simulink® model used to test the DAC08 is shown in Fig. 8 and the output signal, resulting from digital to analogue conversion, is shown in Fig. 9.

Fig. 7. R-2R ladder network implemented as *subsystem*

**/Simulink®**

The AD7524 consists of an R-2R ladder and an input latch with 8 bits. The load cycle of the AD7524 is similar to the write cycle of a random access memory (RAM). In the *write mode*, when *CSb* (chip select – active low) and *WRb* (write – active low) are both LOW, the AD7524 is in the WRITE mode and the AD7524 analogue output responds to data activity at the

In all tests undertaken in the study, the *CSb* and *WRb* were set to a low level. In this mode, the AD7524 acts like a non-latched input D/A converter. Fig. 10 shows the AD7524 model implemented in Simulink®. We can see that the model is comprised of an R-2R ladder network block, a product block and a subtraction block to generate the outputs (*Out1* and *Out2*). This model also has *Data Latch* blocks with 8 bits of input. The R-2R ladder block is similar to the one used by the ADC08 converter presented in Fig. 5. In this model there is a

**3.2 AD7524 modelled in MATLAB®**

*DB0-DB7* data bus inputs.

Fig. 8. Simulink® model used to test the DAC08

Fig. 9. Output analogue waveform resulting from the digital-to-analogue conversion

*subsystem* named *convertor* generated to implement the data type conversion that was used outside of the R-2R ladder block. The 8-bit input for the *Data Latch* block was constructed by two *latch* blocks with 4 input bits and by a single control circuit. This modularity is for easy construction of latch with 12, 16 and 20 bits. Fig. 11 shows the *Data Latch* block.

A Methodology and Tool to

4 bits is shown in Fig. 13.

Fig. 13. Latch block with 4 bits

Fig. 14. Simulink® model to the ADC 7524 converter

Translate MATLAB®/Simulink® Models of Mixed-Signal Circuits to VHDL-AMS 393

The block for control of writing is a simple NOR gate. When the inputs labelled *CS\_b* and *W\_b* have low level logic applied, the output of the block is high level logic. In other words: *out = not (CS\_b + W\_b)*. This enables the latch to transfer the input to the output. Fig. 12 shows the control circuit to enable the transfer of data between the input and the output. The block labelled *latch 8* was created using two *latch* blocks with 4 bits. The *latch* block with

The same waveform was used for testing the Simulink® model as shown in Fig. 14.

Fig. 10. AD7524 model implemented in Simulink®

Fig. 11. 8-bit data latch block

Fig. 12. Block for control of writing

The block for control of writing is a simple NOR gate. When the inputs labelled *CS\_b* and *W\_b* have low level logic applied, the output of the block is high level logic. In other words: *out = not (CS\_b + W\_b)*. This enables the latch to transfer the input to the output. Fig. 12 shows the control circuit to enable the transfer of data between the input and the output. The block labelled *latch 8* was created using two *latch* blocks with 4 bits. The *latch* block with 4 bits is shown in Fig. 13.

Fig. 13. Latch block with 4 bits

392 Engineering Education and Research Using MATLAB

Fig. 10. AD7524 model implemented in Simulink®

Fig. 11. 8-bit data latch block

Fig. 12. Block for control of writing

The same waveform was used for testing the Simulink® model as shown in Fig. 14.

Fig. 14. Simulink® model to the ADC 7524 converter

A Methodology and Tool to

Fig. 16. AD5450 *data latch* block

Fig. 17. The *Control Load* Block

pulses.

Translate MATLAB®/Simulink® Models of Mixed-Signal Circuits to VHDL-AMS 395

*Pulses* block loads data into the latch to be converted after the count circuit receives 16 clock

#### **3.3 AD5450 modelled in MATLAB® /Simulink®**

The AD545n series (n = 0, 1, 2, 3) DAC has a 3-wire interface that is compatible with the majority of interface standards. Data is written to the device in 16-bit words. This word consists of two control bits and 14 data bits. The AD5450 uses 8 data bits and ignores the 6 LSBs. The controls bits *C1* and *C0* are used to load and update the code and change the active clock edge. Table 3 shows the control bits of the AD5450.

Table 3. AD5450 control bits