M^2 fitting for a Gaussian beam with measured data

I have modified the following code originally written by #Ka Shun Wu# (for copyrigth purposes) to try and fit my measured spatial beam profile data. The code runs into a fitting and insufficient data points errors! ( I have attached a sample data .csv files for beam data at different positions along the propagation direction and I have 31 similar files. How should I resolve this issue?

clear all
close all
%profile_loc=30:10:160;%distance of @profile from arb point (mm)
profile_loc=550:10:850;%distance of @profile from arb point (mm)
wavelenght=9.2e-3;%wavelenght in m
pixel_width=80;%cameras pixel width in um
pixel_height=80;% camears pixel width in um
blur_factor=5;%pixel averaging_to estimate beam center
back_av=2;
suppress_plot=1;%1=only plot beam propagation and final M^2
do_M2_fit=0;%1=do_M^2 fit; 0=just do Gaussian fit to profiles

%Define the directory where the files are located
%directory = 'C:\examples\Tempertaure gradient matlab files\Tempertaure gradient matlab files';

%Get a list of all CSV files in the directory
%file_list = dir(fullfile(directory, '*.csv'));

%Loop through each file and rename it
%for i = 1:numel(file_list)
%old_name = fullfile(directory, file_list(i).name);
%new_name = fullfile(directory, ['profile_10', num2str(i), '.csv']);
% movefile(old_name, new_name);
%end

%filename = fullfile(pwd, 'profile_.csv');

%file name for exporting of beam width and height data
output_file='beam_width_data.xls';

num_loop=length(profile_loc);%min 3 for an M^2 fit
x_width=zeros(num_loop, 1);
y_width=zeros(num_loop, 1);
distance=zeros(num_loop, 1);
num_fig=ceil(num_loop/6);%number of figures required

for a=1:num_loop
    dist=profile_loc(a);
    %filename=['56_cm', num2str(dist) '.csv'];
    filename = ['profile_', num2str(dist), '.csv'];
    
    %filename=['spreadsheet.csv',1,1];
    RAW=csvread(filename);
    RAW=transpose(RAW);%depends upon camera orientation
    
    %remove background from 4 samples in blurred image by taking the
    %average of 4 squares in the corners
    ii=size(RAW);
    back=mean(mean([RAW(1:back_av,1:back_av)...
        RAW(ii(1)-back_av+1:ii(1),1:back_av)...
        RAW(1:back_av,ii(2)-back_av+1:ii(2))...
        RAW(ii(1)-back_av+1:ii(1),ii(2)-back_av+1:ii(2))]));
    RAW=RAW-back;
    %blurring image and finding max of profile
    %adjust blur_factor above if necessary. Large number=more blur
    window=(1:blur_factor)*0+1/blur_factor;
    blurred=convn(RAW,window,'same');
    blurred=transpose(convn(transpose(blurred),window,'same'));  
    [colMax max_loc_y]=max(blurred);
    [peak max_loc_x]=max(colMax);
    max_loc_y=max_loc_y(max_loc_x);
    %taking the row and column data centered at the beam centre
    x_mode=transpose(RAW(max_loc_y,:));
    width=transpose(1:length(x_mode))*pixel_width;
    y_mode=RAW(:,max_loc_x);
    height=transpose(1:length(y_mode))*pixel_height;
    %Gaussian curve fitting and extracting the beam width and heigth
    f=fittype('gauss1');
    x_fit=fit(width,x_mode,f);
    y_fit=fit(height,y_mode,f);
    x_coeff=coeffvalues(x_fit);
    y_coeff=coeffvalues(y_fit);
    distance(a)=dist;
    x_width(a)=x_coeff(3)*sqrt(2);
    y_width(a)=y_coeff(3)*sqrt(2);
    %plotting each image with the center line and centroid in an array
    if (suppress_plot==0)
        rows=(num_loop>2)+1;
        cols=rows+(num_loop==2)+(num_loop>4);
        figure(ceil(a/6))
        subplot(rows,cols,(a-6*floor(a-1/6)))
        image(RAW,'CDataMapping','scaled')
        line([max_loc_x,max_loc_x],[0,length(RAW)],'color','w')
        line([0,length(RAW)],[max_loc_y,max_loc_y],'color','w')
        rectangle('position',[max_loc_x-(x_width(a)/pixel_width),...
            max_loc_y(y_width(a)/pixel_heigth),2*(x_width(a)/pixel_width),...
            2*(y_width(a)/pixel_heigth)],'Curvature',[1,1],'EdgeColor','w')
        title(['Profile at' num2str(dist) 'mm'],'FontWeigth','bold')
        set(gca,'DataAspectRatio',[1 1 1],'xtick',[],'ytick',[])
        %plotting each fit (x-axis and y-axis) in an array
        figure(ceil(a/6)+num_fig)
        subplot(rows,cols*2,(a*2-1)-12*floor((a*2-2)/12))
        plot(width,x_mode)
        hold on
        plot(x_fit,'m')
        legend off
        hold off
        title(['x-axis fit' num2str(dist) 'mm'],'FontWeight','bold')
        set(gca, 'xtick',[], 'ytick',[])
        axis ([min(width) max(width)...
            min(x_mode)-abs(min(x_mode)*1.2) max(x_mode)*1.2])
        subplot(rows,cols*2,(a*2)-12*floor((a*2-2)/12))
        plot(height, y_mode)
        hold on
        plot (y_fit, 'm')
        legend off
        hold off
        title(['y-axis fit' num2str(dist) 'mm'],'FontWeight','bold')
        set(gca, 'xtick',[], 'ytick',[])
        axis ([min(width) max(width)...
            min(y_mode)-abs(min(y_mode)*1.2) max(y_mode)*1.2])
    end
end

Error using fit>iFit (line 273)
Insufficient data.
You need at least 3 data points to fit this model.

Error in fit (line 117)
[fitobj, goodness, output, convmsg] = iFit( xdatain, ydatain, fittypeobj, ...

%call M2_fitting to do the M2 fit
if (num_loop>=3)&&(do_M2_fit==1)
    [x_M2,x_M2_fit]=M2_fitting(distance,x_width/1000,wavelength);
    [y_M2,y_M2_fit]=M2_fitting(distance,y_width/1000,wavelength);
end
%write beam width and height data to file
fileID=fopen(output_filename,'w');
fprint(fileID, 'location[mm]\tbeam height [mm]\tbeam width [mm]\n');
for a=1:num_loop
    fprintf(fileID,[num2str(distance(a)) '\t' num2str(y_width(a)/1000)...
        '\t' num2str(x_width(a)/1000) '\n']);
end
if (num_loop>=3)&&(do_M2_fit==1)
    fprintf(fileID,['\nEstimated M^2:\n\t' num2str(y_M2) '\t'...
        num2str(x_M2)]);
end
fclose(fileID);
%plotting profiles of the beam propagation
if (num_loop>=3)
    b=str2double(num2str(max(x_width),0))*2;
    c=str2double(num2str(max(y_width),0))*2;
    xw_axis=linspace(-b,b,100);
    xh_axis=linspace(-c,c,100);
    func=zeros(num_loop,100);
    func2=zeros(num_loop,100);
    peakx=1./(x_width*sqrt(pi()));
    peaky=1./(y_width*sqrt(pi()));
    [X1,Y1]=meshgrid(distance,xw_axis/1000);
    [X2,Y2]=meshgrid(distance,xh_axis/1000);
    for i=1:num_loop
        func(i,:)=peakx(i)*exp(-(xw_axis./x_width(i)).^2);
        func2(i,:)=peaky(i)*exp(-(xh_axis./y_width(i)).^2);
    end
    func=transpose (func);
    func2=transpose (func2);
    figure
    subplot(1,2,1)
    plot3(Y1,X1,func)
    hold on
    contour (Y1,X1,func)
    title('horizontal Gaussian fits', 'FontWeight', 'bold')
    xlabel('Width [mm]')
    ylabel('Beam propagation [mm]')
    zlabel('Intensity [A.U]')
    grid on
    subplot (1,2,2)
    plot3(Y2,X2,func2)
    hold on
    contour (Y2,X2,func2)
    title('vertical Gaussian fits', 'FontWeight', 'bold')
    xlabel('Height [mm]')
    ylabel('Beam propagation [mm]')
    zlabel('Intensity [A.U]')
    grid on
end
%plotting the M^2 fit to beam widths
if (num_loop>=3)&&(do_M2_fit==1)
    figure
    subplot(1,2,1)
    hold on
    plot(distance,x_width/1000,'o')
    plot(x_M2_fit)
    title(['horizontal M^2 fit: M^2= ' num2str(x_M2)], 'FontWeight', 'bold')
    xlabel('Beam propagation [mm]')
    ylabel('Beam width [mm]')
    subplot (1,2,2)
    hold on
    
    plot(distance,y_width/1000,'o')
    plot(y_M2_fit)
    title(['vertical M^2 fit: M^2= ' num2str(y_M2)], 'FontWeight', 'bold')
    xlabel('Beam propagation [mm]')
    ylabel('Beam height [mm]')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

[Matt J]

The code runs into a fitting and insufficient data points errors!

Show us. We can't see it.

[Dismas hoge]

@Matt J

>> M_squared_fitting

Error using fit>iFit (line 261)

Insufficient data.

You need at least 3 data points to fit this model.

Error in fit (line 108)

[fitobj, goodness, output, convmsg] = iFit( xdatain, ydatain, fittypeobj, ...

Error in M_squared_fitting (line 70)

y_fit=fit(height,y_mode,f);

[Catalytic]

You need at least 3 data points to fit this model.

The error message makes it look obvious. Have you even checked whether you've given at least 3 data points to fit?

[Dismas hoge]

@Catalytic all I have is several .csv files having series of x and y values. I expect that the data points are enough!

[Catalytic]

But Matlab is telling you (or hinting heavily) that your expectations are not fulfilled. You don't seem to have verified how many data points are present in the data set that is raising the error.

Answer 1

[Matt J]

Your post doesn't give much detail on what you are trying to do, but if you are just trying to fit a 1D Gaussian to the data in your .csv files, you can download gaussfitn from,

https://www.mathworks.com/matlabcentral/fileexchange/69116-gaussfitn

Example:

[x,y]=readvars('profile_550.csv');

mu0=x(1);
sig0=x(find(y<=y(1)/2,1))*2/2.355; %initial gueses of mean and standard dev

p=gaussfitn(x,y,{min(y), max(y), mu0,sig0^2});

[D,A,mu,var]=deal(p{:}); %estimated parameters

f=@(x)D+A*exp(-0.5*(x-mu).^2/var);

plot(x,y,'x'); xlim([0,50]); hold on; 
fplot(f,xlim); hold off; 
legend Data 'Gaussian Fit'

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

[Dismas hoge]

Thanks @Matt J, actually I have no problem with Gausssian fit to individual .csv data but rather a plot of radii at 1/e^2 from each each profile plot vs the propagation positions. I expect a caustic curve!

[Matt J]

Based on the error message you've shown us, you do have a problem with the Gaussian fitting. The error message is coming from fit() and the only places in your code where fit() is used is when doing a Gaussian fit.

[Dismas hoge]

Gaussian form for such a beam takes the form I=I0exp(-2r^2/w^2) which for this case, r are the x values and I are the y values so that w should be the r value at 1/e^2 of the max I.

[Matt J]

So your questions are all resolved? I assume so since you Accept-clicked my answer.

[Dismas hoge]

Not yet resolved!

[Matt J]

Well the error message you posted should be resolved based on @Catalytic's comments and if that error message is no longer the issue, it is not clear what is.

[Dismas hoge]

I still can't figure out why the error. I have enough data points form the file!

[Catalytic]

I still can't figure out why the error. I have enough data points form the file!

It is clear that you have done nothing to check whether the variables height and y_mode contain more than two points, even though the error message told you they don't. We can run the code right here online and easily show that they contain only two points -

clear all
close all
%profile_loc=30:10:160;%distance of @profile from arb point (mm)
profile_loc=550:10:850;%distance of @profile from arb point (mm)
wavelenght=9.2e-3;%wavelenght in m
pixel_width=80;%cameras pixel width in um
pixel_height=80;% camears pixel width in um
blur_factor=5;%pixel averaging_to estimate beam center
back_av=2;
suppress_plot=1;%1=only plot beam propagation and final M^2
do_M2_fit=0;%1=do_M^2 fit; 0=just do Gaussian fit to profiles

%Define the directory where the files are located
%directory = 'C:\examples\Tempertaure gradient matlab files\Tempertaure gradient matlab files';

%Get a list of all CSV files in the directory
%file_list = dir(fullfile(directory, '*.csv'));

%Loop through each file and rename it
%for i = 1:numel(file_list)
%old_name = fullfile(directory, file_list(i).name);
%new_name = fullfile(directory, ['profile_10', num2str(i), '.csv']);
% movefile(old_name, new_name);
%end

%filename = fullfile(pwd, 'profile_.csv');

%file name for exporting of beam width and height data
output_file='beam_width_data.xls';

num_loop=length(profile_loc);%min 3 for an M^2 fit
x_width=zeros(num_loop, 1);
y_width=zeros(num_loop, 1);
distance=zeros(num_loop, 1);
num_fig=ceil(num_loop/6);%number of figures required

for a=1:num_loop
    dist=profile_loc(a);
    %filename=['56_cm', num2str(dist) '.csv'];
    filename = ['profile_', num2str(dist), '.csv'];
    
    %filename=['spreadsheet.csv',1,1];
    RAW=csvread(filename);
    RAW=transpose(RAW);%depends upon camera orientation
    
    %remove background from 4 samples in blurred image by taking the
    %average of 4 squares in the corners
    ii=size(RAW);
    back=mean(mean([RAW(1:back_av,1:back_av)...
        RAW(ii(1)-back_av+1:ii(1),1:back_av)...
        RAW(1:back_av,ii(2)-back_av+1:ii(2))...
        RAW(ii(1)-back_av+1:ii(1),ii(2)-back_av+1:ii(2))]));
    RAW=RAW-back;
    %blurring image and finding max of profile
    %adjust blur_factor above if necessary. Large number=more blur
    window=(1:blur_factor)*0+1/blur_factor;
    blurred=convn(RAW,window,'same');
    blurred=transpose(convn(transpose(blurred),window,'same'));  
    [colMax max_loc_y]=max(blurred);
    [peak max_loc_x]=max(colMax);
    max_loc_y=max_loc_y(max_loc_x);
    %taking the row and column data centered at the beam centre
    x_mode=transpose(RAW(max_loc_y,:));
    width=transpose(1:length(x_mode))*pixel_width;
    y_mode=RAW(:,max_loc_x);
    height=transpose(1:length(y_mode))*pixel_height;
    %Gaussian curve fitting and extracting the beam width and heigth
    f=fittype('gauss1');
    x_fit=fit(width,x_mode,f);


    whos height y_mode ; break

end

  Name        Size            Bytes  Class     Attributes

  height      2x1                16  double              
  y_mode      2x1                16  double              

[Matt J]

@Catalytic is right. However, it would have been easier just to run the code with "Pause On Errors" selected,

and then examine height and y_mode in the Editor when the code stops,

[Dismas hoge]

So, if I make a change to back_av value to 3, I get a matrix deminsions error!

clear all

close all

%projection_loc=550:10:600;%distance of @profile from arb point (mm)

profile_loc=550:10:850;%distance of @profile from arb point (mm)

wavelenght=9.2e-3;%wavelenght in m

pixel_width=80;%cameras pixel width in um

pixel_height=80;% camears pixel width in um

blur_factor=5;%pixel averaging_to estimate beam center

back_av=3;

Index exceeds matrix dimensions.

Error in M_squared_fitting (line 52)

back=mean(mean([RAW(1:back_av,1:back_av)...

This is limiting the y_mode and height!