speeding up a gaussian fitting....

3 ビュー (過去 30 日間)
Chris E.
Chris E. 2012 年 10 月 28 日
コメント済み: Emma 2016 年 4 月 8 日
This is a function that will output the full width half max of a Gaussian or SuperGaussian. This function works slowly. I was wondering if I can speed it up. This takes 0.3-0.5 seconds to run an individual Gaussian and 294 second to complete the full data set (or about 49 min). This is a little long for what I need to do. Any time cut from this would be amazing! I have a lot of Gaussian that are being looked through to extract the data in the end. PLEASE help out if you know of some solution. Thank you!
Here is the example code of what I have working and the output in the end:
function testing()
YProjY = [7.3800 19.7800 20.4200 20.8800 20.7400 21.2000 23.1000 21.8600 21.2600 20.1800 20.8000 21.4800 21.8200 21.6000 21.6000 20.5800 20.5200 20.7400 21.3800 21.6000 21.2000 21.4000 21.5400 21.8200 21.9600 22.5600 25.3000 25.9200 25.5000 23.4600 23.2800 23.2600 23.9800 24.3600 23.8600 23.3000 23.0200 24.2200 25.1600 24.1600 25.0800 23.9600 24.3400 23.1800 23.3400 24.3200 25.0000 25.2800 26.4400 26.5600 27.8400 27.0000 27.2200 27.4000 26.6200 26.6200 26.7000 27.6800 27.3600 26.9000 26.9400 26.2200 25.9000 26.4800 27.1600 26.2600 27.1200 28.3600 30.0600 29.1800 28.4800 28.7800 29.5400 29.4400 29.4600 29.9600 31.1600 30.0400 29.3600 28.7800 29.4800 30.3800 31.1000 31.3600 33.1800 34.8400 35.7200 34.7200 32.9600 33.1800 34.3400 35.4600 35.8400 35.0600 35.3600 34.5000 35.4400 34.7200 35.4800 35.4400 37.0600 37.1600 38.5600 39.3800 40.4200 38.4800 38.5400 38.2200 39.0000 38.7800 38.8600 39.0800 40.7600 40.4400 41.2600 39.6600 40.5800 42.2800 44.2400 45.5600 46.3600 44.8400 45.0200 45.0200 45.7000 45.1000 45.7400 45.7800 46.6400 47.0800 47.3400 46.8600 48.9400 48.8200 49.6200 50.1800 52.0000 52.6200 53.5000 51.9200 53.7000 54.5800 56.8400 58.0800 57.9200 57.0000 57.4200 59.6800 59.6600 61.5600 61.9600 63.8600 62.8600 61.2400 62.4200 63.6600 64.2800 64.5600 65.2400 67.6400 66.8400 67.8000 67.6400 68.3200 68.8200 68.4800 68.9200 69.6200 71.4400 72.2800 73.4600 75.1200 76.6800 76.9000 76.4200 77.6400 79.4400 80.3400 80.6200 80.9200 83.1600 84.7200 87.2000 88.3000 89.4600 88.9000 90.3400 91.5600 95.4200 96.7200 97.3200 97.5200 99.4800 101.1400 102.8200 104.4600 106.3400 107.0800 109.0800 110.9600 112.6800 114.1200 116.3400 118.0800 120.5200 120.7400 121.9200 124.4200 128.9200 131.7600 133.8400 134.6400 135.4600 137.1800 139.1600 141.8600 144.6000 147.5800 149.0200 150.6600 151.6800 153.2800 155.2800 157.9000 161.9800 162.7000 163.5400 163.2600 166.2600 168.1200 168.2800 169.4200 171.0200 172.3600 174.0600 175.1200 176.9200 177.6000 178.2200 177.3800 176.4800 176.2400 176.0600 174.8400 175.7000 175.1000 173.2800 170.7800 168.6600 170.3400 170.7600 171.2600 169.7800 168.4200 166.7200 165.3200 163.1400 161.9400 160.9400 157.4200 154.6600 152.8600 150.6800 149.3400 146.0800 143.6200 140.2200 136.1200 132.1600 129.8200 127.8000 125.3600 122.4600 119.8000 117.2400 116.2000 114.3000 111.4400 110.5800 108.0600 105.7800 100.3400 97.0200 95.1600 94.6800 94.5000 92.6400 91.0200 89.2400 87.6400 86.6000 84.1800 81.0400 79.8200 79.3600 78.2800 76.4200 76.2400 76.9200 76.4400 73.4800 72.3400 71.5000 70.7200 68.9000 66.9000 65.3800 65.2600 65.5200 64.9000 64.1600 62.3200 60.4400 58.4200 58.4600 59.2800 60.7800 60.7000 58.6400 59.6200 61.0600 61.7000 58.7000 55.8400 54.7200 53.4600 54.5600 55.1400 55.0200 52.9800 51.3800 49.7800 49.8600 50.9600 51.0200 52.2000 52.6600 51.6800 51.7200 51.1000 52.4400 52.3800 50.8800 50.9600 49.4800 48.1200 47.4200 48.0600 48.8200 47.8400 46.2200 44.9600 44.1800 45.2600 44.8200 44.5200 44.1000 45.1000 45.5600 44.3000 42.7000 42.5200 43.1000 44.4800 44.9200 43.4400 41.6800 41.5800 40.0200 39.0200 37.9600 38.3400 38.2400 37.4400 36.8800 37.4000 39.1200 39.1200 39.2200 38.1200 38.4400 37.8800 36.9000 35.8200 36.1400 36.5800 36.7000 36.2800 35.7600 34.1000 33.5800 32.9800 33.8800 33.7600 34.5800 34.0600 35.1800 34.7600 36.2400 35.4600 35.0800 33.0800 32.6600 33.8200 33.0000 33.0200 31.0800 29.6600 29.0600 29.4600 29.7800 30.1400 31.5800 31.4400 30.4000 29.3400 29.5800 28.9400 29.0800 28.7000 28.5600 30.5600 30.3600 30.0000 29.1600 29.2400 28.4600 28.3800 28.0200 27.3000 29.2800 29.6400 30.5000 30.9200 31.3000 30.6000 28.3200 28.2200 27.0600 26.4400 25.6600 25.4600 24.9600 23.8800 24.0400 24.8200 26.1000 28.1200 29.5800 30.1000 28.3200 26.0800 25.4800 24.5600 25.7400 27.1000 26.5000 25.5400 26.0800 27.8000 28.3600 29.1000 28.7400 29.2800 28.2800 28.8400 27.2800 26.4600 25.7600 25.6200 26.0200 25.3000 25.9000 26.3000 25.0000 24.7800 25.5000 27.4400 29.5600 27.8600 25.3400 23.5600 22.7800 23.6000 24.7000 26.3000 27.4400 27.0200 25.6000 24.4600 23.8800 23.9200 23.5800 22.5800 22.5200 21.6200 20.8000 19.5600 19.3800 19.5400 21.0600 22.5000 21.6600 21.5600 21.0600 22.1200 21.7800 22.1800 21.8400 22.5800 21.6000 22.3400 22.6400 22.8600 21.7800 20.5800 21.1200 20.0200 19.0400 19.2400 20.0400 20.1400 20.0800 20.5000 19.2800 18.5400 18.1200 20.3600 19.8600 18.8000 17.4400 18.4600 18.6800 18.7600 16.8400 18.1200 17.8400 18.9000 19.0200 19.3200 18.8000 19.1400 19.8800 20.8600 19.5600 19.3000 19.0600 19.0000 19.3000 19.2000 19.5000 18.6400 19.0800 19.3600 20.3400 20.2400 18.6000 17.4200 16.9400 17.0600 19.1600 21.0000 21.5400 18.2600 16.4200 16.9200 16.7000 17.7000 17.1800 17.2200 16.0800 16.8400 16.0600 15.5200 14.7000 3.5400];
XProjY = 1:576;
%---Y Projection starts here---
% a(1) = base
% a(2) = A; amplitude
% a(3) = x0; center
% a(4) = D; sigma_0
% a(5) = N;
a1Y = 100;
a2Y = 200;
a3Y = 280;
a4Y = 50;
a5Y = 1.5;
a0Y = [a1Y,a2Y,a3Y,a4Y,a5Y];
optsY = optimset('TolX',1e-4,'MaxFunEvals',10000,'MaxIter',10000,'Display','on'); [fitpara] = fminsearch(@SuperGaussianY,a0Y,optsY,XProjY,YProjY);
function chisq = SuperGaussianY(a,X,Y)
SupGau = a(1)+a(2)*exp( -0.5*(abs(X-a(3))/(a(4))).^a(5) );
csq = (Y - SupGau).^2;
chisq = sum(csq);
return
The output returns:
fitpara =
21.2322 168.4825 239.5092 31.9488 1.0939
  2 件のコメント
Jan
Jan 2012 年 10 月 30 日
Please learn how to format code in this forum. It is really easy, but increases the readability substantially. Thanks.
Emma
Emma 2016 年 4 月 8 日
Hi, I was actually looking for ideas for a code to do exactly this! I realise you posted this quite a while ago, but was just wondering if you could give me any pointers on what you've done? (Sorry, your formatting is a bit confusing and I'm also quite rusty with MATLAB at the moment, so more commenting of the code would help me).

サインインしてコメントする。

採用された回答

Eric
Eric 2012 年 10 月 30 日
編集済み: Eric 2012 年 10 月 30 日
I don't have time to work this in detail now, but here are my thoughts:
1. One way to speed the optimization would be to calculate the gradients of your error metric with respect to the fit parameters. Then you could use another optimization routine that uses gradient information (e.g., fminunc if you have the Optimization Toolbox).
2. You might look for ways to algorithmically determine the initial parameter values. If you can find better initial values the optimization will be faster. One way to do this may be to break the operation up into a couple steps. You may have better success first performing an optimization to a Gaussian and then fitting to a super-Gaussian. Alternatively, if you know your data sets are similar you might use an average of previous fit parameters as your initial guess.
Of these two approaches I would start with the first. The difficulty here will be in analytically calculating derivatives of your error metric. If you can do this, supplying gradient information to a good optimization routine should significantly speed convergence.
Good luck,
Eric
  2 件のコメント
Chris E.
Chris E. 2012 年 11 月 1 日
Thanks Eric! At least it gives me a direction to go to. I will try to see what I can do.
Eric
Eric 2012 年 11 月 1 日
If you don't have the Optimization Toolbox, you might check out the software at:
described in the report at
Good luck,
Eric

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeSolver Outputs and Iterative Display についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by