Deconvolution results a single value

Hi There!

I am implementing a MATLAB code for the deconvolution of a spectrum. Initially, I am testing my code on the data provided in a published paper.

According to the text mentioned above, by deconvoluting the Fermi function (represented by the dash line, 'Fermi_func' in the MATLAB code), from the spectrum (represented by the black open circles, 'Spectra' in the code), I should obtain the instrumental function, which is represented by a solid line and is a Gaussian function with a width of 0.55 eV.

I digitized the spectrum in the paper (data.txt file attached) and wrote a MATLAB code for its deconvolution (deconvolution.m file attached). However, I am now getting a single value instead of an array for the instrumental function. I was expecting a Gaussian function. I'm not sure what's going wrong.

I am new to MATLAB and any advice regarding this issue would be highly appreciated.

Thank You.

Lakshitha

===============================

Implemented MATLAB code is given below (the .m file is also attached):

clear all;

close all;

clc

A = readmatrix('Data.txt');

Energy = A(:,1);

Spectrum = A(:,2);

Fermi_func = 1-(1./(1+exp((Energy-8.18)/0.025852)));

Fermi_func(1)= 1e-6; % I got an error from deconv that "First coefficient must be non-zero". So I replaced it with a very small value.

Instrumental_Function = deconv(Spectrum,Fermi_func);

% Instrumental_Function =

% deconv(Counts(find(Counts,1):end),Fermi_func(find(Fermi_func,1):end)) %

% This was also implemeted for removing the error of "First coefficient must be non-zero".

display(Instrumental_Function)

f1 = figure('Name','Inverse Photoemission Spectroscopy Fitting');

h1 = plot(Energy, Spectrum,'ko');

hold on;

xlabel('Energy (eV)')

ylabel('Counts (arb. units)')

h2 = plot(Energy, Fermi_func,'k--');

% h3 = plot(Energy, Instrumental_Function,'k-');

hold off;

savefig(f1,'IPES_Fitting.fig')

Answer 1

[John D'Errico]

Convolution of two vectors yields a new vector that is longer than either of them. So if the vectors have length N and M, then the convolved result will have length N+M-1. For example,

X = rand(1,10);
Y = [1 2 1]/4;
Z = conv(X,Y);
numel(Z)

ans = 12

Deconvolution does the opposite. It produces a shorter vector.

Xhat = deconv(Z,Y);
plot(X,'-rx')
hold on
plot(Xhat,'-bo')

As you can see, the deconvolution yields a result that is shorter than the vector Z. The deconvolved result of deconv(Z,Y) will have length

length(Z) - length(Y) + 1

ans = 10

This is as expected. However, IF you perform a deconvolution, where the two vectors have the same lengths, then the result will have length 1. For example...

deconv(randn(1,10),randn(1,10))

ans = -0.0278

So it will be a scalar result. And that makes complete sense. Did you not get a scalar result in what you did?

[Walter Roberson]

Note that it is common to use the "same" or "valid" options for conv() to ignore the tails. But the tails are needed to reconstruct accurately.

[John D'Errico]

Another way of thinking about convolution is in terms of a polynomial multiplication. That is, we can use convolution to multiply the coefficients of two polynomials, with the result being a higher order polynomial. I use this frequently. For example:

syms x
expand((x^2 - x + 2)*(2*x^2 + 3*x + 1))

ans = 

We can emulate that computation with conv:

conv([1 -1 2],[2 3 1])

ans = 1×5

2 1 2 5 2

The coefficient are identical, as they must be.

In the same sense, deconv can be emulated by polynomial division. So

deconv([2 1 2 5 2],[1 -1 2])

ans = 1×3

2 3 1

Divide one polynomial by another and (if we don't worry here about the remainder) the result will be a lower order polynomial. This is expected. So the ratio of two cubic polynomials will yield a lower order quotient. deconv even yields a remainder term, if you take the second argument from deconv.

[Q,R] = deconv([1 2 3 5],[1 1 2 1])

Q = 1

R = 1×4

0 1 1 4

The result is equivalent to this symbolic equality:

Q + dot(R,[x^3,x^2,x,1])/dot([1 1 2 1],[x^3,x^2,x,1]) == dot([1 2 3 5],[x^3,x^2,x,1])/dot([1 1 2 1],[x^3,x^2,x,1])

ans = 

So when we deconvolve two vectors of the same length, we completely expect to see a scalar result for the quotient. This is fully expected.

And Walter's comment about the tails being necessary? This is important. That is, we cannot just pad the one vector with zeros at each end, and expect an accurate result. You should see why that is true if you think about the polynomial division metaphor.

[lakshitha Lathpandura]

Thank you very much, John and Walter, for your feedback regarding my post. I understand now that deconvoluting two functions with the same length results in a single value. John, your explanation was very helpful. Based on your suggestion, I worked for the last 24 hours (without a good sleep :) ), and I came to the conclusion that the 'deconv' function is not going to help me since deconvoluting a function from experimental noise data does not result in a valid function.

So, I started using the 'conv' function to solve it in the other direction, and it works pretty well. However, now I am facing a new issue with the 'lsqcurvefit' function for the convoluted function. I have posted a new question in the link below, and you can see how things went after your suggestions as well.

Thank You.

https://www.mathworks.com/matlabcentral/answers/1907700-using-lsqcurvefit-for-fitting-the-convoluted-function-to-the-data-set

[John D'Errico]

Deconvolution is a difficult problem. It can be ill-posed, which suggests it will be difficult to solve.