# Implementation of Mellin transform of exponential decay

I'm trying to understand this paper: 10.1016/j.jmr.2010.05.015. It is about using a Mellin transform of curves that contain multiple exponential decays of varying contributions (CPMG data from Nuclear Magnetic Resonance experiments). I'm going to copy here the important equations here, and then show the code I wrote in Matlab trying to implement them, and how my graphs are different from theirs, and then ask a few questions. I must preface that I'm totally out of my field, so bear with me.

To summarize what I want, I'm trying to plot $$G(\omega)$$ by $$\omega$$ (Fig2B). Here, $$G(\omega)=\ln\langle T_2^\omega \rangle$$ (Eq. 9), and Eq. 5 is

$$\displaystyle \langle T_2^\omega \rangle = \frac{(-1)^n}{\Gamma(\mu)\phi} \int_0^\infty t^{(\mu-1)} \left[ \frac{d^n M(t)}{dt^n} \right]dt$$

where $$\phi$$ is porosity, $$t$$ is time, $$M(t)$$ is the data (I'm simulating the simplest case, a monoexponential decay of type $$\exp(-t/0.1)$$), $$\omega = \mu - n$$ with $$n=0$$ (no derivative) if $$\omega > 0$$ or $$n=\left[-\omega \right] + 1$$. $$\Gamma$$ is the typical gamma function (factorial for non integers) and $$\phi$$ is just a constant (porosity), which I'm setting to 1 (by their own words, it's the unmeasured datum at t=0, so for this simple decay, it's 1). In this paper, they used $$\left[-\omega\right]$$ to mean the integer part of $$\omega$$.

Later on, they give some equations to calculate these values using discrete sums, which is what I'm after. My final objective is to perform a measurement, get $$M(t)$$ and use this Mellin transform they proposed to get some information on the distribution of exponential decays.

They divided the calculation in there's three cases:

1. If $$\omega = 0$$ then, by definition, $$\langle T_2^{\omega=0} \rangle = 1$$
2. If $$\omega > 0$$, which leads to their Eq. 18

$$\langle T_2^\omega \rangle = k + \frac{1}{\Gamma(\omega+1) \phi} \sum_{i=1}^N \Delta_i M(it_E)$$

where $$k= \frac{\tau_{min}}{\Gamma(\omega + 1)}$$, $$\tau_{min} = t_E^\omega$$ ($$\tau_E$$ is the spacing between the measurements) and the deltas are $$\Delta_1 = 0.5 t_E^\omega [2^\omega - 1^\omega]$$, $$\Delta_i = 0.5 t_E^\omega [(i+1)^\omega - (i-1)^\omega]$$, $$\Delta_N = 0.5 t_E^\omega [N^\omega - (N-1)^\omega]$$. I'm assuming here that $$[\text{number}]$$ still means the integral part of a number.

1. If $$-1 < \omega < 0$$, leads to their Eq. 22

$$\langle T_2^\omega \rangle = k + \frac{1}{\Gamma(\omega + 1)\phi} \left[ \left(\frac{a_1 \omega}{\omega + 1} \right) \tau_{min}^{(\omega+1)/\omega} + \sum_{i=1}^N \Delta_i M(it_E) \right]$$

where the Deltas and k were given above and $$a_1$$ is the derivative of $$M(t)$$ at $$t=0$$ (for my toy exponential, this is -10).

I tried implementing this, but my resulting graph is not at all the straight line they showed. It has a rounded part close to $$\omega=1$$, then it goes straight down, then up with some jagged edges (resulting from using only the integral parts?), then down again. They state in a later paper that this function should be continuous. I've checked the equations dozens of times I can't find what I'm doing wrong. Here's where I think I might have a problem.

1. I used Matlab's fix to calculate $$[\text{number}]$$. Is this the right function to calculate the "integer part" of a number? Are some the instances of $$[\text{number}]$$ the application of this, or are some of them just following the sequence $$\{[()]\}$$?
2. If I set $$\phi$$ to a large value (10000), the curve becomes straight in the range I'm looking at, but this doesn't make any sense, either by their own description, or physically (porosities should be maxed at 100%, of course). If I set $$\phi$$ to between 0 and 1, I still get the jagged edges I mentioned before.
3. I'm overthinking and those jagged edges are unimportant, because I only need the first and second derivatives of $$G(\omega)$$ at $$\omega=0$$ (their own conclusion), so $$G(\omega)$$ only has to be well behaved in this region. Still doesn't explain why their graphs do not contain these edges.
4. I am missing something obvious and can't get find it.

Here's my matlab code. I tried to comment it so it's easy to follow.

% Setting up stuff
dead_time = 0.01; %seconds
num_points = 1000;
t = linspace(dead_time, 0.8, num_points); % Possible decay range, Fig 3A
t_echo = t(2) - t(1); % Spacing
t2 = 0.1; % Fig 2A
exp_decay = exp(-t/t2); % single decay, Fig 2A, trying to get Fig2B. This is M(t)
phi = 1; % The equation at t=0

% Calculating the estimated moments, Sec 3.3
num_omegas = 500;
omegas = linspace(-1, 4, num_omegas); % From Fig2B
results = zeros(num_omegas, 1); % <T_2^omega>

for i = 1:num_omegas
omega = omegas(i);
if omega == 0 % Case 1
results(i) = 1;
elseif omega > 0 % Case 2, Eq. 18
tau_min = t_echo ^ omega; % Eq. 19b
k = tau_min / gamma(omega + 1); % Eq. 19a
% Calculate the sum using deltas from 19c-19e)
sum = 0;
for j = 1:num_points
% [omega] according to them is "the integral part of number omega". Interpreted this as fix(omega).
if j == 1; DeltaJ = 0.5 * t_echo ^ omega * fix(2 ^ omega - 1 ^ omega); % Eq. 19c.
elseif j == num_points; DeltaJ = 0.5 * t_echo ^ omega * fix(num_points ^ omega - (num_points - 1) ^ omega);
else; DeltaJ = 0.5 * t_echo ^ omega * fix((j + 1) ^ omega - (j - 1) ^ omega);
end
sum = sum + DeltaJ * exp_decay(j); % i t_E is the i-th echo, i.e., the i-th measurement
end
results(i) = k + 1 / (gamma(omega + 1) * phi) * sum;
elseif (omega > -1) && (omega < 0)  % Case 3, Eq. 22
tau_min = t_echo ^ omega; % Eq. 19b
k = tau_min / gamma(omega + 1); % Eq. 19a
a1 = -10; % From Eq. 22. d/dt e(-t/0.1) = -10e^(-10t) |t=0 = -10.
% Calculate the sum using deltas from Eq. 19c-19e
sum = 0;
for j = 1:num_points
if j == 1; DeltaJ = 0.5 * t_echo ^ omega * fix(2 ^ omega - 1 ^ omega);
elseif j == num_points; DeltaJ = 0.5 * t_echo ^ omega * fix(num_points ^ omega - (num_points - 1) ^ omega);
else; DeltaJ = 0.5 * t_echo ^ omega * fix((j+1) ^ omega - (j - 1) ^ omega);
end
sum = sum + DeltaJ * exp_decay(j);
end
% Eq. 22
results(i) = k + 1 / (gamma(omega+1) * phi) * fix( (a1 * omega)/(omega+1) * tau_min ^ ((omega+1)/omega) + sum );

end
end

% Eq. 9
G_omega = log(results);
plot(omegas, G_omega); xlabel('\omega'); ylabel('G(\omega)')


Here's the resulting graph. As you can see, it's very different from the one in the paper. I don't know if the data I'm feeding is wrong, or if I've implemented the algorithm wrong. I would appreciate a fresh set of eyes on this, very much.