2017-05-23 4 views
0

나는 계량 시스템을 가지고 있는데, 지수계 감소 (파란색 선, 측정 된 데이터)가있는 단계 (녹색 선)에 응답합니다. enter image description heredeconvolution에 대한 단계 응답에 충분한 정보가 있습니까?

나는 deconvolution을 사용하여 파란색 선에서 녹색 선으로 돌아가고 싶습니다. 이 단계 응답은 이미 deconvolution에 충분한 정보입니까 아니면 임펄스 응답이 필요합니까?

도움을 주셔서 감사합니다.

답변

0

동일한 문제가있었습니다. Dirac delta가 derivative of Heaviside function이라는 사실을 사용하여 해결할 수 있다고 생각합니다. 당신은 단지 스텝 반응의 수치 적 미분을 취하여 이것을 데콘 볼 루션에 대한 임펄스 응답으로 사용할 필요가 있습니다. 다른 related posts 및 확인에서

Heaviside and Dirac

Signal and Filter Estimate

또한 수혜가 예상 :

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

이것은이 두 수치를 얻을 수 있습니다 : 여기

은 예입니다 내 코드에서 부분적으로 사용하는 example of Wiener deconvolution.

도움이 될지 알려주세요.

관련 문제