Re: [Gnuplot-info] sin(x)/x smoothing option?

[Reginald B. <pul...@ya...>]

For regularly sampled data this is trivial.  There is an excellent sinc(x) interpolator in the CWP/SU corpus which is BSD licensed that Dave Hale wrote.  I supported both CWP/SU and SLT/SU for many years, so I know the code base very well.    I'll extract that and remove any gratuitous dependencies later this week. I *wasn't * busy when I posted the question and was working on a project when suddenly I had higher priority tasks confront me.

Hale's interpolator is zero phase as that is what is used in seismic processing.  However, for many applications such as electrical engineering I think one needs a minimum phase interpolator.  I need to study the subject a bit and verify that is what is being done in digital storage oscilloscopes.  However, that's a relatively simple change to Hale's code.

Irregularly sampled data requires a different approach.  There is a relatively new (2004-2008) development in applied mathematics called "compressive sensing" developed by Emmanuel Candes and David Donoho.  The mathematical justification is *very* complex.  However, it's quite easy to use and is extremely general.  I refer to the field as "sparse L1 pursuits".  Applications include the "Netflix prize problem", predict what someone would like based on previous selections, inverse problems where you know the form of an equation that fits the data, but don't know the coefficients, "blind source separation" in which from a small number of recordings you select a single conversation in a room full of people all talking at the same time and a slew of other things.  It's a general solution to Ax=y under the constraint that most of the terms in x are zero and the columns of A are incoherent.  Mathematically the constraint is called the "restricted isometry property", any two combinations of columns has near zero cross correlation.  That is NP hard to solve, so in practice one attempts to solve Ax=y and if you get an answer it's been proven to be the optimal L0 solution, but obtained in L1 time.

The Marquardt-Levenberger L2 solver in gnuplot was how I always did things until I stumbled across sparse L1 pursuits.  I did a lot of data analysis with gnuplot using that.  Addling sparse L1 pursuits would make gnuplot a *very* powerful data analysis tool.

In "Numerical Methods that Work" Foreman Acton has sections on problems not to try to solve because they are illposed.  One of these is solving d = a*exp(A*t) + b*exp(B*t) which is the form of the radioactive decay of two isotopes.  Acton's advice is to send the problem to a chemistry lab to get A & B.  Sparse L1 pursuit will solve for a, A, b and B without any problems.

I had been solving inverse problems in fluid flow in porous media with gnuplot which struggled as it is a more complex variant of  the problem referred to by Acton.  Following Mallat's 3rd ed of "A Signal Processing Tour of Wavelets" I decided to try basis pursuit.  I was quite stunned by the results.  And then I remembered Acton's advice to the reader.  So I got a copy of "A Mathematical Introduction to Compressive Sensing" by Foucart and Rauhut and began a 3 year odyssey into the most complex and bizarre mathematical logic I have ever encountered.

So to summarize:

I am going to implement both sin(x)/x interpolation and sparse L1 pursuit inversions for gnuplot.  I've used gnuplot for about 25 years and consider it one of the best pieces of software I have ever used.  So I think it high time I made a significant contribution.

I don't know the code base for gnuplot at all and am only marginally familiar with the GLPK code which is needed for sparse L1 pursuits.  Normally I use the GMPL languagewhen using GPLK.  However, I've maintained several million lines of other people's code and know all the relevant background material very well.  So it's just a matter of putting in the time.  Sinc(x) interpolation is a couple of days work, so I'll get that done fairly quickly.  Sparse L1 pursuits are quite a bit more work, probably several weeks full time, so that will take a couple of months to complete as being "retired" has made me busier than I ever was before.