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)

[Matt J]

You seem to have assumed that T and Z1 are the same and can therefore be used interchangeably in,

Z2 = df.* T./G;
Z3 = df.* T./Q;

You can easily check, though, that they are not the same. It should really be,

Z2 = df.* Z1./G;
Z3 = df.* Z1./Q;

The reason T and Z1 are not the same is that in the actual continuous-time convolution of g(t) with a step, you get a temperature profile that increases on the interval , but then gradually decays to zero on . Conversely, in your construction of temppad, and hence T, there is no gradual decay. You simply truncate abruptly to zero once is reached.

[Paul]

@Matt J

"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. "

I don't understand that statement.

Let's say we have a causal, finite-duration, discrete-time signal, x[n], i.e., x[n] = 0 for n < 0 and n >= N, that is defined by uniform sampling of a continuous-time signal, for which we want to approximate samples of its Continuous Time Fourier Transform.

We take the DFT of x[n], which only requires the sample vector xs = x[0 <= n <= N-1]. All elements of xs are drawn from x[n] on the positive time axis. Now we have the DFT of x[n] based on what the DFT does in accordance with its definition.

Is that DFT operation equivalent to some other operation that treats the right half of the sample vector of xs as being drawn from the negative time axis of some other signal? There must be another signal involved, since x[n] = 0 for n < 0, but those right half elements of xs are non-zero (or can be, in general). What is this other, equivalent operation?

Can this equivalency be demonstrated with a simple example? I'm genuienly curious.

[Matt J]

@Paul

Can this equivalency be demonstrated with a simple example? I'm genuienly curious.

Sure. Here's an example showing that the FFT of a causal triangle pulse is the same as the Fourier Transform integral (discretized as a Riemann sum) of an anticausal ramp:

%time-frequency discretization parameters
N=100;
dt=2/N;
df=1/N/dt;

tcaus=linspace(0,2,N+1); %causal signal  (triangle pulse)
  tcaus(end)=[];
xcaus=1-abs(tcaus-1);

t=linspace(-1,1,N+1);   %equivalent anticausal signal (ramp)
  t(end)=[];
x=abs(t);

plot(tcaus,xcaus,'r-o'   ,t,x,'b--' ); legend Causal Anticausal
xlabel 'time (sec)'

%Calculate spectrum of causal signal, using FFT
Fcaus=( fftshift(fft(xcaus)))*dt;


%Calculate spectrum of anticausal signal, using Riemann sum
f=t/dt*df;
F=( sum(x.*exp(-2j*pi.*f(:).*t) ,2) ).'*dt;


%Visualize spectra
plot(f , Fcaus ,'r-o'   ,f,F,'b--' ); legend Causal Anticausal

Warning: Imaginary parts of complex X and/or Y arguments ignored.

xlabel 'Freq. (Hz)'; xlim([-10,10]);

●

percentDiff = norm(F-Fcaus)/norm(F)*100

percentDiff = 3.9755e-13

[Paul]

Ok. Now I see where you're coming from. The other (in the context of my previous comment) signal is formed by windowing the central period of the N-periodic extension of x[n] and the other operation is sampling the DTFT of that that other signal.

[Vijayananda]

@Matt J Thank you so much Matt for your insightful answers. regarding the problem of inverse heat conduction, I only have the output T(t) and impulse response g(t). I need to find q(t). so the only way I am going to get q is to use

where the the output exists only upto 1 seconds. when I padd it it abruptly comes to zero without a decay.

And yes. in this code:

Z2 = df.* T./G;

Z3 = df.* T./Q; if I use Z1 instead of T, I get back the impulse response. But in reality, I dont have the convolved T from G and Q. I want to find Q from other two functions. I hope I am not overcomplicating things.

Or maybe I am not understanding it correctly. You said, "The reason T and Z1 are not the same is that in the actual continuous-time convolution of g(t) with a step, you get a temperature profile that increases on the interval , but then gradually decays to zero on . Conversely, in your construction of temppad, and hence T, there is no gradual decay. You simply truncate abruptly to zero once is reached"

But the time domain signals of both Z1 and T looks the same. They extend only upto 1 seconds. I am not able to see any decay. Also analytically, for the step input, the temperetaure keeps on increasing with time.

[Matt J]

@Vijayananda

But the time domain signals of both Z1 and T looks the same.

No, they do not. If you apply the inverse FFT to Z1 without truncating to the interval , you will see that the convolution result has a decaying part in . For good measure, I also demonstrate the agreement of this with time-domain convolution in the code below.

I understand that the interval is not of physical interest for you, but unfortunately you are not free to ignore or change to zero the values that reside there if you still want the Fourier convolution duality theorem to hold. Those values contribute to the Fourier transform of conv(g,q).

To put it another way, recall that the Fourier Transform is originally an integral from to . It only reduces to an integral on a finite interval when the values of the signal are zero outside that interval, but conv(g,q) is not zero outside of . Ultimately, i think you will need to consider non-Fourier deconvolution methods, like I proposed in my other answer.

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


%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;

% Back to time domain
z = (ifft(Z1)); 
%znew = z(1:length(t));

tlong=(0:numel(z)-1) *dt;
convgq=conv(g,q,'full')*dt;

skip=round(numel(z)/30);
plot( tlong , z, 'b-', tlong(1:skip:end),convgq(1:skip:end),'ro')
legend('Fourier Convolution','Time Domain Convolution',Location="northeast")

[Vijayananda]

@Matt J Yes true Matt. Thank you. Yes there is a decay in temperature outside the interval [0 1] for the analytical temperature i used. However, I do not have any measured temperature data outside the interval [0 1] which means I have to use the measured temperature data with zero padding. Even if I take data beyond 1, the temperature keeps going up.

So does this mean, I have to use another method for finding the heat flux ?

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.

[Matt J]

In addition to my suggestions above, here is a non-Fourier method which you might consider. It seems to do a fair job for dt=1e-3. It requires the Optimization Toolbox and that you download,

https://www.mathworks.com/matlabcentral/fileexchange/52795-constrained-minimum-l1-norm-solutions-of-linear-equations

    dt = 1e-3;%sampling period

    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)';


    qrec=deconvolver(DeltaT, g*dt);

Optimal solution found.

    grec=deconvolver(DeltaT, q*dt);

Optimal solution found.

    tiledlayout('vertical')

    nexttile
    skip=round(numel(t)/30);
    plot(t,qrec,'b--',t(1:skip:end),q(1:skip:end),'ro');
    legend( 'Recovered q','True q', Location="east"); axis padded
    ylim([0,1.1]*1e5)

    nexttile
    skip=round(numel(t)/20);
    w=50;
    plot(t,grec,'b--',t([1:w,w:skip:end]),g([1:w,w:skip:end]),'ro');
    legend( 'Recovered g','True q'); axis padded


    nexttile
    f=@(a)a(1:ceil(end/2));
    T1=f(conv(g,qrec,'full'))*dt;
    T2=f(conv(grec,q,'full'))*dt;

    T=DeltaT;
    plot(t,T1,'b--',t,T2,'g--', t(1:skip:end),T(1:skip:end),'ro');
    legend( 'g*qrec','grec*q', 'True Temp', Location="east"); axis padded

    

function [q,G]=deconvolver(y,g)
arguments
    y (:,1)
    g (:,1)
end

    N=numel(y);
    assert(numel(g)==N,'Wrong length')

    [c,r]=deal(g,zeros(N,1));
    r(1)=g(1);

    G = toeplitz(c, r);

    D=diff(speye(N));
    A=[D;-D]; err=repelem(0.1,N-1)';
    b=[err;err];

    Aeq=mean(G,1); beq=mean(y,1);

    q=minL1lin(G,y,A,b,Aeq,beq,0*y,0*y+1e6);

end

[Vijayananda]

Yes @Matt J. The function wont go upto t = 1e-6 since it consumes a lot of memory

[Matt J]

It is not clear to me why you would need dt=1e-6. Your functions are very smooth and slow varying. In any case, try the other suggestions.

[Vijayananda]

I am sorry, the data acquisition is done at 50kHz. so thats why I am looking for a small dt. Thank you som much @Matt J. I will see what I can do about this.

[Matt J]

I am sorry, the data acquisition is done at 50kHz. so thats why I am looking for a small dt.

That doesn't really explain it. Even if you are forced to acquire at 50 kHz, you are free to discard samples, if you don't need them.

[Vijayananda]

@Matt J. I sample at 50 KHz because the phenomena is short lasting. Also, looking at the frequency distribution, its centered around 0 Hz. Can anything be done in this case?

[Matt J]

All real signals have their (magnitude) spectrum is centered at 0.