I’ve updated my automatic 2d Gaussian surface fitting function, available in Matlab Central, to add a function to fit Gabors to noisy data. As I’ve discussed previously, fitting a parametric surface to noisy data is pretty trivial whether it’s a Gabor, Gaussian, or otherwise — it’s a straightforward application of numerical optimization that can be done straightforwardly with lsqcurvefit. The annoying bit is choosing the initial parameters when the error surface is non-convex, especially when you have a bunch of data to fit.

The Gaussian surface fitting function does this by an exhaustive search through parameter space, which is made efficient via FFT-based convolution. The new Gabor fitting function uses the fact that the absolute value of the Fourier transform of a Gabor is a Gaussian. It gets a good guess of the orientation, spatial frequency and bandwidth of the Gabor using this fact, locates the approximate position of this Gabor, and then fills in the details with lsqcurvefit.

I’ve removed the function for Gibbs sampling the Gaussian fit; it didn’t work too well in practice with data that does not conform well to the model assumptions. I think this is partially because the posterior is often composed of islands of high-probability separated by a sea of low probability, and the slice sampler is not well adapted to this configuration.

Update: I’ve now added in a Metropolis-Hastings sampling to replace to original Gibbs method. Not sure how much more reliable with will be but this scheme is at least less messy to implement. Let me know.