47 lines
1.8 KiB
Matlab
47 lines
1.8 KiB
Matlab
function [ModulationFlag] = determineDensityModulation(psi, Params, Transf)
|
|
% Axes scaling and coordinates in micrometers
|
|
x = Transf.x * Params.l0 * 1e6;
|
|
dx = x(2)-x(1);
|
|
% Compute probability density |psi|^2
|
|
n = abs(psi).^2;
|
|
nyz = squeeze(trapz(n*dx,1));
|
|
|
|
densityProfile = sum(nyz,2);
|
|
|
|
% We need to first smoothen the density profile and subtract the smoothened profile from the
|
|
% original density profile to -
|
|
% 1. Remove the dominant low-frequency content.
|
|
% 2. Isolate the high-frequency components (like a sinusoidal modulation).
|
|
% FFT of the residual then highlights the periodic part, making the
|
|
% modulation easy to detect.
|
|
|
|
% Step 1: Smooth the density profile (Gaussian smoothing)
|
|
smoothedProfile = smooth(densityProfile, 10);
|
|
|
|
% Step 2: Compute the residual (original - smoothed)
|
|
residual = densityProfile - smoothedProfile; % We do this
|
|
|
|
% Step 3: Compute the Fourier Transform of the residual
|
|
N = length(residual);
|
|
Y = fft(residual);
|
|
P2 = abs(Y/N); % Two-sided spectrum
|
|
P1 = P2(1:N/2+1); % Single-sided spectrum
|
|
P1(2:end-1) = 2*P1(2:end-1); % Correct for the energy in the negative frequencies
|
|
|
|
% Step 4: Check for significant peaks in the Fourier spectrum
|
|
% We check if the peak frequency is above a certain threshold
|
|
threshold = 1E-3; % This can be adjusted based on the expected modulation strength
|
|
peakValue = max(P1);
|
|
|
|
if peakValue > threshold
|
|
ModulationFlag = true; % Indicates sinusoidal modulation
|
|
else
|
|
ModulationFlag = false; % Indicates otherwise
|
|
end
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
|