非线性光纤光学中分步傅里叶算法（SSFFT）的matlab代码实现

非线性光纤光学中分步傅里叶算法（SSFFT）的matlab代码实现

SSFFT（分布傅里叶算法）函数代码如下：

function [waveform, f_spectrum] = SSFFT_array(step_num, sum_distance, recard_step, GVD_array, N, input_waveform, tau, omega)% split-step fast Fourier transform
    input_waveform = input_waveform.';
    % tau = tau.';
    omega = omega.';
    nt = length(input_waveform);
    waveform = zeros(nt, recard_step + 1);
    f_spectrum = zeros(nt, recard_step + 1);
    uu_array = zeros(nt, recard_step + 1);
    inverl = step_num / recard_step;
    deltaz = sum_distance / step_num; % step size in z

    % ---set the initialization parameters
    beta2 = GVD_array(1); %The second order dispersion
    beta3 = GVD_array(2); %The third order dispersion

    uu = input_waveform ./ abs(max(input_waveform)); % normalized amplitude

    dispersion = exp(deltaz * (...
        (1i * beta2 / 2 * omega.^2) + ...
        (1i * beta3 / 6 * omega.^3))); %phase factor

    hhz = 1i * N^2 * deltaz; % nonlinear phase factor

    % **********[Beginning of MAIN Loop]**********
    % scheme:1/2N->D->1/2N; first half step nonliner
    x = 1;
    uu_array(:, x) = uu;

    for n = 1:step_num

        temp = fft(ifft(uu) .* dispersion);
        uu = temp .* exp(hhz .* abs(temp).^2);

        if mod(n, inverl) == 0
            x = x + 1;
            uu_array(:, x) = uu;
        end

    end

    waveform = abs(uu_array).^2;
    f_spectrum = fft(uu_array, [], 1);
    f_max = max(abs(f_spectrum(:, 1)));
    f_spectrum = abs(fftshift(f_spectrum, 1) / f_max).^2;

    waveform = waveform.';
    f_spectrum = f_spectrum.';
    % **********[End of MAIN Loop]**********
end

通过一个例子调用一下SSFFT：

clear
clc

step_num = 2000;
recard_step = 1000;
sum_distance = 10;
z_grid = linspace(0, sum_distance, recard_step + 1);

nt = 2^12; Tmax = 200; % FFT points and window size
dtau = (2 * Tmax) / nt; % step size in tau

%---tau and omega arrays
tau = (-nt / 2:nt / 2 - 1) * dtau; % temporal grid
omega = (pi ./ Tmax) .* [(0:nt / 2 - 1) (-nt / 2:-1)]; % frequency grid

figure(1);
U0 = exp(-1/2 * (tau).^2); % 输入为高斯脉冲
plot(tau, U0, 'b', 'LineWidth', 4);
axis([-20 20 -0.3 1.1]);
xlabel('\tau');
ylabel('U(\tau)');

[A, B] = SSFFT_array(step_num, sum_distance, recard_step, [1, 1], 2, U0, tau, omega);

figure(2);
subplot(1,2,1);
mesh(tau, z_grid, A);
axis([-20 20 0 sum_distance -inf inf]); view([0 90]);
xlabel('归一化时间');
ylabel('距离');
zlabel('|A(\tau)|^2'); % |A(τ)|^2
colormap(jet);

subplot(1,2,2);
mesh(fftshift(omega), z_grid, B / 2);
axis([-10 10 0 sum_distance -inf inf]); view([0 90]);
xlabel('归一化频率');
ylabel('距离');
zlabel('|A(\omega)|^2'); % |A(ω)|^2
colormap(jet);

输入高斯脉冲如图：

输入高斯脉冲如图：