How to calculate the correct deconvoluition?

We have 'step type' model data 'M'.

We calculated the convolution of the model data ('C') using the transfer function h(t).

Why can't get correct solution when we deconvolute the convoluted data ('C')?

r equals to C using this code.

If calculated correctly, r equals to M.

Please tell me the correct way to use the function 'deconv'.

close all

%% Step1　Definition of transfer function h(t)

% Time range
n1 = 100;
t = [1:100]';

% Input A and B
A = 5.12
B = 24.8
h = sqrt(A/(pi*(B^2)))*exp(-1/(4*A))*exp(-1*A*((log(t/B).^2)))

%% Step2　Model data 'M'（Solution of deconvolution）

% Model data
M = zeros(1,n1);
M = M';

M(1) = 0;
M(2) = 1;
M(3) = 2;
M(4) = 3;
M(5) = 4;
M(6) = 5;
M(7) = 4;
M(8) = 3;
M(9) = 2;
M(10) = 1;
M(11) = 0;

%% Step3　Model data 'C'（Convolution of 'M'）

% Convoluted data
C = conv(M, h, "full");
C = C(1:n1);

%% Step4　Calcurate deconvolution (Deconvolution of 'C')

% Calcurate deconvolution of 'C'
[q,r] = deconv(C,h);
plot(r)

Answer 1

[Raag]

Hi A.O.,

As per my understanding, the issue you are encountering is due to the numerical instability of the classic ‘deconv’ function when applied to your specific kernel and signal.

• This is a common occurrence, especially when the kernel contains very small values or zeros in its frequency domain representation.
• In such cases, direct deconvolution can amplify noise and computational errors, resulting in extremely large, infinite, or undefined (NaN) values.
• The classic ‘deconv’ function performs direct division in either the time or frequency domain.If the kernel contains values close to zero, the division becomes unstable and can lead to overflow or undefined results.

A more robust approach in these scenarios is to use ‘Wiener deconvolution’ with a regularization parameter. Wiener deconvolution adds a small constant (epsilon) to the denominator, which prevents division by zero and suppresses the amplification of noise. This makes the deconvolution process much more stable and reliable.

Here is how we can do that:

% Define regularization parameter 
epsilon = 1e-8; 
 
% Pad kernel to match signal length 
h_padded = [h; zeros(length(C_full) - length(h), 1)]; 
 
% Perform Wiener deconvolution in frequency domain 
C_fft = fft(C_full); 
H_fft = fft(h_padded); 
Q_wiener_fft = C_fft .* conj(H_fft) ./ (H_fft .* conj(H_fft) + epsilon); 
q_wiener = real(ifft(Q_wiener_fft)); 
 
% Truncate to original signal length 
q_wiener = q_wiener(1:n1);

By introducing a small regularization term (epsilon), the Wiener deconvolution method mitigates the effects of near-zero denominators and suppresses noise, resulting in a stable and accurate recovery of the original signal as you can see below: