Matching spectral data to known spectral lines

Preamble:

I designed some code to help me take spectral images from a spectrometer and convert pixel number into wavelength. It does reasonably well at matching a calibration spectrum of a neon lamp to known spectral lines of neon, but only if I make sure to include all actual neon peaks from the calibration spectrum. Sometimes it's hard to tell if a peak in the spectrum is just noise/contamination or if it's real. If I exclude real peaks, my algorithm fails miserably. It's not a huge problem, but it'd be neat to automate this even further.

The Problem:

The essence of the problem is this: I have a series of known spectral line of neon. I take a spectrum of a neon lamp, and then want to match the lines I see in the image up to the known lines. However, I may not have exactly the same number of calibration lines as I do reference lines.

My reference data is this:

refdat = {{585.25, 27}, {588.19, 11}, {594.48, 29}, {597.55,
7}, {603.00, 6}, {607.43, 20}, {609.62, 35}, {614.31,
49.5}, {616.36, 16}, {621.73, 10.5}, {626.65, 24}, {630.48,
11}, {633.44, 29}, {638.30, 43}, {640.22, 78.5}, {650.65,
60}, {653.29, 19}, {659.90, 25.5}, {667.83, 39.5}, {671.70,
30}, {692.95, 35.5}, {703.24, 100}, {717.39, 4.5}, {724.52, 50}, {743.89, 6}};

The x-values (wavelengths in nm) are the most reliable, so I would prefer to use these only. The intensity response (y-values) are pretty non-linear so I'm not sure they're of much use.

A good set of data looks like this:

good = {{11, 0.1032706710176553}, {44, 0.1146552044370179}, {115,
0.2113968543791952}, {150, 0.05051515444593649},
{213, 0.052753899062279425}, {264, 0.13868347448613108}, {290,
0.22283763808371357},
{344, 0.4252736853859361}, {368, 0.12755375457964316}, {430,
0.10566687181148425},
{487, 0.2319607538664892}, {531, 0.08104101568883736}, {565,
0.2742127219013794},
{622, 0.4568371231579498}, {644, 0.7877989366604489}, {765,
0.36850130722650787},
{795, 0.13327028393091367}, {872, 0.1449608891795199}, {964,
0.20420376201460919},
{1009, 0.13586396031680412}, {1256, 0.18521400598076912}, {1376,
1.}} and a "bad" set of data might look like this:

{{11, 0.1032706710176553}, {44, 0.1146552044370179}, {115,
0.2113968543791952}, {264, 0.13868347448613108},
{290, 0.22283763808371357}, {344, 0.4252736853859361}, {368,
0.12755375457964316},
{430, 0.10566687181148425}, {487, 0.2319607538664892}, {565,
0.2742127219013794},
{622, 0.4568371231579498}, {644, 0.7877989366604489}, {765,
0.36850130722650787},
{795, 0.13327028393091367}, {872, 0.1449608891795199}, {964,
0.20420376201460919},
{1009, 0.13586396031680412}, {1256, 0.18521400598076912}, {1376,
1.}} Basically, the threshold on my FindPeaks call just missed out on 3 peaks that were real. My algorithm is just using FindMinimum to do a linear fit of the pixel values to the reference data as in scale * pixel_values + offset:

results = Table[
FindMinimum[
Mean[(scale*pks[[All, 1]] - offset -
refdat[[i ;; i + Length@pks - 1, 1]])^2],
{{scale, 22.8}, {offset, 617.4}}],
{i, Length@refdat - Length@pks + 1}];

Where pks is either good or bad. If I only have 22 peaks, but refdat has 25 points, I naively try refdat[[1;;22]], refdat[[2;;23]], refdat[[3;;24]], and refdat[[4;;25]]. This is fine if I happen to have selected every real peak but my spectrometer only collected the first or last $$n$$ peaks; it fails if I miss one in the middle somewhere.

Here is what happens if I have "good" data (it's pretty clear in this data what is a peak and what isn't, but that's not always the case): But if I miss those 3 points: My Question:

Is there a reasonably efficient algorithm for linearly scaling one vector to match another with (a few; maybe 2-4) possible missed values?

Simply trying every combination of values seems prohibitive. My current code isn't too bad as I can usually pick out all the real values after a couple of tries with moving the FindPeaks threshold, though not always. As a last resort I can always pick out values manually, but it would be neat to know if there's a better solution. This seems like a problem that would probably be well-studied from a mathematical perspective, but I'm not well-versed enough to find it.

If most points in refdat do have corresponding points in good then you should be able to find a stretch of $$n$$ points that exist in both sets. Let's say that $$n=5$$. You can then find the scaling like this:

ref = Partition[First /@ refdat, 5, 1];
data = Partition[First /@ good, 5, 1];

findLinearModel[ref_, data_] := LinearModelFit[Transpose[{data, ref}], x, x]
findLinearModelError[ref_, data_] := Module[{lm},
lm = findLinearModel[ref, data];
Norm[ref - lm /@ data]
]

errors = Outer[findLinearModelError, ref, data, 1];
minPos = Position[errors, Min[errors]];

lms = Outer[findLinearModel, ref, data, 1];
lm = First@Extract[lms, minPos];

lines = InfiniteLine[{#, 0}, {0, 1}] & /@ (lm@*First /@ good);
intersections = Flatten[{x, y} /. Quiet@Solve[{x, y} ∈ #, {x, y} ∈ Line[refdat]] & /@ lines, 1];

ListLinePlot[
refdat,
Epilog -> {
Red, PointSize[Large], Point[intersections]
}] With this estimate based on five points, you can now exclude points you don't have correspondences for and find a new transform based on all of the ones you have, if you want.

x1 = refdat[[All, 1]];
x2 = lm /@ good[[All, 1]];
filteredRefDat = First /@ Nearest[x1, x2];
lm = LinearModelFit[Transpose[{First /@ good, filteredRefDat}], x, x];

lines = InfiniteLine[{#, 0}, {0, 1}] & /@ (lm@*First /@ good);
intersections = Flatten[{x, y} /. Quiet@Solve[{x, y} ∈ #, {x, y} ∈ Line[refdat]] & /@ lines, 1];

ListLinePlot[
refdat,
Epilog -> {
Red, PointSize[Large], Point[intersections]
}] • Thanks, this works pretty well! It's usually (though not always) possible to get 4 or 5 peaks in a row. I changed x1 and x2 in your filteredRefDat to refdat[[All, 1]] and intersection[[All, 1]] - I assume that was what was supposed to go there. – MassDefect Aug 17 at 19:20
• @MassDefect Thanks, updated to fix that. – C. E. Aug 17 at 19:31