How to fit a noise signal in the shape of Lorentzian to a Lorentzian curve.

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));

I tried to fit the above curve using lorentzfit command but couldnt do it properly.

Answer 1

[Umar]

Hi @Yogesh,

First download this function from link provided below.

https://www.mathworks.com/matlabcentral/fileexchange/33775-lorentzfit-x-y-varargin?s_tid=srchtitle_support_results_1_lorentzfit

Once, you have downloaded the function, you can name the script as “ test function.m”. To fit a Lorentzian model to the data generated from your simulation, I will provide an updated code snippet that incorporates the necessary modifications. This will ensure you can effectively fit your data and visualize the results. The provided MATLAB code by you simulates optical fiber dynamics, including parameters such as fiber length, pump wavelength, and noise modulation. The end goal is to fit a Lorentzian function to the power spectrum obtained from your signal analysis. The lorentzfit function you included is designed for this purpose but requires proper integration into your simulation workflow. Here's how you can integrate the lorentzfit function into your existing code:

% Add this after calculating Power_FFt
figure; 
plot(nu, 2*n*c*eps0*A*(FFt_EL0t/Nt).^2, 'b-', 'LineWidth', 2);
title('Power Spectrum');
xlabel('Frequency (Hz)');
ylabel('Power Spectrum');

% Fit the Lorentzian function
[xData, yData] = deal(nu, 2*n*c*eps0*A*(FFt_EL0t/Nt).^2); % x and y data for   fitting

% Initial parameter guess: [P1, P2 (center frequency), P3 (width), C 
(constant)]
initialParams = [max(yData), mean(xData), 1e-6, min(yData)];

% Perform the fitting
[yFit, params, resnorm, residual] = lorentzfit(xData, yData, initialParams);

% Plotting the fit
hold on;
plot(xData, yFit, 'r-', 'LineWidth', 2);
legend('Original Data', 'Lorentzian Fit');
hold off;

% Displaying parameters
fprintf('Fitted Parameters:\n');
fprintf('P1: %.4f\n', params(1));
fprintf('P2: %.4f\n', params(2));
fprintf('P3: %.4f\n', params(3));
fprintf('C: %.4f\n', params(4));

Explanation of Changes

Data Preparation: The xData and yData variables are created to store the frequency and corresponding power spectrum data.

Initial Parameter Guess: An initial guess for the parameters[P1, P2 (center frequency), P3 (width), C (constant)] is provided. These are critical for ensuring that the optimization converges correctly.

Lorentzian Fitting: The lorentzfit function is called with the prepared data and initial parameters. It returns fitted values (yFit) optimized parameters (params) residual norm (resnorm), and residuals.

Visualization: A plot is generated to compare the original power spectrum with the fitted Lorentzian curve. This provides a visual representation of how well the model fits the data.

Parameter Output: The fitted parameters are printed to the console for review.

Please see attached.

Here are some additional insights that I would like to share with you.

Choosing Initial Parameters: Properly choosing initial guesses can significantly affect convergence in non-linear fitting problems. Adjust these values based on your specific dataset characteristics.

Error Handling: Consider implementing error handling around your fitting process to manage cases where fitting might fail or return non-physical parameters.

Validation of Fit: After fitting, evaluate how well your model describes the data using statistical measures like R² or visual inspection of residuals.

By integrating these changes into your MATLAB script, you should be able to fit a Lorentzian function to your simulated power spectrum effectively.

If you encounter issues during fitting or need further assistance with specific aspects of your simulation or analysis, feel free to ask!

[Yogesh]

Hey thanks for the reply @Umar but I cant help but notice the plot where lorentz fit is a straight line at the bottom of the curve. Shouldnt it follow the contour of the plot??..

[Umar]

Hi @Yogesh,

Please see attached. Is this what you were looking for.

[Yogesh]

Hey this was exactly what I was looking for... Thank you!!.. I'll try to replicate this for my data...

[Umar]

Hi @Yogesh,

Please see the attached updated code. If this helped resolve your issue, please don’t forget to click “Accept Answer”. Glad to know if this was the solution you were looking for.

[Yogesh]

Hey one more doubt, is using command lorentzfit enough or should I call the function as a seperate script!!?..

[Umar]

Hi @Yogesh,

My suggestion would be defining lorentzfit function defined per your original code at the end of the script to ensure it fits the generated sample data properly. Again, refer to my attached code snippets for reference. Hope this helps.