MATLAB (Octave) code from "Solving self-mixing equations for arbitrary feedback levels: a concise algorithm" (https://scihub.bban.top/10.1364/AO.53.003723)

absolute_distance.m

% Example self-mixing signal generation for a modulated laser - Kliese et al., 2014
c = 299792458; % Speed of light in vacuum (m/s)
C = 2; % Feedback parameter
alpha = 4.6; % Linewidth enhancement factor
T = 29.3e-3; % Simulation time (s)
N = 1000; % Number of samples
beta = 0.1; % Laser power modulation coefficient
lambda0 = 845e-9; % Laser wavelength (m)
L = 0.024 + 0.5e-6; % Target distance (m)
deltaf = - 46e9; % Frequency modulation coefficient (Hz)
pm = 2.9; % Power modulation coefficient
t = 0:T/(N - 1): T; % Sample times
triper = T; % Triangle period (s)
tri = 1 + sign (mod (t/triper, 1) - 0.5).*(1 - 2* mod (t/triper, 1)); % Triangle waveform
phi0 = 4*pi*L*(1/lambda0 + deltaf * tri/c); % Round-trip phase samples
p = zeros (1, N); % Self-mixing signal samples
for i = 1 : N % Generate the synthetic self-mixing signal
  p(i) = beta * selmixpower (C, phi0 (i), alpha) + pm * tri(i);
end
plot (t, p); xlabel ('Time (s)'); ylabel ('Power Variation '); % Plot the results

absolute_distance_and_velocity.m

% Example self-mixing signal generation for a modulated laser with the target
% moving with a constant velocity - Kliese et al., 2014
c = 299792458; % Speed of light in vacuum (m/s)
C = 2; % Feedback parameter
alpha = 4.6; % Linewidth enhancement factor
T = 29.3e-3; % Simulation time (s)
N = 1000; % Number of samples
beta = 0.1; % Laser power modulation coefficient
lambda0 = 845e-9; % Laser wavelength (m)
L0 = 0.024; % Target start distance (m)
v = 300e-6; % Target velocity (m/s)
deltaf = - 46e9; % Frequency modulation coefficient (Hz)
pm = 2.9; % Power modulation coefficient
t = 0 : T/(N - 1): T; % Sample times
triper = T; % Triangle period (s)
tri = 1 + sign (mod (t/triper, 1) - 0.5).*(1 - 2* mod (t/triper, 1)); % Triangle waveform
phi0 = 4* pi *(L0 + v*t).*(1/lambda0 + deltaf * tri/c); % Round trip phase samples
p = zeros (1, N); % Self-mixing signal samples
for i = 1 : N % Generate the synthetic self-mixing signal
  p(i) = beta * selmixpower (C, phi0 (i), alpha) + pm*tri(i);
end
plot (t, p); xlabel ('Time (s)'); ylabel ('Power Variation '); % Plot the results

harmonic_motion.m

% Example self-mixing signal generation for harmonic motion - Kliese et al., 2014
C = 2; % Feedback parameter
alpha = 4.6; % Linewidth enhancement factor
T = 20e-3; % Simulation time (s)
N = 2000; % Number of samples
beta = 1; % Laser power modulation coefficient
lambda0 = 850e-9; % Laser wavelength (m)
L0 = 0.1; % Target nominal distance (m)
A = 2.5e-6; % Motion amplitude (m)
f = 100; % Motion frequency (Hz)
t = 0 : T/(N - 1): T; % Sample times
d = A*cos (2* pi*f*t); % Displacement samples
phi0 = 4* pi/lambda0 *(L0 + d); % Round-trip phase samples
p = zeros (1, N); % Self-mixing signal samples
for i = 1 : N % Generate the synthetic self-mixing signal
  p(i) = beta * selmixpower (C, phi0 (i), alpha);
end
plot (t, p); xlabel ('Time (s)'); ylabel ('Power Variation '); % Plot the results

selmixpower.m

% Functions for solving self-mixing equations - Kliese et al., 2014
function power = selmixpower (C, phi0, alpha) % Power level at a sample in time
    if (C<= 1.0)
        [phimin, phimax] = boundsweak (C, phi0);
    else
        [phimin, phimax] = boundsstrong (C, phi0, alpha);
    end

    excessphase = @(x)x - phi0 + C*sin(x + atan (alpha));
    % If the value at the left bound positive, then it will be very close to the solution.
    % If the value at the upper bound is negative, it will be very close to the solution.
    if (excessphase (phimin) > 0)
        excessphase (phimin)
        phi = phimin;
    elseif (excessphase (phimax) < 0)
        excessphase (phimax)
        phi = phimax;
    else
        phi = fzero (excessphase, [phimin, phimax]);
    end

    power = cos (phi);
end

function [phimin, phimax] = boundsweak (C, phi0) % Find search region when C <= 1
    phimin = phi0 - C;
    phimax = phi0 + C;
end

function [phimin, phimax] = boundsstrong (C, phi0, alpha) % Find search region when C >= 1
    persistent m; % Solution region number
    
    if isempty (m); m = 0; end

    % Calculate upper & lower values of m where solutions exist then ensure m is between them
    mlower = ceil ((phi0 + atan (alpha) + acos (1/C) - sqrt (C*C - 1))/(2*pi) - 1.5);
    mupper = floor ((phi0 + atan (alpha) - acos (1/C) + sqrt (C*C - 1))/(2*pi) - 0.5);

    if (m < mlower); m = mlower; end
    
    if (m > mupper); m = mupper; end
    
    phimin = (2*m+1)*pi + acos (1/C) - atan (alpha); % Trough
    phimax = (2*m+3)* pi - acos (1/C) - atan (alpha); % Peak
end

velocimetry.m

% Example self-mixing velocimetry signal generation - Kliese et al., 2014
C0 = 0.36; % Feedback parameter, nominal value
alpha = 4.6; % Linewidth enhancement factor
T = 0.8; % Simulation time (s)
N = 2 ^ 16; % Number of samples
beta = 1; % Laser power modulation coefficient
lambda0 = 850e-9; % Laser wavelength (m)
L0 = 0.1; % Target nominal distance (m)
fwhm = 1800; % Doppler signal FWHM (Hz)
fD = 10000; % Doppler peak frequency (Hz)
noise = 0.60; % Noise (RMS value)
t = 0:T/(N - 1): T; % Sample times
p = zeros (1, N); % Self-mixing signal samples
% Generate a realization of the random process, psi, from the power spectral
% density, S, corresponding to the random speckle process (scaled to have
% an RMS value of 1 after generating the realization of the random process).
f = 0 : 1/T:(N - 1)/T;
S = N* sqrt (2* sqrt (2* log (2)) / (sqrt (pi)*T* fwhm))* exp ( - 4* log (2)*(f-fD).^ 2/fwhm ^ 2);
psi = ifft (S.* exp (sqrt ( - 1)*2* pi* rand (1, N)));
phi0 = 4* pi*L0/lambda0 + angle (psi); % Round-trip phase samples
for n = 1:N % Generate the synthetic self-mixing signal
  p(n) = beta *abs(psi(n))* selmixpower (C0*abs(psi(n)), phi0 (n), alpha);
end
noisy_signal = p + noise/2 * randn (1, N); % Add synthetic noise to the self-mixing
fft_size = 256; [spectra, f2] = pwelch (noisy_signal, fft_size, [], fft_size, N/T); % PSD
plot (f2/1000, 10* log10 (spectra)); xlabel ('Frequency (kHz)'); ylabel ('P.S.D. (dBV ^ 2/Hz)');