Deconvolution using FFT - a classical problem

Hello friends, I am new to signal processing and I am trying to achive deconvolution using FFT. I have an input step function u(t) applied to an impulse response given by . The output function is . I am trying to convolve g and u to get y as well as deconvolve y and g to get u. However, I quite cannot get the right answers. I understand that the deconvolution process is ill-posed and I have to use some kind of normalization process but I am lost. I also apply zero padding to twice the length of the input signals. Any sort of guidance will be appreciated.

After using deconvolution in the fourier domain:

Y = fft(y)

G = fft(g)

X = Y./G

x = ifft(X)

I am getting an output shown below:

Which is not the expected outcome. Can someone shead light on what is happening here? Thank you.

Answer 1

[Matt J]

dt=0.001;
N=20/dt;

t= ( (0:N-1)-ceil((N-1)/2) )*dt; %t-axis

u=(t>=0);
g=3*exp(-t).*u;
y=conv(g,u,'same')*dt;

Y = fft(y);
G = fft(g);
X = Y./G;
x = fftshift(ifft(X,'symmetric')/dt);

figure;
sub=1:0.3/dt:N;
plot(t,3*(1-exp(-t)).*u,'r.' , t(sub), y(sub),'-bo'); 
xlabel t
legend Theoretical Numeric Location northwest
title 'Output y(t)'

figure;
plot(t, u,'r.' , t(sub), x(sub),'-bo'); ylim([-1,4])
title Deconvolution
xlabel t
legend Theoretical Numeric Location northwest

[Vijayananda]

Thank you so much Matt. This was insightfull. When I just change the time axis to t= (0:N-1)*dt this, my whole result changes. I want the time to start from 0 since negative time is not physical. However, the results should not vary right? What could be the reason?I changed nothing else other than the time.

I also pad my data to twice the length of the array and time interval dt is aorund e-6. Will this create any issues?

[Matt J]

FFTs and IFFTs inherently integrate over both positive and negative times and frequencies. So, no, you cannot define the t-axis arbitrarily. But you can just delete the unwanted time points from the final result if you want.

As an aside, it is strange that you would be processing causal signals with FFTs. Usually FFTs are for non-causal, finite duration signal processing, whereas infinite, causal signals are usually processed with z-transforms.

[Vijayananda]

You are right. I am trying to deconvolve a heat transfer problem. The time is finite from 0 to 1 with an interval of 1 microseconds. I also zeropad the impulse response and the output data to twice their length to negate the periodicity of teh FFT.

[Matt J]

I also zeropad the impulse response and the output data to twice their length to negate the periodicity of teh FFT.

Zero-pad the impulse response and the input data, I think you meant to say.

That is effectively the same as my original answer, although I am zero-padding on the left, instead of the right. In other words, you can interpret the zeros for t<=0 that you say you don't want as an alternative implementation of zero-padding.

[Paul]

@Vijayananda

If you show you're code, it will be easier for people to help with whatever specific problem you're having.

@Matt J

"FFTs and IFFTs inherently integrate over both positive and negative times and frequencies."

I'm curious about this statement.

Did you mean "integrate" in the context of summation?

It seems to me that the DFT and IDFT, which are implemented in Matlab with fft and ifft respectively, are inherently defined only for positive times and positive frequencies. The DFT is (typically? usually? certainly in fft) defined as the summation over x[n], n = 0:N-1, where n is the discrete-time independent variable. The IDFT is defined as the summation over X[k], k = 0:N-1, where 2*pi*k/N are the discrete "frequencies" (which are actually angles). Postive times and positive "frequencies" respectively. Of course, if a finite-duration discrete-time signal is non-zero over an interval that includes n < 0, we can map it to another signal defined over n = 0:N-1, take the DFT of the mapped signal, and then manipulate that DFT (if necessary, depending on how we do the manipulation) to get a frequency domain representation of the original signal. But that manipulation is needed precisely because the DFT is only defined/implemented for signals that are non-zero for positive time. An analgous statement could be made regarding the IDFT.

"Usually FFTs are for non-causal, finite duration signal processing, whereas infinite, causal signals are usually processed with z-transforms."

See above regarding "FFTs are for non-causal" (agree on the finite-duration part)

Infinite-duration* signals can be processed (analyzed?) with the Discrete Time Fourier Transform or the z-transform, both causal and non-causal.

* More precisely, certain subsets of infinite-duration signals.

[Matt J]

Did you mean "integrate" in the context of summation?

What I really meant was, when you use a DFT to discretize a continuous Fourier Transform integral of some finite duration signal s(t), what the DFT does is always equivalent to treating the right half of the sample vector as drawn from the negative time axis, and the left half as drawn from the positive axis. Shifting the sampling interval in the continuous domain or redefining where t=0 is won't change that.

[Vijayananda]

Hello friends @Paul, I am trying to solve an inverse heat transfer problem. In a semi-infinite medium which is supplied with a constant heat flux , the temperature distribution is given by:

This is an LTI system. The impulse response of the system is given by:

and the constant heat flux in the step form is q" given the signal starts from t = 0.

So we can write:

In fourier domain:

Z = G*Q

I have analytical forms for all three functions as shown below.

and

I am getting temperature by convolving g and q. No issues. However, when I try to get either g, or q back from the other signals, I am running into error. Division in the fourier domain is ill posed and there are zeros in the denominator. I am not sure how to rectify these. If you add noise to the signals, solving becomes more difficult. The code is attached to this comment.

clc
clearvars
dt = 1e-6;%sampling period
df = 1/dt;%sampling frequency
t = 0:dt:1-dt;%Time vector
qmax = 1e5;%W/m2 maximum step heat flux
alpha = 3.56E-6;%m2/s diffusivity
k = 15; %W/m-k thermal conductivity
x = 1e-5;%m junction depth
DeltaT = (((2*qmax)/k).*sqrt((alpha*t)/pi).*exp((-x*x)./(4*alpha*t)))-(((qmax*x)/k).*(1-erf(x./(2*sqrt(alpha*t)))));
g = sqrt(alpha./(pi*k*k.*t)).*exp((-x*x)./(4*alpha.*t));
g(1) = 1e-9;%to replace NaN at t = 0
q = qmax*ones(length(t),1)';
n = length(DeltaT);
m = length(g);

%Analytical function plots
% figure
% subplot(3,1,1)
% plot(t,DeltaT,LineWidth=2)
% title("Temperature")
% subplot(3,1,2)
% plot(t,g,LineWidth=2)
% title("Impulse response")
% subplot(3,1,3)
% plot(t,q,LineWidth=2)
% title("Heat flux")

%Zero padding to get linear convolution from circular convolution
qpad = paddata(q,n+m-1);%padded heat flux
gpad = paddata(g,n+m-1);%padded impulse response
temppad = paddata(DeltaT,n+m-1);%padded temperature

% Perform fourier transform
G = fft(gpad);
T = fft(temppad);
Q = fft(qpad);

% Convolve heat flux and impulse response
Z1 = dt.* G.*Q;

% Convolve temperature and impulse response
Z2 = df.* T./G;

% Convolve temperature and heat flux
Z3 = df.* T./Q;

% Back to time domain
z = (ifft(Z3)); %change Z3 to Z2 and Z1 for other outputs
znew = z(1:length(t));

plot(t,znew)

Answer 2

[Matt J]

Since you are trying to deconvolve in the presence of noise, it would make sense to use a regularized deconvolver like deconvreg or deconvwnr. These are from the Image Processing Toolbox, but there is no reason they wouldn't apply to one-dimensional signals as well.