Deconvolution of signal 1 from a known signal 2

25 views (last 30 days)
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.
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);
hold on
The convolution should only starts from time = 0, so some shifting or rescalling needs to be done.
Many thanks in advance.
ChemPhysicist on 21 Aug 2022
Edited: ChemPhysicist on 21 Aug 2022
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.

Sign in to comment.

Answers (1)

Paul on 24 Aug 2022
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('', '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
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
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.
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)) = [];
Not too bad after the initial transient, which I'm nearly certain arises from the contradiction discussed above.
Zoom in after the transient
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.
Paul on 24 Aug 2022
Edited: Paul on 25 Aug 2022
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.

Sign in to comment.




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!