阶跃响应是否足以提供反卷积信息？

Question

我有一个测量系统，它以指数下降（蓝线，这将是测量数据）响应步长（绿线）。

我想使用反卷积从蓝线回到绿线。这个阶跃响应是否已经为反卷积提供了足够的信息，还是需要脉冲响应？

感谢您的帮助，

Answer 1

我遇到了同样的问题。我认为可以使用 Dirac delta 是derivative of Heaviside function 的事实来解决这个问题。您只需对阶跃响应进行数值导数，并将其用作反卷积的脉冲响应。

这是一个例子：

import numpy as np
from scipy.special import erfinv, erf
from scipy.signal import deconvolve, convolve, resample, decimate, resample_poly
from numpy.fft import fft, ifft, ifftshift

def deconvolve_fun(obs, signal):
    """Find convolution filter
    
    Finds convolution filter from observation and impulse response.
    Noise-free signal is assumed.
    """
    signal = np.hstack((signal, np.zeros(len(obs) - len(signal)))) 
    Fobs = np.fft.fft(obs)
    Fsignal = np.fft.fft(signal)
    filt = np.fft.ifft(Fobs/Fsignal)
    return filt
    
def wiener_deconvolution(signal, kernel, lambd=1e-3):
    """Applies Wiener deconvolution to find true observation from signal and filter
    
    The function can be also used to estimate filter from true signal and observation
    """
    # zero pad the kernel to same length
    kernel = np.hstack((kernel, np.zeros(len(signal) - len(kernel)))) 
    H = fft(kernel)
    deconvolved = np.real(ifft(fft(signal)*np.conj(H)/(H*np.conj(H) + lambd**2)))
    return deconvolved

def get_signal(time, offset_x, offset_y, reps=4, lambd=1e-3):
    """Model step response as error function
    """
    ramp_up = erf(time * multiplier)
    ramp_down = 1 - ramp_up
    if (reps % 1) == 0.5:
        signal =  np.hstack(( np.zeros(offset_x), 
                                ramp_up)) + offset_y
    else:
        signal =  np.hstack(( np.zeros(offset_x),
                                np.tile(np.hstack((ramp_up, ramp_down)), reps),
                                np.zeros(offset_x))) + offset_y 
      
    signal += np.random.randn(*signal.shape) * lambd  
    return signal
    
def make_filter(signal, offset_x):
    """Obtain filter from response to step function
    
    Takes derivative of Heaviside to get Dirac. Avoid zeros at both ends.
    """
    # impulse response. Step function is integration of dirac delta
    hvsd = signal[(offset_x):]
    dirac = np.gradient(hvsd)# + offset_y
    dirac = dirac[dirac > 0.0001]
    return dirac, hvsd

def get_step(time, offset_x, offset_y, reps=4):
    """"Creates true step response
    """
    ramp_up = np.heaviside(time, 0)
    ramp_down = 1 - ramp_up
    step =  np.hstack(( np.zeros(offset_x),
                        np.tile(np.hstack((ramp_up, ramp_down)), reps),
                        np.zeros(offset_x))) + offset_y          
    return step

# Worst case scenario from specs : signal Time t98%  < 60 s at 25 °C
multiplier = erfinv(0.98) / 60
offset_y = .01
offset_x = 300
reps = 1
time = np.arange(301)
lambd = 0
sampling_time = 3 #s

signal = get_step(time, offset_x, offset_y, reps=reps)
filter = get_signal(  time, offset_x, offset_y, reps=0.5, lambd=lambd)
filter, hvsd = make_filter(filter, offset_x)
observation = get_signal(time, offset_x, offset_y, reps=reps, lambd=lambd)
assert len(signal) == len(observation)
observation_est = convolve(signal, filter, mode="full")[:len(observation)]

signal_est = wiener_deconvolution(observation, filter, lambd)[:len(observation)]
filt_est = wiener_deconvolution(observation, signal, lambd)[:len(filter)]

这将允许您获得这两个数字：

Heaviside 和狄拉克

信号和滤波器估计

您还应该从检查我在代码中部分使用的其他 related posts 和 example of Wiener deconvolution 中受益。

如果这有帮助，请告诉我。