Is a step response enough information for deconvolution?

Question

I have a measurement system, which responds to a step (green line) with an exponential decline (blue line, which would be the measured data).

I want to go back from the blue line to the green line using deconvolution. Is this step-response already sufficient information for the deconvolution or would it be necessary to have the impulse response?

Thanks for your help,

Answer 1

I had the same problem. I think that it can be addressed using of the fact that Dirac delta is a derivative of Heaviside function. You just need to take numerical derivative of your step response and use it as impulse response for the deconvolution.

Here is an example:

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

This will allow you to obtain these two figures:

Heaviside and Dirac

Signal and Filter Estimate

You should also benefit from checking other related posts and the example of Wiener deconvolution that I partly use in my code.

Let me know if this helps.