5
$\begingroup$

Forewords

This question is about methodology references and numerical application. I am posting on Signal Processing because I think this question belong to this place. I am new to the stack, feel free to edit if you feel the Original Post should be improved.

Experiment Setup

Basic experiment setup: inject a dye into a reactor to analyze its dynamic. There are two detector cells (sensitive and selective to dye substance) at in entrance and another at the exit of the reactor. No reaction occurs, the dye is a passive tracer in this experiment.

Two signals are recorded:

  • $x(t)$, a finite impulse (eg.: we realized a gaussian-like impulse, at the entrance);
  • and $y(t)$, the response of a system to this impulse (at the exit);

Goal is to compute $g(t)$: the system Retention Time Distribution (RTD).

Signal Processing

For the MCVE sake, I created a synthetic dataset to perform computation on well controlled scenario without noise.

enter image description here

Where $y(t)$ is a convolution such as:

$$y(t) = g(t)*x(t) \xrightarrow{\mathcal{F}} Y(\nu) = G(\nu)\cdot X(\nu)$$

Now, it might resume to get the inverse transform of the transfer function:

$$ G(\nu) = \frac{Y(\nu)}{X(\nu)} $$

But numerically, even for this smooth example, it does not work well. Mainly because $X(\nu)$ has zeros and small numbers that make the operation unstable.

Alternative

There is a processing (adapted from this example) that seems to solve my problem:

$$ G'(\nu) = \frac{X^*(\nu)\cdot Y(\nu)}{X^*(\nu)\cdot X(\nu) + \lambda^2} $$

Where $\lambda$ is related to Signal-to-Noise Ratio:

$$ \lambda^2 = \frac{1}{\mathrm{SNR}} $$

I feel comfortable with this solution, because the operation tends to the former when $\lambda$ goes to zero:

$$\lim\limits_{\lambda \rightarrow 0} G'(\nu) = G(\nu)$$

Then I can estimate RTD by taking the inverse transform of $G'(\nu)$:

$$\hat{g}(t) = \mathcal{F}^{-1}[G'(\nu)]$$

Numerical Example

For the MCVE dataset, it renders as follow:

enter image description here

Derived RTD (after inverse transformation) quite fits the original function:

enter image description here

And error on this method is acceptable for the MCVE:

enter image description here

Even if there are artifacts at the very beginning and end of the recovered signal.

Observations

Setting $\lambda$ to higher values change the error behavior, but the RTD looks fine:

enter image description here

Setting $\lambda$ to smaller value may affect the shape of the RTD:

enter image description here

enter image description here

Questions

Reading on the internet I had difficulties to find references where the transfer function is the unknown. Most of the articles speaks about deconvolution and estimation of $\hat{x}(t)$ knowing $y(t)$ and $g(t)$.

I also found the Wikipedia page of Wiener Deconvolution which shows a formulae that is very similar to $G'(\nu)$ in my OP. But I could not relate them together. I miss something to get the full understanding on why this methodology is correct, if it is, when performed as I did.

My questions are:

  • How can I formally justify the expression of $G'(\nu)$? Can it be shown there is a relation with Wiener Filter?
  • Is there a good reference to apprehend this class of problem?
  • How is called the oscillating effect at the ends of the signal?

Your help is greatly appreciated.

Code

Below the Python code of the MCVE:

import numpy as np
from scipy import stats, signal, fftpack

# Setup:
fs = 10   # Sampling Frequency [Hz]
dt = 600  # Experiment runtime [s]
l = 1e-4  # 1/SNR in term of amplitude [-]
t = np.linspace(0, dt, dt*fs+1) # Time [s]

# Signals [AU]:
xt = stats.norm(loc=20, scale=5).pdf(t)
gt = stats.lognorm(0.45, loc=-20, scale=100).pdf(t)
yt = signal.convolve(xt, gt, mode='full')
yt = yt[:t.size]/np.sum(yt[:t.size]/fs) # Normalization

# Fourier Transform:
Xs = fftpack.fft(xt)
Ys = fftpack.fft(yt)
fq = fs*fftpack.fftfreq(t.size)
Gs = (np.conj(Xs)*Ys)/(Xs*np.conj(Xs) + l**2) # Transfer function

# Inverse Transform:
gti = fftpack.ifft(Gs)
np.allclose(np.imag(gti), 0) # True
gti = np.real(gti)
gti = gti/(gti.sum()/fs) # Normalization

# Error:
ea = gti-gt
er = ea/gt
$\endgroup$
  • 3
    $\begingroup$ Well-detailed question, filled with personal experience, code and results. I hope you will get useful answers $\endgroup$ – Laurent Duval Nov 7 at 21:26
1
$\begingroup$

How can I justify the expression of G′(ν)?

That's easy enough. Denominator is the sum of the signal energy $|X(\omega)|^2$ and the noise energy $\lambda ^2$. If the signal energy is significantly larger, then the whole expression simplifies to $G(\omega)$. If the noise is larger, we can't do a anything useful with the information and the $1/ \lambda ^2$ term will drive the whole expression towards zero and some continuous way.

Is there a good reference to apprehend this class of problem?

Deconvolution does work fine here. Convolution is commutative, so it doesn't really matter if you are trying to calculate the impulse response or the input signal (at least in your case).

There are a few other methods that can also work well which are typically based on a least square error optimization in either the time or the frequency domain.

Any of these methods can be greatly improved by using any a-priori knowledge of the system. For example, if you that it's a relatively well behaved physical system with a some known degrees of freedom, you can simply do a least square error fit trying to determine a small number of poles and zeros.

Each application has it's own specific requirement and the best method depends a lot on what you have and what exactly you need.

How is called the oscillating effect at the ends of the signal?

Not sure. I think your problem here may be that your noise estimate is not a function of frequency and that somehow a single or a few discrete frequencies get overal amplified in the error signal.

$\endgroup$
  • $\begingroup$ Thank you for answering, you brought insights. Yes it is clear that $G'(\nu)$ tends to $G(\nu)$ when noise vanishes. Justification in term of energy was enough to me to perform the computation. Then when I saw it worked, I tried to formally justify this formulae, sorry if I was not clear about that specific point. I am intrigued by the similarity with Wiener Filter expression but I could not relate it. Any chance you know if there is a relation between those expressions. $\endgroup$ – jlandercy Nov 9 at 8:13
  • $\begingroup$ I created the MCVE without adding noise to make it as simple as possible. Regions where oscillations occur is when signals (which are distributions) are close to machine precision. It might be a float precision error combined with your suggestion. I will investigate it further. $\endgroup$ – jlandercy Nov 9 at 8:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.