signal-processingdeconvolution

Is a step response enough information for deconvolution?


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

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,


Solution

  • 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 Heaviside and Dirac

    Signal and Filter Estimate 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.