Curve fitting

Shared Code
#1

I am doing some experiments with optics and ray tracing. This involves calculating the amplitude at various points across the image plane, given a known frequency and an unknown phase and amplitude. Currently I just generate a dozen points on the curve and pick the one with the highest amplitude, but not only is that lazy and inaccurate, it is also very expensive to generate a dozen samples per pixel of the image plane (it involves iterating over hundreds of thousands of samples of the source wavefront).

it would be better to generate two or three samples and then use some kind of curve fitting to get the amplitude. There is probably an algebraic method, but that would involve maths.

So here is my first stab at a curve fitting algorithm. It involves 200 iterations to get a curve match with an error of ~0.01%, although the method occasionally glitches out with an error of 1%.

Original: http://pastebin.java-gaming.org/f023f14592214

Edit: Latest: http://pastebin.java-gaming.org/23f4932542610

When I run it I get output like this:

Error by curve fitting : 0.00263
Error by rough curve sampling : 1.17135
Error by fine curve sampling : 0.28717

#2

Sounds a bit to me like you want to emulate a fourier transform? But I’m not sure if I understood you right. If you want to transform an image into frequencies, you can apply a 2D fourier transform, but your application seems to mean to transform a line of pixels, which is a 1D fourier transform.

#3

Yes, I am exactly emulating a fourier transform. I struggle to understand the FT maths but I do understand the physics, so I am building a simulation instead.

#4

Maybe this page can help you - it has some pseudocode:

http://www.nayuki.io/page/how-to-implement-the-discrete-fourier-transform

Similar, with C code:

http://www.scratchapixel.com/old/lessons/mathematics-physics/discrete-fourier-transform-dft-part-1/1d-forward-and-inverse-dft-in-c/

#5

Thanks, I will have a look at these. I would like to understand FT, but I would need to understand it to a level where I knew the equations represent a physical reality. With ray tracing and simple addition of wavefronts, I know what I am doing because I avoid the mathematical sophistication. Whereas with FT, I can google some code and use it, but i can’t check my own work yet.

For example, with my program I can compare the Airy disc of my current telescope to a larger one I am thinking of buying:

With the curve fitting however, all I am trying to do is get back to a sine curve from a minimum number of samples with least effort. At the moment, I am just getting amplitude by brute force by calculating lots of samples.

#6

The idea of FT is not so different from yours, just that the basic functions are sin and cos waves of raising frequency. FT calculates coefficients for those waves, in order that the sum of those waves matches the original function.

You can simplify (or compress) the results if you leave out higher frequency coefficients, and remore coefficients which are near zero (e.g. JPEG compression does this.)

#7

I find the easiest way to understand the DFT is as least squares: http://en.wikipedia.org/wiki/Linear_least_squares_(mathematics).

You’re finding the best coefficients for a set of sinusoids so their sum is closest to your data. By using sufficient sinusoids of different frequencies these form a basis for all 1D signals, so you can get a perfect answer.

The X in that page are all the sinusoids as a matrix (one per column).
beta are the coefficients – they’re complex so you get magnitude and phase.
y is the data (your 1D signal).

#8

And once you have DFT working there is the FFT to make it fast.

#9

I’m not entirely sure that Fourier analysis is what is needed here as the best way to do the particular curve fitting you are trying to accomplish. But I did spend some time last year working through the math of the DFT, and collecting links. I managed to get the basic concepts of DFT and for programming convolution. I didn’t make it as far as FFT.

Maybe something here will be useful.

I think this is the main book I worked through. “The Scientist and Engineer’s Guide to
Digital Signal Processing”
http://www.dspguide.com/ch8/1.htm

ccrma is a fantastic resource, with sections for some of the background math.
https://ccrma.stanford.edu/~jos/mdft/DFT_Math_Outline.html

I have a link to some Java implementations of FFT, from StackOverflow, that I haven’t followed up on yet.

StackOverlow now has a DSP site:


I haven’t tried using it yet. Just signed up prior to this post.

#10

I’m really not keen on learning a bunch of fancy maths just to fit a curve. Fourier transforms are more interesting for the more general problem of transforming a telescope aperture to its Fraunhofer diffraction pattern, but ray tracing works fine for now.

Regarding the curve fitting problem that leaves me with… Firstly my curve fitting algo is not a general solution, it glitches severely for samples with a y value close to zero (amplitude errors up to 95%). I solved that problem by taking four samples and picking the pair of consecutive values furthest from zero. Changing from 12 samples per pixel to 4 samples per pixel obviously helps memory and reduces the size of my innermost loop by a factor of 3.

Secondly, the process needs to be quick enough for processing an image plane of 300x300 pixels. I find I can get a max error of ~0.04% in two seconds with just over 200 iterations per pixel. Two seconds is fast enough, given that the ray tracing itself is going to take between five minutes and five hours depending on the aperture being calculated.

I updated the pastebin with the small changes I’ve made to solve the glitchouts.


Avg Error by curve fitting        : 0.00116%
Avg Error by rough curve sampling : 1.13904%
Avg Error by fine curve sampling  : 0.28740%

Max Error by curve fitting        : 0.03782%
Max Error by rough curve sampling : 3.40740%
Max Error by fine curve sampling  : 0.85551%

Steps for curve fitting           : 221
Millis taken for ccd              : 1994.1
Millis taken for program          : 2166.0