If I create two very simple signals, one an exponential decay and the other a gaussian, I can convolve this two signals very easily using scipy’s fft and ifft functions:
convolved = ifft(fft(decay) * fft(gaussian)).real
However, if I then try do recover the original exponential decay signal by doing:
decay_recovered = ifft(fft(convolved) / fft(gaussian)).real
The recovered decay is much more noisy then the original decay. Why? Is there any processing that can be applied to the fft of the concolved signal to prevent this?
[![enter image description here][1]][1]
Please see the code and the results below.
[![# Importing the necessary libraries.
import numpy as np
import pandas as pd
from scipy.fftpack import fft, ifft
import random
import matplotlib.pyplot as plt
def create_decay(x,I,t,q):
"""This function creates a power-law fluorescence decay from the x-axis
values (x), the initial fluorescence intensity value (I), the centre
of the lifetime distribution (t), and the heterogeneity parameter (q)."""
decay = ((2-q)/t)*(1-((1-q)*(x/t)))**(1/(1-q))
decay = decay * I
return decay
def convolve_IRF(irf,decay):
"""This function takes and Instrument Response Function and a fluorescence
decay and performs a convolution of the two."""
conv = ifft(fft(decay) * fft(irf)).real
return conv
def create_gaussian_curve(x):
"""This function creates a Gaussian curve using the x-axis values (x),
the center value (mu), and the sigma value (sig)."""
sig = round(random.uniform(0.07,0.10),2)
mu = round(random.uniform(1.10,1.40),2)
gauss = np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
return gauss
def add_Poisson_noise(decay):
"""This function adds Poison noise to a curve."""
noise_mask = np.random.poisson(decay)
decay = decay + noise_mask
return decay
def create_decay_real(x,irf):
"""This function creates a realistic fluorescence decay. The inputs are
the x-axis values (x), and an IRF (irf). The initial fluorescence value
(I), the centre value of the lifetime distribution, and the heterogeneity
parameter (q) are obtained randomly by generating values between 100 and
500 thousand (I), between 0.5 and 8 (t), and between 0.51 and 1.8 (q).
Gaussian noise is added to the decays. The decays are returned
normalised."""
I = random.randint(200,1000)
t = round(random.uniform(1.250,8.000),3)
q = round(random.uniform(0.950,1.990),3)
if q != 1.000:
decay_raw = create_decay(x,I,t,q)
decay_conv = convolve_IRF(irf,decay_raw)
decay_noisy = add_Poisson_noise(decay_conv)
decay_noisy = decay_noisy/np.max(decay_noisy)
return decay_raw, decay_noisy, t, q
# We create the x-axis values.
x = np.arange(0.04848485,49.6,0.09696969)
# Create IRF.
irf = create_gaussian_curve(x)
# We create the decay.
decay_raw, decay, t, q = create_decay_real(x,irf)
# We calculate what we need to come back to the pure decay.
fft_decay = fft(decay)
fft_irf = fft(irf)
# We obtain the pure decay.
pure_decay = ifft(fft_decay / fft_irf).real
# We normalise both decays.
decay_raw = decay_raw / np.max(decay_raw)
pure_decay = pure_decay / np.max(pure_decay)
# We plot this decay.
plt.plot(x,pure_decay,'k--',x,decay_raw,'r')
plt.savefig('test.png')