How to properly fit the data to lorentzian curve as I am getting a few errors...

clear all
% close all
clc
tic
RMS=[];
R=[];
P_threshold=[];
L=5; %[m] Fiber length
V1=[1];
P0=4.7960;
F1=0;
count=0;
count1=0;
% for V=V1
%     count1=count1+1;
%     for fsine=0.50e9:150e6:1.5e9
%         count=count+1;
%    while 1
%         P0=P0+10*F1;
%[W] Input power
%iterations=131072; %number of iterations for the calculation of the dynamics 
%% general constants
n=1.45; %Index of refraction
eps0=8.854e-12; % [F/m] Vacuum permittivity
mu0 = 4*pi*1e-7;%[H/m] Vacuum permeability
c=2.9979e8; % [m/sec] Speed of light
Z0=sqrt(mu0/eps0);   %[Ohm] Vacuum impedance
dt = 6e-12; dz=dt*c/n;    %Spacial and Temporal step sizes.
% dz=0.05;dt=n/c*dz;
N=round(L/dz); % Fiber length discretization
z=0:dz:(dz*N-dz);
T=20*2*L*n/c;   %time taken for 10 round trips
Nt=round(T/dt);

%% material characteristics
lambdaL=1.064e-6;            % [m] pump wavelength
WL=2*pi*c/lambdaL;       %[rad/sec]  pump frequency
GAMMAb=2*pi/17.5e-9; %[1/t] phonons decay rate(zeringue)
A=7.85e-11;   %[m^2] fiber's effective area
Va=5960;    %[m/s] Velocity of sound
%Omega=1.01e11;%[rad/sec] acoustic frequency
Omega=2*n*Va*WL/c; %[rad/sec] (zeringue)  
WS = WL-Omega;
% gammaE=0.8967;  %electrostrictive coefficient
gammaE=1.95;  %electrostrictive coefficient


kT=300*1.3806504e-23;       %[Joule]        at room temprature. 
roh0=2201;       % [kg/m^3]  mean density for SiO2, bulk (http://www.virginiasemi.com/pdf/generalpropertiesSi62002.pdf) . 
Q=2*kT*roh0*GAMMAb/(Va^2*A);     %Noise std.
g0=(gammaE^2*WL^2)/n/c^3/Va/roh0/GAMMAb; %[m/W] (Jenkins)
sigma=(WL*gammaE)/2/n/roh0/c;    %(derivation)
q=2*WL*n/c;                       %derivation
I1_0=P0/A;  %[W/m^2] Pump intensity

LAMBDA=gammaE*Omega/c/Z0/2/Va^2;    %derivation
kappa=g0*GAMMAb*n/2/Z0/LAMBDA;      %derivation
%ten_RT_iter=20*N;%corresponds to 10 * roundtrip time
Nz=N;
%% Coefficients
G=g0*I1_0*L;                         %Gain factor

Gth=21; %Agrawal, pp.361

%Elt=repelem(sqrt(P0/A/2/n/c/eps0),Nz);
Elt=repelem(0,Nz);
Elt(1)=sqrt(P0*Z0/A/2/n);
%Est=repelem(sqrt(P0/A/2/n/c/eps0),Nz);
Est=repelem(0,Nz);
%%
t = (-Nt/2:1:Nt/2-1)*dt;
% nu=(-Nt/2/T/dt:1/Nt/dt:Nt/2/T/dt-1);
nu = (-Nt/2:1:Nt/2-1)*1/T;
vpi = 3;
 RL=50;
 V=1;
 fsine=2e9;
  sine = V*sin(2*pi*fsine*t);
 %% White noise modulation
 % randn('state', seed_phaseNoise)
 %  randn('state',0);
  % for i=0:2
 np = 22; % power of RF White noise signal in dBm
 np_l = 10^(np/10)*1e-3; % power of RF White noise signal in Watts
nbw = 0.5e9; % bandiwdth in Hz of the RF WNS

  rng('default');
  rng(2);
 n1 =randn(1, Nt); % generate a white gaussian noise in time domain which is delta correlated in time domain
  %  n_lpf1 = zeros(1, Nt);
  % n_FT = fftshift((fft(n1-mean(n1))));
    % SGF =exp(-0.5*(nu/(nbw)).^10); % Create a super-gaussian low pass filter with bandwidth = noise bandwidth
 % lor=(nbw^2)./(nu.^2+nbw^2);
 % x_axis=linspace(-1e9,1e9,1e4);
 % s=(sinc(x_axis));
s=nbw*sinc(pi*t*nbw);
  n_LPF_T=conv(n1,s,'same');
   % n_LPF = ((abs(n_FT).^2).*SGF); % Filter the White gaussian PSD with the PSD of LPF
   % n_LPF = ((abs(n_FT)).*SGF);
 % n_nor=n_LPF./max(n_LPF);
 % figure;plot(nu/1e9,n_nor);
 % xlim([-13 13]);
 % 
    % n_LPF_T = real(ifftshift(ifft(n_LPF-mean(n_LPF)))); % Filtered white gaussian noise in time domain
    % n_LPF_T = (ifftshift((ifft(n_LPF-mean(n_LPF))))); % Filtered white gaussian noise in time domain
  % n_LPF_T = n_LPF_T-mean(n_LPF_T);
 % nn = conv(n1-mean(n1),real(ifftshift(ifft(SGF-mean(SGF)))),'same');
 % n_LPF_T = nn;
    n_lpf = sqrt(np_l*T*RL)*(n_LPF_T).*(1./sqrt(trapz(t, abs(n_LPF_T).^2))); %normalize the noise volatge signal
   % n_lpf = sqrt(np_l*T)*(n_LPF_T).*(1./sqrt(trapz(t, abs(n_LPF_T).^2))); %normalize the noise volatge signal
% n_lpf = sqrt(np_l*T*RL)*(n_LPF).*(1./sqrt(trapz(t, abs(n_LPF).^2))); %normalize the noise volatge signal
%%

 phi=n_lpf;
 ES_0t=sqrt(I1_0/2/n/c/eps0)*exp((1i*phi));  % Original signal in time 
 % figure;plot(nu,fftshift(2*n*c*eps0*A*(abs(fft((ES_0t.^2))))));
 % xlim([-13e9 13e9]);

Power=trapz(t,2*n*c*eps0*A*abs(ES_0t).^2)/T;        % Area under the curve in time domain 

FFt_EL0t=fftshift(abs(fft(ES_0t)));             % fourier transform of the original signal 
Power_FFt=T*trapz(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2);   % Area under the curve in frequency domain 
 figure;plot(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2)
 title('dt=',num2str(dt));

 % Assuming x and y are defined as above
 [xData, yData] = deal(nu, 2*n*c*eps0*A*(FFt_EL0t/Nt).^2);

[yprime, params, resnorm, residual, jacobian] = lorentzfit(xData, yData);

Unrecognized function or variable 'lorentzfit'.

figure;plot(nu,yprime);

How to properly fit the curve the lorentzian curve to the data...

[Sam Chak]

@Yogesh,

Rather than using the special function, lorentzfit(), right from the beginning, please describe your ideal or desired mathamtical model of the lorentzian curve. Can you do so?

[Yogesh]

Hii , thank you for answering but even if I define my ideal model of Lorentzian , How do I find FWHM or width of the curve with just data???...

[Sam Chak]

How does your Lorentzian model look like? Your data is very noisy.

[Yogesh]

Hey I tried using the standard equation for lorentzian a/((x-b)^2+a^2) but it's not a good fit... I'm confused at the moment. ... Should I try smoothening the curve for better results ..

Answer 1

[Sam Chak]

Hi @Yogesh

Here is how I fit the Lorentzian model:

RMS=[];
R=[];
P_threshold=[];
L=5; %[m] Fiber length
V1=[1];
P0=4.7960;
F1=0;
count=0;
count1=0;

%% general constants
n=1.45; %Index of refraction
eps0=8.854e-12; % [F/m] Vacuum permittivity
mu0 = 4*pi*1e-7;%[H/m] Vacuum permeability
c=2.9979e8; % [m/sec] Speed of light
Z0=sqrt(mu0/eps0);   %[Ohm] Vacuum impedance
dt = 6e-12; dz=dt*c/n;    %Spacial and Temporal step sizes.
N=round(L/dz); % Fiber length discretization
z=0:dz:(dz*N-dz);
T=20*2*L*n/c;   %time taken for 10 round trips
Nt=round(T/dt);

%% material characteristics
lambdaL=1.064e-6;            % [m] pump wavelength
WL=2*pi*c/lambdaL;       %[rad/sec]  pump frequency
GAMMAb=2*pi/17.5e-9; %[1/t] phonons decay rate(zeringue)
A=7.85e-11;   %[m^2] fiber's effective area
Va=5960;    %[m/s] Velocity of sound
Omega=2*n*Va*WL/c; %[rad/sec] (zeringue)  
WS = WL-Omega;
gammaE=1.95;  %electrostrictive coefficient
kT=300*1.3806504e-23;       %[Joule]        at room temprature. 
roh0=2201;       % [kg/m^3]  mean density for SiO2, bulk (http://www.virginiasemi.com/pdf/generalpropertiesSi62002.pdf) . 
Q=2*kT*roh0*GAMMAb/(Va^2*A);     %Noise std.
g0=(gammaE^2*WL^2)/n/c^3/Va/roh0/GAMMAb; %[m/W] (Jenkins)
sigma=(WL*gammaE)/2/n/roh0/c;    %(derivation)
q=2*WL*n/c;                       %derivation
I1_0=P0/A;  %[W/m^2] Pump intensity
LAMBDA=gammaE*Omega/c/Z0/2/Va^2;    %derivation
kappa=g0*GAMMAb*n/2/Z0/LAMBDA;      %derivation
Nz=N;

%% Coefficients
G=g0*I1_0*L;                         %Gain factor
Gth=21; %Agrawal, pp.361
Elt=repelem(0,Nz);
Elt(1)=sqrt(P0*Z0/A/2/n);
Est=repelem(0,Nz);

%%
t = (-Nt/2:1:Nt/2-1)*dt;
nu = (-Nt/2:1:Nt/2-1)*1/T;
vpi = 3;
RL=50;
V=1;
fsine=2e9;
sine = V*sin(2*pi*fsine*t);

%% White noise modulation
np = 22; % power of RF White noise signal in dBm
np_l = 10^(np/10)*1e-3; % power of RF White noise signal in Watts
nbw = 0.5e9; % bandiwdth in Hz of the RF WNS

rng('default');
rng(2);
n1 = randn(1, Nt); % generate a white gaussian noise in time domain which is delta correlated in time domain
s=nbw*sinc(pi*t*nbw);
n_LPF_T=conv(n1,s,'same');
n_lpf = sqrt(np_l*T*RL)*(n_LPF_T).*(1./sqrt(trapz(t, abs(n_LPF_T).^2))); %normalize the noise volatge signal

%%
phi=n_lpf;
ES_0t=sqrt(I1_0/2/n/c/eps0)*exp((1i*phi));  % Original signal in time 
Power=trapz(t,2*n*c*eps0*A*abs(ES_0t).^2)/T;        % Area under the curve in time domain 
FFt_EL0t=fftshift(abs(fft(ES_0t)));             % fourier transform of the original signal 
Power_FFt=T*trapz(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2);   % Area under the curve in frequency domain 
% figure;plot(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2)

%% ---------------------------------
%% Fitting with Lorentzian model
%% ---------------------------------

y   = 2*n*c*eps0*A*(FFt_EL0t/Nt).^2;
[up, ~] = envelope(y, round(0.01*size(nu/max(nu), 2)), 'peak');

format long g
x   = nu;
xm  = max(nu)

xm = 
          83332316172.4138

% Fitting model
fo  = fitoptions('Method', 'NonlinearLeastSquares',...
               'Lower', [0.00008, 0.01, -0.002],...
               'Upper', [0.00032, 0.02,  0    ]);
ft  = fittype('a*b/((x/83332316172.4138 - c)^2 + b^2)', ... 
    'independent', {'x'}, ...
    'coefficients', {'a', 'b', 'c'}, ...
    'options', fo);

%% Fit curve to data
[yfit, gof] = fit(x', up', ft, 'StartPoint', [0.00015, 0.0075, 0.001752])

yfit = 
     General model:
     yfit(x) = a*b/((x/83332316172.4138 - c)^2 + b^2)
     Coefficients (with 95% confidence bounds):
       a =    0.000166  (0.0001658, 0.0001662)
       b =     0.01318  (0.01315, 0.0132)
       c =   -0.001571  (-0.001586, -0.001555)

gof = struct with fields:
           sse: 0.00746975520339178
       rsquare: 0.971847764622057
           dfe: 161221
    adjrsquare: 0.971847415384236
          rmse: 0.000215249613065384

figure
plot(yfit, x, y), grid on, grid minor, xlim([-2e10, 2e10])
xlabel('\nu'), ylabel('y')

[Yogesh]

Hey thanks for the answer but i had one doubt , why are you tryingh to fit the curve to envelope ...

Shouldn't we try to fit it ti the data 'y' instead of 'up'???

[Torsten]

@Yogesh

Your data form a point cloud, not a curve. Thus usual curve fitting is inappropriate,

[Sam Chak]

Hi @Yogesh

Because the scattered data exhibits a Gaussian distribution (up to R² = 0.9995), your specific request for a Lorentzian distribution (R² = 0.9718) is used to fit the distribution envelope of the data.

In fact, the data are distributed irregularly across the domain. For example, how would you fit this mass of data points within the red-circled area? There is no identifiable curve, other than the distribution envelope.

If this is not what you are looking for, please share your comments and provide a sketch here so that other data scientists 👨🏻‍💻 in the forum can assist you. 👍

[Yogesh]

Oh so if I use gaussian equation then it would be a better fit...

Answer 2

[M.]

Are you referring to this lorentzfit function ? If so, make sure it is in your matlab path.

You can also have a look at these others functions :

https://mathworks.com/matlabcentral/fileexchange/63748-fitnlorentzian-xydata-parameters

https://mathworks.com/matlabcentral/fileexchange/13648-lorentzian-fit

Answer 3

[Sam Chak]

Hi @Yogesh

Here is the fitting with a Gaussian model.

RMS=[];
R=[];
P_threshold=[];
L=5; %[m] Fiber length
V1=[1];
P0=4.7960;
F1=0;
count=0;
count1=0;

%% general constants
n=1.45; %Index of refraction
eps0=8.854e-12; % [F/m] Vacuum permittivity
mu0 = 4*pi*1e-7;%[H/m] Vacuum permeability
c=2.9979e8; % [m/sec] Speed of light
Z0=sqrt(mu0/eps0);   %[Ohm] Vacuum impedance
dt = 6e-12; dz=dt*c/n;    %Spacial and Temporal step sizes.
N=round(L/dz); % Fiber length discretization
z=0:dz:(dz*N-dz);
T=20*2*L*n/c;   %time taken for 10 round trips
Nt=round(T/dt);

%% material characteristics
lambdaL=1.064e-6;            % [m] pump wavelength
WL=2*pi*c/lambdaL;       %[rad/sec]  pump frequency
GAMMAb=2*pi/17.5e-9; %[1/t] phonons decay rate(zeringue)
A=7.85e-11;   %[m^2] fiber's effective area
Va=5960;    %[m/s] Velocity of sound
Omega=2*n*Va*WL/c; %[rad/sec] (zeringue)  
WS = WL-Omega;
gammaE=1.95;  %electrostrictive coefficient
kT=300*1.3806504e-23;       %[Joule]        at room temprature. 
roh0=2201;       % [kg/m^3]  mean density for SiO2, bulk (http://www.virginiasemi.com/pdf/generalpropertiesSi62002.pdf) . 
Q=2*kT*roh0*GAMMAb/(Va^2*A);     %Noise std.
g0=(gammaE^2*WL^2)/n/c^3/Va/roh0/GAMMAb; %[m/W] (Jenkins)
sigma=(WL*gammaE)/2/n/roh0/c;    %(derivation)
q=2*WL*n/c;                       %derivation
I1_0=P0/A;  %[W/m^2] Pump intensity
LAMBDA=gammaE*Omega/c/Z0/2/Va^2;    %derivation
kappa=g0*GAMMAb*n/2/Z0/LAMBDA;      %derivation
Nz=N;

%% Coefficients
G=g0*I1_0*L;                         %Gain factor
Gth=21; %Agrawal, pp.361
Elt=repelem(0,Nz);
Elt(1)=sqrt(P0*Z0/A/2/n);
Est=repelem(0,Nz);

%%
t = (-Nt/2:1:Nt/2-1)*dt;
nu = (-Nt/2:1:Nt/2-1)*1/T;
vpi = 3;
RL=50;
V=1;
fsine=2e9;
sine = V*sin(2*pi*fsine*t);

%% White noise modulation
np = 22; % power of RF White noise signal in dBm
np_l = 10^(np/10)*1e-3; % power of RF White noise signal in Watts
nbw = 0.5e9; % bandiwdth in Hz of the RF WNS

rng('default');
rng(2);
n1 = randn(1, Nt); % generate a white gaussian noise in time domain which is delta correlated in time domain
s=nbw*sinc(pi*t*nbw);
n_LPF_T=conv(n1,s,'same');
n_lpf = sqrt(np_l*T*RL)*(n_LPF_T).*(1./sqrt(trapz(t, abs(n_LPF_T).^2))); %normalize the noise volatge signal

%%
phi=n_lpf;
ES_0t=sqrt(I1_0/2/n/c/eps0)*exp((1i*phi));  % Original signal in time 
Power=trapz(t,2*n*c*eps0*A*abs(ES_0t).^2)/T;        % Area under the curve in time domain 
FFt_EL0t=fftshift(abs(fft(ES_0t)));             % fourier transform of the original signal 
Power_FFt=T*trapz(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2);   % Area under the curve in frequency domain 
% figure;plot(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2)

%% ---------------------------------
%% Fitting with Gaussian model
%% ---------------------------------

y   = 2*n*c*eps0*A*(FFt_EL0t/Nt).^2;
[up, ~] = envelope(y, round(0.01*size(nu/max(nu), 2)), 'peak');

format long g
x   = nu;
xm  = max(nu)

xm = 
          83332316172.4138

% Fitting model 
fo  = fitoptions('Method', 'NonlinearLeastSquares',...
               'Lower', [0.01, 0.01, -0.002],...
               'Upper', [0.02, 0.02,  0    ]);
ft  = fittype('a*exp(- 1/2*((x/83332316172.4138 - c)/b)^2)', ... 
    'independent', {'x'}, ...
    'coefficients', {'a', 'b', 'c'}, ...
    'options', fo);

%% Fit curve to data
[yfit, gof] = fit(x', up', ft, 'StartPoint', [0.015, 0.015, -0.001])

yfit = 
     General model:
     yfit(x) = a*exp(- 1/2*((x/83332316172.4138 - c)/b)^2)
     Coefficients (with 95% confidence bounds):
       a =     0.01156  (0.01155, 0.01156)
       b =     0.01426  (0.01426, 0.01426)
       c =   -0.001754  (-0.001756, -0.001751)

gof = struct with fields:
           sse: 0.00020534180419802
       rsquare: 0.999226101706507
           dfe: 161221
    adjrsquare: 0.999226092106042
          rmse: 3.56884660350626e-05

figure
plot(yfit, x, y), grid on, grid minor, xlim([-1e10, 1e10])
title('Fitting Point Cloud distribution with Gaussian model')
xlabel('\nu'), ylabel('y')