Deconvolution of signal 1 from a known signal 2

Could you please help deconvoluting a signal 1 (a real signal) from a signal 2 (an intrinsic response of the detector). Please note the times of signals are different, and the response intensity can be set a variable in the deconvolution.

data: https://docs.google.com/spreadsheets/d/1PoDzJcYLgEiTNHuzjy5fwZdI21XJTNKlI6HNE6NH4HE/edit?usp=sharing

clear all

close all

Data = xlsread('signals.xlsx','Sheet1','A1:D25000');

% response of the detector (signal to be deconvoluted from the real signal), also called impulse response

Time_resp = Data(:,1);

response = Data(:,2);

% real signal

Time_sig = Data(:,3);

signal = Data(:,4);

figure

plot(Time_resp,response)

hold on

plot(Time_sig,signal)

The convolution should only starts from time = 0, so some shifting or rescalling needs to be done.

Many thanks in advance.

[Paul]

@ChemPhysicist

The blue curve is the impulse response, h, the red curve is the output, y, and the goal is fiind the input, u, such that

y = h cxc u, where cxc is the convolution operator. Is that correct?

If so, what should be assumed for the output for t > 45 (or so). Does it drop to zero instantananeously? Decay to zero with a time constant? Something else?

[ChemPhysicist]

@Star Strider

If we get a function that reproduces the impulse function, then we can use the number of data points we need .. that could be an alternative to fix the sampling frequency?

I tried to fit the impulse function with this general function that looks applicable to its shape:

x = linspace(0,8,800);

mdl = fittype('a*exp(-b*x) .* sin(c*x)','indep','x')

fittedmdl = fit(Time_resp,response,mdl)

and got the following parameters. It reproduces quite good the function with these outputs:

fittedmdl =

General model:

fittedmdl(x) = a*exp(-b*x) .* sin(c*x)

Coefficients (with 95% confidence bounds):

a = 0.9318 (0.929, 0.9347)

b = 1.359 (1.354, 1.364)

c = 2.725 (2.722, 2.728)

We can therefore use this function instead of the response signal. Could you please retry with this function and these outputs and use the number of data points required for having the right sampling frequency.

Answer 1

[Paul]

I'm pretty sure that any DFT (fft()) based method will have limitiations, but we can try as follows.

Assume that the system model is

y(t) = convolutionintegral(h(t),u(t))

We'll use samples of y(t) and h(t) and approximate the convolution integral with a convolution sum.

Read in the data (readtable command courtesy of @Star Strider)

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1103615/Signals.xlsx', 'VariableNamingRule','preserve');

Assign the time and response vectors for the impulse response

th = T1{:,1}.';
hval = T1{:,2}.';

Trim off the leading data for t < 0 and compute a nice time vector

hval = hval(th >= -1/10000);
th = (0:numel(hval)-1)/2500;

Do the same for the output

ty = T1{:,3}.';
yval = T1{:,4}.';
yval = yval(ty >= -1/10000);
ty = (0:numel(yval)-1)/500;

Plot the impulse response and the output

figure;
plot(th,hval,ty,yval)
legend('h','y')

Interpolate the output to the sampling period common of the impulse response (could also try decimating the impulse repsponse to the sampline period of the output)

dt = 1/2500;
yval5  = interp(yval,5);
ty5 = (0:numel(yval5)-1)*dt;

Assume that the impulse response is identically zero for t > th(end) and that the input is the same length as the output. The latter assumption seems resasonable for a physical process. Also, the output doesn't display any serious transient behavior towards the end of the data. We can try to solve for the input using the relationship between linear convolution and the DFT. However, we have a problem insofar as that if y = conv(h,u), then numel(y) > numel(u), which contradicts an assumption. I think the result of this contradiction will be be apparent below.

uval = ifft(fft(yval5,numel(yval5)+numel(hval)-1) ./ fft(hval,numel(yval5)+numel(hval)-1))/dt;
tu = (0:numel(uval)-1)*dt;

uval should be real

max(abs(imag(uval)))

ans = 0

Plot the reconstructed input. It looks like noise. But this shouldn't be too surprising becasuse of the noise on the output and our assumption that output is the response to only the input.

figure
plot(tu,uval)

Check to see how closely the predicted output based on the reconstructed input matches the actual output

ypredict = conv(hval,uval)*dt;
typredict = (0:numel(ypredict)-1)*dt;
ypredict(typredict > ty(end)) = [];
typredict(typredict > ty(end)) = [];

plot(ty,yval,typredict,ypredict)
legend('actual','predicted')

Not too bad after the initial transient, which I'm nearly certain arises from the contradiction discussed above.

Zoom in after the transient

figure
plot(ty,yval,typredict,ypredict);
legend('actual','predicted')
xlim([10 50])

The analysis would be different if we instead assumed a system model of the form

y(t) = convolutionintegral(h(t),u(t)) + n(t)

If interested, there may be an alternative approach that gets a better answer.

[ChemPhysicist]

@Paul

I really like your approach. The problem indeed is that the first part of the predicted signal is so bad, and in fact it's mainly the part that I wamt to solve in order to perform a good fitting on my signal .. and that's why I have measured the response (impulse) in order to get the information on that spike and remove it.

I think it seems not an easy task to do this. I therefore just decided to perform a fitting on the measured signal with my physical model of 3 exponential terms, then used the fit parameters to build a function that fit the signal (but does not fit the impulse), and I used this function in a convolution with the impulse to reproduce the measured signal. It gave a result not too bad but needs a bit more work.

mdl = fittype(' -c*exp(-5.7E+6*x) -e*exp(-f*x)+ (c+e)*exp(-b*x) ','indep','x')

fittedmdl = fit(x,y,mdl)

coefficientValues = coeffvalues(fittedmdl);

Fittfunction = -coefficientValues(2)*exp(-5.7E+6*x) -coefficientValues(3)*exp(-coefficientValues(4)*x)+ (coefficientValues(2)+coefficientValues(3))*exp(-coefficientValues(1)*x);

figure

plot(x,Fittfunction)

hold on

plot(x,y,'b.')

having the fit function and the impulse function:

Convolute the two functions (signals):

As you can see, this did a bit the job .. but if instead we consider that the impulse is an "undesirable" signal which has no physical meaning rather than a a signal/ringing which adds up to the real signal, in this case the total signal which is measured is the sum of the signal coming from the process we study + the response signal. Here we can just add the impulse to the fit to reproduce well the measured signal:

The pure sugnal is just a subtraction of the raw signal minus the response signal:

[Paul]

@ChemPhysicist

I'm now confused about the end goal. Based on a comment to another now-deleted Answer, I thought the goal was to find the input signal, u(t), that generates an output signal, y(t), based on an impulse response, h(t), where h(t) and y(t) are given (or measured). Not sure if there was going to be a second step after u(t) was found, assuming it could be (even approximately).

Now it looks like the real goal is to find something called the "pure signal," but I don't know what that term means. As close as I can tell, the red curve in the last plot in the comment above is just the raw signal minus the impulse response. I'm not sure about the mathematical justification for that operation. To be sure, for an LTI system, the output due to any input is a weighted sum of the (generally complex) exponentials that form the impulse response and the input. Maybe the goal is to subtract out the contribution of the weighted impulse response exponentials to get the output due to just the input? Of course, one can always subtract the impulse response from the output; the mathematical meaning of the result is murky, at best, at least to me.

A few other thoughts:

It looks like the raw signal in the comment is different that the raw signal in the original Question (and in the attached spreadsheet)?

Not sure what the intent is of convolving the impulse with the (fitted) output. I mean, the output is already formed from the convolution of the impulse response and the input, so convolving the output again with the impulse response is just an additional filtering operation on the output. In this case, the impulse response has a low pass characteristic, which results in that transient in the third plot because the fitted output is discontinuous at t = 0. Of course, any other low-pass filter could have been used.

Is the discontinuity at t = 0 in the fitted output realistic? That implies the input inlcudes the second derivative of the Dirac delta function (assuming the impulse response is of the form h(t) = exp(-a*t)*sin(w*t), which I thought was already shown elsewhere in this thread). Is that a reasonable expectation for the input?

Finally, I attacked the problem a different way that yieled an input signal, u(t), that resulted in ypredicted being a much better match to y(t) for that initial transient, but now I don't know if that result is useful. If it is, I can post it was well.