2017-05-23 3 views
0

私は、指数関数的な低下(測定データとなる青い線)を伴うステップ(緑色の線)に応答する測定システムを持っています。 enter image description hereデコンボリューションのステップ応答で十分な情報ですか?

デコンボリューションを使用して青い線から緑の線に戻っていきたいです。このステップ応答は既にデコンボリューションのための十分な情報か、インパルス応答を持つ必要がありますか?

ご協力ありがとうございます。

答えて

0

私は同じ問題がありました。 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)] 

これは、あなたがこれらの二つの数字得ることができるようになります:あなたはまた、他のrelated postsとチェックの恩恵を受けられるはずです

Heaviside and Dirac

Signal and Filter Estimate

example of Wiener deconvolution私は自分のコードで部分的に使用しています。

これが役立つかどうか教えてください。

関連する問題