Shifting Peaks in Curve Fitting
This blog post describes algorithm(s) to handle shifting peaks in signal separation and curve fitting problems.
All the related code can be found from my GitHub repo.
Background
The basic problem of signal separation problem here can be formulated as follows: given the pure component signals and mixture signal that is a linear combination of the pure component signals, one needs to solve the contributions of individual components. This problem appears in many different signal processing applications.
In many cases, the problem is easy to solve using good old classical least squares fit method. The method works pretty well in many cases. To be more specific, it can handle random amplitude noise well. The problems appear when the mixture is not a linear combination of pure components (when it was assumed to be so), the mixture signal contains unknown components or there is error also in independent variable (also called x).
This blog post concentrates on one type of problem where independent variable contains an error that is systematic by nature. This means that errors change from one sample to sample another but they follow some certain model.
This blog post presents a method to handle this kind of errors. This method requires the user to specify a model for x-axis errors. In the ideal case, the model is known beforehand or it can be approximated somehow, using either theory or measurements. If the model is not known, the method can still be used by using some flexible error model, e.g. polynomial.
Shifting of peaks causes problems in many applications, e.g. in NMR spectroscopy and chromatography.
Algorithm
Parameters
- Error model for independent variable (e.g. polynomial)
- Pure component signals
- Curve fitting method used in conjunction with x-axis correction (e.g. regular non-negative least squares, NNLS)
Input
- Mixture signal
- Initial guess for model parameters (optional)
Output
- Estimated pure component contributions
- Estimated error model parameters
Description
Make an initial guess for error model parameters
While True:
- Generate new set of independent variable values using the current estimate for error model parameters
- Interpolate measured signal to new x-axis formed by new variable values
- Calculate solution using regular curve fitting method and interpolated signal
- Calculate RMSE for solution
- Check termination condition and return solution if true
- Update error model parameters
The update of error model parameters can be done in different ways. In this project, three different algorithms were implemented.
Grid search
All candidates are tested and the best option is returned as the final solution.
Gauss-Newton
The parameter estimates are updated using the Gauss-Newton optimization method. This method uses gradient to find direction and step size for the update.
Evolutionary algorithm
Parameter candidates are updated using an evolutionary algorithm. This method uses a population of random parameter combinations, evaluates candidates in population, and then generates a new population, based on results of the previous population.
Combinations
These three different update algorithms can be used separately or they can be combined to produce better results. For example, one can use the grid search to make the first rough estimate for parameters and then use more accurate Gauss-Newton to optimize the solution from there. Using the grid search first will reduce the risk that gradient-based Gauss-Newton optimization ends up to the local minimum instead of the global minimum. The evolutionary algorithm is much slower than Gauss-Newton but can handle local minima better.
Demonstration with simulated data
Data generation
For testing and demonstration purposes, a synthetic data test set was generated. This sample set contained 100 samples with 3 components. Errors to the independent variable were generated using a quadratic polynomial model. To make the simulation more realistic, also amplitude noise was added to mixture samples. The pure components and generate mixtures are plotted below.
Analysis
For the analysis, the same error model (quadratic polynomial) was used as for the data generation. All the 100 samples were analyzed with and without correction using non-negative least squares as the fitting method. For x-axis correction, the parameter update was done using Gauss-Newton.
In the figure below, fitting of one sample is illustrated. At the first iteration, the algorithm does regular NNLS fit with the initial guess (in this case, this means no x-axis correction at all). The error is large, as indicated by the large residual between measured and estimated signal. After that, the algorithm makes error model parameter update, using Gauss-Newton, and proceeds to the next iteration. In the next iteration, the x-axis correction is made with updated parameters and the residual is considerably smaller. This process is repeated until RMSE cannot be decreased anymore. For this sample, the algorithm converged after five iterations.
The figure below shows analysis results for all the 100 generated samples. As can be seen, the errors without x-axis correction are large. The errors with x-axis correction are due to generated amplitude noise.
Conclusions
This blog post described the problem that is related to shifting peaks in curve fitting and signal separation. The post described one type of algorithm that can be used to solve problems where shifting of peaks (errors in independent variable) follow some model.
To use the algorithm, the user needs to provide an error model that is followed by errors of independent variable. The algorithm works by looking minimum RMSE by changing error model parameters and doing a regular curve fit after this. This process is repeated iteratively. The accepted solution is the one that produces minimum RMSE. Error model parameter updates can be made using different algorithms, e.g. Gauss-Newton, grid search, or evolutionary algorithm. In this algorithm, correction of the x-axis and the actual curve fitting are separated; this means that x-axis correction can be used in conjunction with any curve fitting method.
