27

私は理解しようとしていますscipy.signal.deconvolve。

数学的な観点からは、畳み込みはフーリエ空間での乗算にすぎないため、2 つの関数fとについては次のようになりgます。
Deconvolve(Convolve(f,g) , g) == f

numpy/scipy では、これは当てはまらないか、重要な点が抜けています。SO のデコンボリューションに関連するいくつかの質問がありますが (ここやここのように)、この点に対処していません。他の質問は不明 (これ) または未回答 (ここ) のままです。また、SignalProcessing SE に関する 2 つの質問 ( thisおよびthis ) に対する回答は、scipy の deconvolve 関数がどのように機能するかを理解するのに役立ちません。

問題は次のとおりです。

  • 畳み込み関数 g を知っていると仮定して、畳み込み信号から元の信号をどのように再構築fしますか?
  • 言い換えれば、この疑似コードはどのようDeconvolve(Convolve(f,g) , g) == fに numpy / scipy に変換されるのでしょうか?

編集：この質問は、数値の不正確さを防ぐことを目的としているのではなく（これも未解決の質問です）、scipyで畳み込み/デコンボリューションがどのように連携するかを理解することを目的としていることに注意してください。

次のコードは、Heaviside 関数とガウス フィルターを使用してこれを実行しようとしています。画像でわかるように、畳み込みのデコンボリューションの結果は元のヘビサイド関数ではありません。誰かがこの問題に光を当てることができれば幸いです。

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

ここに画像の説明を入力

編集:を使用して矩形信号を畳み込む/デコンボリューションする方法を示す、 matlab の例があることに注意してください

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

この質問の精神で、誰かがこの例を python に翻訳できればそれも役に立ちます。

4

2 に答える 2

19

いくつかの試行錯誤の後、結果を解釈する方法を見つけ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)

ここに画像の説明を入力

いくつかの重要な調査結果は次のとおりです。

  • フィルタは信号よりも短くする必要があります
  • フィルターはどこでもゼロよりもはるかに大きくする必要があります (ここでは > 0.013 で十分です)
  • 畳み込みにキーワード引数mode = 'same'を使用すると、信号と同じ配列形状に存在することが保証されます。
  • デコンボリューションにはn = len(signal) - len(gauss) + 1ポイントがあります。したがって、同じ元の配列形状にも存在させるにはs = (len(signal)-n)/2、両側で拡張する必要があります。

もちろん、この質問に対するさらなる調査結果、コメント、提案は大歓迎です。

于 2016-11-17T13:32:48.170 に答える
7

コメントに書かれているように、最初に投稿した例についてはお手伝いできません。@Stelios が指摘したように、数値の問題により、デコンボリューションがうまくいかない可能性があります。

ただし、編集で投稿した例を再現できます。

ここに画像の説明を入力

これは、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()
于 2016-11-17T13:36:02.283 に答える