समझना scipy deconvolve

Question

मैं scipy.signal.deconvolve को समझने की कोशिश कर रहा हूं।

देखने के गणितीय दृष्टि से एक घुमाव के बस फूरियर अंतरिक्ष में गुणन तो मैं उम्मीद होती है दो कार्यों f और g के लिए कि:
Deconvolve(Convolve(f,g) , g) == f

numpy में/SciPy यह या तो मामला नहीं या है मुझे एक महत्वपूर्ण बात याद आ रही है। हालांकि एसओ पर deconvolve से संबंधित कुछ प्रश्न हैं (जैसे here और here) वे इस बिंदु को संबोधित नहीं करते हैं, अन्य अस्पष्ट रहते हैं (this) या अनुत्तरित (here)। सिग्नलप्रोसेसिंग एसई (this और this) पर दो प्रश्न भी हैं जिनके जवाब समझने में सहायक नहीं हैं कि कैसे scipy deconvolve फ़ंक्शन काम करता है।

सवाल होगा:

• आप कैसे एक जटिल संकेत से मूल संकेत f फिर से संगठित है, संभालने क्या आप जानते हैं convolving समारोह छ।?
• या दूसरे शब्दों में: यह छद्म कोड Deconvolve(Convolve(f,g) , g) == f कैसे numpy/scipy में अनुवाद करता है?

संपादित: ध्यान दें कि यह प्रश्न संख्यात्मक अशुद्धियों को रोकने में लक्षित नहीं है (हालांकि यह भी एक open question है), लेकिन scipy में एक साथ कैसे समझ convolve/deconvolve काम पर।

निम्न कोड हेवीसाइड फ़ंक्शन और एक गाऊशियन फ़िल्टर के साथ ऐसा करने का प्रयास करता है। छवि में देखा जा सकता है, संकल्प के deconvolution का परिणाम सभी मूल हेवीसाइड समारोह पर नहीं है। अगर कोई इस मुद्दे में कुछ प्रकाश डाल सकता है तो मुझे खुशी होगी।

import numpy as np 
import scipy.signal 
import matplotlib.pyplot as plt 

# Define heaviside function 
H = lambda x: 0.5 * (np.sign(x) + 1.) 
#define gaussian 
gauss = lambda x, sig: np.exp(-(x/float(sig))**2) 

X = np.linspace(-5, 30, num=3501) 
X2 = np.linspace(-5,5, num=1001) 

# convolute a heaviside with a gaussian 
H_c = np.convolve(H(X), gauss(X2, 1), mode="same" ) 
# deconvolute a the result 
H_dc, er = scipy.signal.deconvolve(H_c, gauss(X2, 1)) 


#### Plot #### 
fig , ax = plt.subplots(nrows=4, figsize=(6,7)) 
ax[0].plot(H(X),   color="#907700", label="Heaviside", lw=3) 
ax[1].plot(gauss(X2, 1), color="#907700", label="Gauss filter", lw=3) 
ax[2].plot(H_c/H_c.max(), color="#325cab", label="convoluted" , lw=3) 
ax[3].plot(H_dc,   color="#ab4232", label="deconvoluted", lw=3) 
for i in range(len(ax)): 
    ax[i].set_xlim([0, len(X)]) 
    ax[i].set_ylim([-0.07, 1.2]) 
    ax[i].legend(loc=4) 
plt.show() 

संपादित: नोट/

yc=conv(y,c,'full')./sum(c); 
ydc=deconv(yc,c).*sum(c); 

का उपयोग कर इसे इस सवाल की भावना में एक आयताकार संकेत deconvolve है कि वहाँ एक matlab example, कैसे convolve को दिखा अगर कोई इस उदाहरण को पाइथन में अनुवाद करने में सक्षम था तो भी मदद करेगा।

Answer 1

कुछ परीक्षण और त्रुटि के बाद मुझे पता चला कि scipy.signal.deconvolve() के परिणामों का व्याख्या कैसे करें और मैं अपने निष्कर्षों को उत्तर के रूप में पोस्ट करता हूं।

:

के एक काम उदाहरण कोड

import numpy as np 
import scipy.signal 
import matplotlib.pyplot as plt 

# let the signal be box-like 
signal = np.repeat([0., 1., 0.], 100) 
# and use a gaussian filter 
# the filter should be shorter than the signal 
# the filter should be such that it's much bigger then zero everywhere 
gauss = np.exp(-((np.linspace(0,50)-25.)/float(12))**2) 
print gauss.min() # = 0.013 >> 0 

# calculate the convolution (np.convolve and scipy.signal.convolve identical) 
# the keywordargument mode="same" ensures that the convolution spans the same 
# shape as the input array. 
#filtered = scipy.signal.convolve(signal, gauss, mode='same') 
filtered = np.convolve(signal, gauss, mode='same') 

deconv, _ = scipy.signal.deconvolve(filtered, gauss) 
#the deconvolution has n = len(signal) - len(gauss) + 1 points 
n = len(signal)-len(gauss)+1 
# so we need to expand it by 
s = (len(signal)-n)/2 
#on both sides. 
deconv_res = np.zeros(len(signal)) 
deconv_res[s:len(signal)-s-1] = deconv 
deconv = deconv_res 
# now deconv contains the deconvolution 
# expanded to the original shape (filled with zeros) 


#### Plot #### 
fig , ax = plt.subplots(nrows=4, figsize=(6,7)) 

ax[0].plot(signal,   color="#907700", label="original",  lw=3) 
ax[1].plot(gauss,   color="#68934e", label="gauss filter", lw=3) 
# we need to divide by the sum of the filter window to get the convolution normalized to 1 
ax[2].plot(filtered/np.sum(gauss), color="#325cab", label="convoluted" , lw=3) 
ax[3].plot(deconv,   color="#ab4232", label="deconvoluted", lw=3) 

for i in range(len(ax)): 
    ax[i].set_xlim([0, len(signal)]) 
    ax[i].set_ylim([-0.07, 1.2]) 
    ax[i].legend(loc=1, fontsize=11) 
    if i != len(ax)-1 : 
     ax[i].set_xticklabels([]) 

plt.savefig(__file__ + ".png") 
plt.show()  

इस कोड को निम्न छवि पैदा करता है, दिखा हम क्या चाहते हैं (Deconvolve(Convolve(signal,gauss) , gauss) == signal)

कुछ महत्वपूर्ण निष्कर्ष हैं साथ शुरू करते हैं

• फ़िल्टर शोर होना चाहिए संकेत
• फिल्टर से Ter हर जगह शून्य की तुलना में काफी बड़ा होना चाहिए (यहाँ> 0.013 काफी अच्छा है)
• कीवर्ड तर्क mode = 'same' घुमाव के लिए सुनिश्चित करता है कि यह संकेत के रूप में एक ही सरणी आकार पर रहता है का उपयोग करना।
• deconvolution n = len(signal) - len(gauss) + 1 अंक है। तो इसे उसी मूल सरणी आकार पर भी रहने के लिए हमें इसे दोनों तरफ s = (len(signal)-n)/2 तक विस्तारित करने की आवश्यकता है।

बेशक, इस प्रश्न के लिए और निष्कर्ष, टिप्पणियां और सुझाव अभी भी स्वागत है।

Answer 2

टिप्पणियों में लिखे गए अनुसार, मैं आपके द्वारा मूल रूप से पोस्ट किए गए उदाहरण के साथ मदद नहीं कर सकता। जैसा कि @ स्टेलियोस ने इंगित किया है, संख्यात्मक मुद्दों के कारण deconvolution काम नहीं कर सकता है।

मैं, फिर भी, उदाहरण आप अपने संपादित में तैनात पुन: पेश कर सकते हैं:

कोड जो matlab स्रोत कोड से एक सीधा अनुवाद है कि:

import numpy as np 
import scipy.signal 
import matplotlib.pyplot as plt 

x = np.arange(0., 20.01, 0.01) 
y = np.zeros(len(x)) 
y[900:1100] = 1. 
y += 0.01 * np.random.randn(len(y)) 
c = np.exp(-(np.arange(len(y)))/30.) 

yc = scipy.signal.convolve(y, c, mode='full')/c.sum() 
ydc, remainder = scipy.signal.deconvolve(yc, c) 
ydc *= c.sum() 

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(4, 4)) 
ax[0][0].plot(x, y, label="original y", lw=3) 
ax[0][1].plot(x, c, label="c", lw=3) 
ax[1][0].plot(x[0:2000], yc[0:2000], label="yc", lw=3) 
ax[1][1].plot(x, ydc, label="recovered y", lw=3) 

plt.show()