Как использовать БПФ для 1D деконволюции? - PullRequest
Новая
       26

Как использовать БПФ для 1D деконволюции?

0 голосов
/ 25 октября 2019

Задача
Я пытаюсь деконволюировать два измеренных данных A и B , используя теорему свертки. Я знаю, что для свертки вы должны обнулить свои данные, чтобы предотвратить циклическую свертку. Однако я запутался, если заполнение нулями также необходимо для деконволюции.

Вопрос
1.Как правильно выполнить деконволюцию на основе теоремы о свертке?
2Почему следующий пример не работает?

Подход
Поскольку A и B измерены, я создал пример для дальнейшегоисследования. Идея состоит в том, чтобы создать B , используя scipy.signal.convolve в режиме same. 

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import convolve
from scipy.fftpack import next_fast_len

# A, in the description above
A = np.array([1, 1, 1, 2, 1, 1])
# The result, I want to get from the deconvolution
kernel = np.array([0, 1, 2, 1, 0, 0]) 
#B, in the description above
B = convolve(kernel, data, mode='same') 

# Using the deconvolution theorem
f_A = np.fft.fft(A)
f_B = np.fft.fft(B)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)

Результат для dk: 

dk = array([ 2.28571429-9.25185854e-18j,  1.28571429+9.25185854e-18j,
       -0.71428571-9.25185854e-18j, -0.71428571+9.25185854e-18j,
        0.28571429-9.25185854e-18j,  1.28571429+9.25185854e-18j])

Ожидаемое значение: 

dk = array([0, 1, 2, 1, 0, 0])

1 Ответ

1 голос
/ 25 октября 2019

Действительно, поскольку ядро ​​1.0 2.0 1.0 центрировано на 2.0 (размытие и раздувание), ширина ядра равна 3. Поскольку массив [1001 * не равен нулю на [0..5], полный извилистый массив paddedB не равно нулю [-1..6]. Тем не менее, функция scipy.signal.convolve(...,'same') возвращает сложный извилистый массив B(0..5)=paddedB(0..5). Следовательно, информация, связанная с paddedB (-1) и paddedB (6), теряется , и восстановление ядра затрудняется , если параметр same из np.convolve() равенused .

Чтобы избежать такой потери информации, вывод paddedB должен быть дополнен, чтобы содержать поддержку извилистого сигнала, вычисленного как сумма Минковского поддержки функции A и поддержки ядра. Опция full из np.convolve() непосредственно вычисляет paddedB без потери информации. 

kernel=[1,2,1]
paddedB = convolve(kernel, A, mode='full')

Для извлечения ядра с использованием теоремы о свертке,входной сигнал A должен дополняться поддержкой функции paddedB

paddedA=np.zeros(paddedB.shape[0])
paddedA[kernel.shape[0]/2: kernel.shape[0]/2+A.shape[0]]=A[:]

# Using the deconvolution theorem
f_A = np.fft.fft(paddedA)
f_B = np.fft.fft(paddedB)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)
# shift to get zero frequency in the middle:
dk=np.fft.fftshift(dk)

Обратите внимание на использование функции np.fft.fftshift() для получения нулевой частоты в середине. 

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import convolve
from scipy.fftpack import next_fast_len

# A, in the description above
A = np.array([1, 1, 1, 2, 1, 1])

kernel=np.asarray([1,2,1])
paddedB = convolve(kernel, A, mode='full')
print paddedB

paddedA=np.zeros(paddedB.shape[0])
paddedA[kernel.shape[0]/2: kernel.shape[0]/2+A.shape[0]]=A[:]
#pad both signal and kernel. Requires the size of the kernel

# Using the deconvolution theorem
f_A = np.fft.fft(paddedA)
f_B = np.fft.fft(paddedB)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)
# shift to get zero abscissa in the middle:
dk=np.fft.fftshift(dk)

print dk

Если получение paddedB невозможно, а B - единственные доступные данные, вы можете попытаться восстановить paddedB, дополнив B нулями, или сгладив последние значения B. Требуется некоторая оценка размера ядра. 

B = convolve(A,kernel, mode='same')
paddedB=np.zeros(A.shape[0]+kernel.shape[0]-1)
paddedB[kernel.shape[0]/2: kernel.shape[0]/2+B.shape[0]]=B[:]
print paddedB

Наконец, окно может быть применено как к paddedA, так и к paddedB, что означает, что значения в середине более значимы, чемядро должно быть оценено. Например, окно Parzen / de la Vallée Poussin: 

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import convolve
from scipy.fftpack import next_fast_len
from scipy.signal import tukey
from scipy.signal import parzen

# A, in the description above
A = np.array([1, 1, 1, 2, 1, 1])

kernel=np.asarray([1,2,1])
paddedB = convolve(kernel, A, mode='full')
print paddedB


B = convolve(A,kernel, mode='same')
estimatedkernelsize=3
paddedB=np.zeros(A.shape[0]+estimatedkernelsize-1)
paddedB[estimatedkernelsize/2: estimatedkernelsize/2+B.shape[0]]=B[:]
print paddedB

paddedA=np.zeros(paddedB.shape[0])
paddedA[estimatedkernelsize/2: estimatedkernelsize/2+A.shape[0]]=A[:]

#applying window
#window=tukey(paddedB.shape[0],alpha=0.1,sym=True) #if longer signals, should be enough.
window=parzen(paddedB.shape[0],sym=True)
windA=np.multiply(paddedA,window)
windB=np.multiply(paddedB,window)


# Using the deconvolution theorem
f_A = np.fft.fft(windA)
f_B = np.fft.fft(windB)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)
# shift to get the zero abscissa in the middle:
dk=np.fft.fftshift(dk)

print dk

Тем не менее, предполагаемое ядро ​​далеко не идеально, так как размер A мал: 

[ 0.08341737-6.93889390e-17j -0.2077029 +0.00000000e+00j
 -0.17500324+0.00000000e+00j  1.18941919-2.77555756e-17j
  2.40994395+6.93889390e-17j  0.66720653+0.00000000e+00j
 -0.15972098+0.00000000e+00j  0.02460791+2.77555756e-17j]
...