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

31 次查看（过去 30 天）
Yogesh
Yogesh 2024-10-10
回答： Sam Chak 2024-10-16
Ran in:
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...
  4 个评论
显示 2更早的评论
Yogesh
Yogesh 2024-10-13
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 ..

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Sam Chak
Sam Chak 2024-10-13
Ran in:
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')
  4 个评论
显示 2更早的评论
Sam Chak
Sam Chak 2024-10-16
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. 👍

请先登录，再进行评论。

更多回答（2 个）

M.
M. 2024-10-10
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
Sam Chak
Sam Chak 2024-10-16
Ran in:
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')

类别

在 Help Center 和 File Exchange 中查找有关 Chemistry 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by