Creating Convolution Matrices

By convention, bold upper case letters are used to denote all matrices. In code we use Unicode UTF-8 U+1D400 to U+1D419. It's easiest to copy paste the character from a UTF table.

Note

UTF-8 Captial Bold Letters: 𝐀, 𝐁, 𝐂, 𝐃, 𝐄, 𝐅, 𝐆, 𝐇, 𝐈, ...

Fredholm Integral

The response in channel i (corresponding to a transmission through the DMA at a single voltage) is given by the Fredholm convolution integral

$R_i = \int_{z_a}^{z_b} \sum_{k=1}^m \Omega(z,Z_{i,k}^s)T_c(D[z,1])T_l(D[z,1])\frac{dN}{d\ln D}\frac{d \ln D}{dz}dz \;\;\;\;\;\;\ i = 1,2...,n$

The integral is performed over the limits $z_a$ and $z_b$, which corresponds to the upper and lower mobility limit set by the voltage range used to operate the DMA. The function $\frac{dN}{d\ln D}\frac{d\ln D}{dz}dz$ evaluates to the number concentration of particles that lie in the interval $[z,z + dz]$. Note that $D[z,1]$ is used in the charge filter and loss function since the integral is applied over the transform of the selected centroid mobility $Z_{i,k}^s$. $Z$ is the mobility vector of the grid and the subscript $i$ denotes the response channel. The convolution integral can be discretized:

$R_i = \sum_{j=1}^n \left [ \sum_{k=1}^m \Omega(Z_j,Z_{i,k}^s)T_c(D[Z_j,1])T_l(D[Z_j,1])N(Z_j) \right]$

$N(Z_j)$ is the the number concentration of particles in the $j^{th}$ bin, $i = 1...n$ are indices the observed instrument channel, $j = 1...n$ are indices of the physical size bins, and $k = 1...m$ are indices of charges carried by the particle. Here it is assumed that $\Omega(Z_j,Z_{i,k}^s)$ can be approximated being constant over the bin $[Ze_{j},Ze_{j+1}]$, which is only true if a sufficiently large number of size bins is used.

The discretized version represents a set of $n$ equations that can be written in matrix form

$R = \mathbf{A}N$

where $R$ is the response vector, $\mathbf{A}$ is the convolution matrix, and $N$ is the discretized true number concentration.

Computing the Convolution Matrix

The terms comprising the convolution matrix

$\Omega(Z_1,Z_{1,k}^s)T_c(D[Z_1,1])T_l(D[Z_1,1])$

can be written as

T(zˢ,k,Λ,δ) = δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp)

The generic solution to constructing the convolution matrix is to define a forward transmission model and collect the terms. The resulting matrix is $n \times n$ square matrix, where $n$ equals to the number of bins.

bins,z₁,z₂ = 60, vtoz(Λ,10000), vtoz(Λ,10)
δ = setupDMA(Λ, z₁, z₂, bins)
T(zˢ,k,Λ,δ) = δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp)
𝐀 = (hcat(map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)...))'

60×60 adjoint(::Matrix{Float64}) with eltype Float64:
 0.131588     0.0540361    1.19965e-16  …  4.19054e-19  0.0
 0.0605105    0.137591     0.0564765       0.0          3.64643e-51
 2.58731e-5   0.0633166    0.143694        2.96024e-33  0.0
 2.19162e-17  3.75145e-5   0.0661762       9.05185e-15  1.81298e-34
 8.031e-18    4.99485e-17  5.3083e-5       0.0          1.9204e-34
 0.0320654    5.75813e-17  9.03059e-17  …  3.57155e-49  0.0
 0.0843208    0.0325587    2.53042e-17     0.0          0.0
 0.0398241    0.0855716    0.0328725       0.0          4.67259e-73
 0.00115426   0.0404367    0.0863451       0.0          0.0
 0.0294711    0.00112585   0.0408265       0.0          2.62136e-73
 ⋮                                      ⋱               
 0.0          0.0          0.0             1.75125e-17  0.0
 0.0          0.0          0.0             0.0          3.71538e-37
 0.0          0.0          0.0             1.96499e-17  0.0
 0.0          0.0          0.0             7.95257e-12  2.70424e-22
 0.0          5.55016e-17  0.0          …  3.61098e-6   4.88025e-11
 0.0          8.74369e-19  0.0             0.00114259   6.02266e-6
 0.0          0.0          0.0             0.00980112   0.00115801
 0.0          0.0          0.0             0.0161758    0.00869253
 0.0          0.0          0.0             0.0107893    0.0139506

See Notebook S2 in the Notebooks section for a step-by-step derivation. For a narrated description check out Session 2 of the Tutorial.

Precomputed Matrices

The following matrices are precomputed for each DMA and stored in δ. The structure is mutable and the matrices can be altered if needed.

Matrix 𝐀

The matrix 𝐀 describes transmission through the DMA with neutralizer and transmission loss function. The matrix 𝐀 is useful for modeling the measured response function of a size distribution passing through a stepping or scanning DMA.

δ = setupDMA(Λ, z₁, z₂, bins)                   
T(zˢ,k,Λ,δ) = δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp).*δ.Tl(Λ,δ.Dp)
𝐀 = (hcat(map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)...))'

The type δ is mutable. You can therefore override the precomputed matrix by creating your own and manually adding it to the DMA grid.

T(zˢ,k,Λ,δ) = δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tc(k,δ.Dp)
my𝐀 = (hcat(map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)...))'
δ.𝐀 .= my𝐀

Matrix 𝐎

The matrix 𝐎 describes transmission through the DMA without neutralizer and with a transmission loss function. The matrix 𝐎 is useful for modeling the measured response function of a known mobility distribution passing through the second DMA in tandem DMA setups.

T(zˢ,k,Λ,δ) = δ.Ω(Λ,δ.Z,zˢ/k).*δ.Tl(Λ,δ.Dp)
𝐎 = (hcat(map(zˢ->Σ(k->T(zˢ,k,Λ,δ),Λ.m),δ.Z)...))'

Matrix 𝐒

Talukdar and Swihart (2003) introduced the matrix 𝐒: sum the rows of 𝐀 and place the results on the diagonal of 𝐒. The matrix 𝐒 is used to compute an initial guess to constrain the Tikhonov inverse.

𝐒 = zeros(bins, bins)
for i = 1:bins
	@inbounds 𝐒[i, i] = sum(𝐀[i, :])
end

Matrix 𝐈

The identity matrix is used as weights matrix when computing the Tikhoniv inverse.

𝐈 = Matrix{Float64}(I, bins, bins)