Une réponse d'étape est-elle suffisante pour la déconvolution?

Question

J'ai un système de mesure, qui répond à une étape (ligne verte) avec un déclin exponentiel (ligne bleue, qui serait les données mesurées). Je veux revenir de la ligne bleue à la ligne verte en utilisant la déconvolution. Cette étape-réponse est-elle déjà une information suffisante pour la déconvolution ou serait-il nécessaire d'avoir une réponse impulsionnelle?

Merci pour votre aide,

Answer 1

J'ai eu le même problème. Je pense qu'il peut être adressé en utilisant le fait que Dirac delta est un derivative of Heaviside function. Vous avez juste besoin de prendre la dérivée numérique de votre réponse d'étape et l'utiliser comme réponse impulsionnelle pour la déconvolution.

Voici un exemple:

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

Cela vous permettra d'obtenir ces deux chiffres:

Heaviside and Dirac

Signal and Filter Estimate

Vous devez également bénéficier de la vérification d'autres related posts et example of Wiener deconvolution que j'utilise en partie dans mon code. Faites-moi savoir si cela aide.