OK, here are my 2 cents. I'll try to outline the logic and then make a conclusion. We'll start with the Gaussian case.
Suppose that we have a Gaussian peak $f(t)=e^{-t^2/2}$. We want to de-convolve it to the Dirac delta-measure. Since the Fourier transform $\widehat f(s)=\int f(t)e^{-ist}\,dx$ satisfies $\widehat{(f*g)}=\widehat f \widehat {g{,}}$ and the Fourier transform of the standard Gaussian is (up to a constant factor) the standard Gaussian again, the naive but natural idea is just to multiply $\widehat f$ by $e^{s^2/2}$ and to take the inverse Fourier transform of the resulting identically $1$ function.
Unfortunately this is unfeasible as written because of the noise, finite computational precision, etc., so we have to truncate the natural Fourier mulitplier $e^{s^2/2}$ somehow. The idea is to multiply it by some even positive cutoff function $\Psi$ with bounded support. Then the noise will be amplified (on the worst frequency) by
$$
M=\max_s e^{s^2/2}\Psi(s)
$$
but instead of the Dirac point mass, we will get the inverse Fourier transform of $\Psi$. We want that inverse Fourier transform to be as concentrated near the origin as possible. We also want to make it a non-negative function rather than some sign-changing wave. Since we must have $\Psi(0)=1$ to preserve the integral of $f$, we come to the problem of choosing an even function $\Psi$ with bounded support (say $[-\pi,\pi]$; we can always scale later) such that $\Psi(0)=1$, the inverse Fourier transform $\Psi^\vee\ge 0$, and $-\Psi''(0)=\int\Psi^\vee(t)t^2\,dt$ is minimized (other measurements of concentration give slightly different optimizers but pretty much the same results).
This problem is solvable and the solution is the (properly normalized) convolution of $\chi_{[-\pi/2,\pi/2]}(s)\cos s$ with itself, i.e.,
$$
\Psi(s)=\frac 1\pi[(\pi-|s|)\cos s+\sin|s|]
$$
(I'll skip the derivation). Thus, after various scalings, the final family of multipliers is
$$
H_{\alpha,K}(s)=e^{\alpha^2 s^2/2}\Psi(K^{-1}\alpha s)\,.
$$
and the "peak sharpening formula" is $f\mapsto (H_{\alpha,K}\widehat f)^\vee$.
It turns out that $K=1$ does not noticeably change the picture (but filters out the high frequency noise from the measurements). The noise amplification grows pretty fast with $K$. The corresponding maxima are about
$$
\begin{matrix}
K &\max
\\
1.0 & 1.09
\\
1.1 & 1.55
\\
1.2 & 2.65
\\
1.3 & 5.33
\\
1.4 & 12.4
\\
1.5 & 33.2
\end{matrix}
$$
(with noise at $1$ to $10\%$ of the signal, the higher values of $K$ are pretty useless), so I would recommend setting $K$ between $1.3$ and $1.4$. The whole game is choosing $\alpha$. The legitimate choice is up to the width of the narrowest peak but somewhat higher values are also sort of admissible except you'll get peaks exchanging energy a bit. The rule of thumb is to watch for noticeable negative values in the transformed signal: once they appear, you certainly went too far. Before that moment you are reasonably safe (but not entirely foolproof).
The only difference for EMG is that you want to add one more factor $1+i\beta s$ to your multiplier, i.e., to play with the family
$$
H_{\alpha,\beta, K}(s)=e^{\alpha^2 s^2/2}(1+i\beta s)\Psi(K^{-1}\alpha s)\,.
$$
If $\beta>0$ gets large, it may force you to make $K$ smaller. The effect of going too far in $\beta$ is the same: you'll see noticeable negative values in the transform plus the peaks will noticeably shift to the left. Playing with $2$ parameters is harder than with one and I have no clear advice on in which order to adjust them to get the best results. Just try and let me know what you think.
.
The original signal is in green, the dots are the actual locations at which EMG's are "centered". Something like that is possible IF you can guess one particular parameter right. I still cannot figure out how to teach the machine that guessing.
You'll, probably, like this one too:
Or this one:
The current version of the code. This is a dirty homemade contraption, of course, but it is totally automatic.
import graph;
access settings;
//settings.outformat="png";
size(500,300,IgnoreAspect);
int sec=seconds();
srand(sec); //Just generates a new picture every time
//srand(1562245957); //was really bad, fixed
//srand(1562264140); //was bad, fixed
//srand(1562265837); //still bad
bool ON=false; //Set to true if you want to see suppressed low values
real beta=0.04; //Noise level (maximum; you need to adjust the code a bit to use the average)
real Beta=0.12; //Cutoff level (approximately 3 times maximal noise)
real K=1.7, L=3, TOL=2.5; //Don't change these UNLESS
int res=8, dirr=5;//you understand what you are doing
int M=4; //The number of peaks. Set to 1 to see the individual response.
int n=4096; //The sampling range (4096 time ticks with no signal before and after)
real h=1/n;
real QQ=0.07*unitrand(); //Average lambda; set to 0 to see pure Gaussians
real AA=0.15*sqrt(unitrand()); //Average sigma; set to zero to see pure exponentials
real det=unitrand(); //Peak spacing from equidistant det=1 to random det=0
real[] A,Q,X,Y;
for(int m=0;m<M;++m)
{
A[m]=AA*(0.2+unitrand());
Q[m]=QQ*(0.5+0.5*unitrand());
X[m]=det*m/M*n/2+(1-det)*(rand()%quotient(n,2))+quotient(n,8);
Y[m]=0.45+0.55*unitrand();
}
X=sort(X); A=sort(A); //sigma grows with time
real qual(real x)
//deconvolution quality control: the bigger the sum of qual(p[k]), the better.
{
return
max(x-2*Beta,0)^3-min(x,0)*Beta^2;
}
//write(A); write(Q); pause();
pair[] p, ptemp;
//p is the (simulated) raw data. It would be interesting to try it on something real.
////////////////Generation of p /////////////////
for(int k=0;k<n;++k) {p[k]=(0,0); ptemp[k]=(0,0);}
p.cyclic=true; ptemp.cyclic=true;
for(int m=0;m<M;++m)
{
pair f(int k)
{
return exp(-pi*A[m]^2*k^2)*expi(-2*pi*k*X[m]*h)/(1,2*pi*Q[m]*k);
}
for(int k=-floor(n/2);k<n/2;++k) ptemp[k]=f(k);
ptemp=fft(ptemp);
ptemp/=max(map(abs,ptemp));
ptemp.cyclic=true;
p+=Y[m]*ptemp;
}
real s=0;
for(int k=0;k<n;++k) {s+=p[k].x;}
for(int k=0;k<n;++k) p[k]+=beta*(2*unitrand()-1);
//////////// End of generation of p ///////////////////
//////////// Algorithm begins ///////////////////
pair[] g=fft(p)/n;
real psi(real t)
{
real s=abs(t);
if(s>pi) return 0;
else
//return sin(s)/s*(1-s/pi);
return ((pi-s)*cos(s)+sin(s))/pi;
}
pair R(pair x)
//Suppression of low values.
{
if (ON) return x;
real y=max(x.x-Beta,0);
if(x.x<0) return (0,0);
else return ((x.x<4*Beta? y*4/3:x.x),x.y)/(1-4*Beta^2);
}
real expp(real x)
{
return 1+x+x^2/2+x^3/12;//A bit more stable than the true exp
return exp(min(x,15));
}
struct T {pair[] p; real u; real q;}
int count=0;
real[] pM; pM[n]=0;
for(int k=n-1; k>-1; --k) pM[k]=max(p[k].x,pM[k+1]);
pair[] pp=copy(p), p1;
T U(pair D, real qq)
{
++count;
write("U "+string(count));
real u=0,v=1;
void proc()
{
for(int kk=0; kk<10;++kk)
{
real w=(u+v)/2;
real a=w*D.x, q=w*D.y;
real b=min(max(a*sqrt(2*pi)/K,2*pi*q/L),1/5);
pair[] G=fft(pp)/n, mult;
p1=copy(pp);
for(int k=0;k<n;++k) mult[k]=(0,0);
mult.cyclic=true;
G.cyclic=true;
for(int k=floor(-n/2);k<n/2;++k)
{
mult[k]=expp(pi*a^2*k^2)*(1,-2*pi*q*k)*psi(b*k);
if(abs(mult[k])>TOL) mult[k]=TOL*unit(mult[k]);
G[k]*=mult[k];
}
pair[] ppp=fft(G); ppp=reverse(ppp);
ppp.cyclic=true;
bool flag=true;
//for(int k=0;k<n;++k) if(abs(mult[k])>TOL) flag=false;
for(int k=floor(n*qq-n/res);k<n+n/2/res && flag;++k)
{
if(
ppp[k].x<-0.8*Beta || (k>-1 && k<n && pM[k]<Beta/2 && ppp[k].x>Beta)
)
flag=false;
}
if(flag && abs(G[floor(n/2)-1])<1/10^10) {u=w; p1=copy(ppp);} else v=w;
}
}
proc();
p1.cyclic=true;
if(p1[n-1].x>Beta) {write(p1[n-1].x,pM[n-1]); pause();}
real q=0;
for(int k=floor(n*qq-n/res);k<n;++k) q+=qual(p1[k].x);
T tt=new T;
tt.p=copy(p1); tt.u=u; tt.q=q;
write(tt.u);
return tt;
}
T[] tt;
for(int k=0; k<res; ++k)
{
draw(box((k/res,0),((k+1)/res,1)),yellow);
pair D=expi(pi/2);
tt[k]=U(D,k/res); real W=tt[k].q;
for(int m=dirr-1; m>-1; --m)
{
D=dir(90*m/dirr);
T ttt=U(D,k/res);
if(ttt.q>W) {W=ttt.q; tt[k]=ttt;}
}
//write (tt[k].u, tt[k].q);
pp=copy(tt[k].p);
}
for(int k=0;k<n;++k)
{
int q=floor(k/n*res);
if(q<1) {pp[k]=tt[0].p[k];}
else
{
real t=q+1-k/n*res;// t=sin(pi*t/2);
pp[k]=((1-t)*tt[q].p[k]+t*tt[q-1].p[k]);
}
}
pp=map(R,pp);
real ss=0;
for(int k=0;k<n;++k) {ss+=pp[k].x;}
write("*************");
path Q=nullpath, P=nullpath;
for(int k=0;k<n;k+=4)
{
Q=Q--(k*h,pp[k].x);
P=P--(k*h,p[k].x);
}
draw(P,green);
draw(Q);
draw((0,1)--(0,0)--(1,0));
draw(box((0,-Beta),(1,Beta)),red);
for(int m=0; m<M ;++m)
{
dot((X[m]*h,0));
pair Z=intersectionpoint((X[m]*h,-20*M)--(X[m]*h,20*M),Q),ZZ=(X[m]*h,Y[m]);
draw((X[m]*h,0)--Z,blue); dot(ZZ,magenta);
}
label(string(sec),(0,0),S);
label("Area ratio="+string(floor(1000*ss/s)/1000),(0.8,1.13),(0,0),fontsize(8));
//pause();
