为什么 scipy 的快速傅立叶变换操作会给数据添加如此多的噪声？

Question

如果我创建两个非常简单的信号，一个是指数衰减，另一个是高斯信号，我可以使用 scipy 的

fft

和

ifft

函数非常轻松地对这两个信号进行卷积：

convolved = ifft(fft(decay) * fft(gaussian)).real

但是，如果我尝试通过执行以下操作来恢复原始指数衰减信号：

decay_recovered = ifft(fft(convolved) / fft(gaussian)).real

恢复的衰减比原始衰减的噪声要大得多。为什么？是否有任何处理可以应用于卷积信号的 fft 来防止这种情况？

请参阅下面的代码和结果。

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

Answer 1

我最终意识到噪声不是来自傅里叶变换，而是来自我添加到信号中的泊松噪声。如果噪声被去除，重建的信号与原始信号相同。