Why do fast fourier transform operations from scipy add so much noise to the data?

  Kiến thức lập trình

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')

LEAVE A COMMENT