function [s_T, A, freqs] = FAR_signal_model_sample(FAR_N, FAR_M)
    global c f_c N N_wide B ranges_wide betas_wide FIGURE numP f_s N_high T_r T_p TEST method;

    pulse_number = FAR_N;
    single_pulse_length = N;
    all_pulse_length = single_pulse_length * pulse_number;
    pulse_width = N_high;
    
    delta_f = B / FAR_M;
    delta_r = c / (2 * B);
    delta_v = c / (2 * f_c * FAR_N * T_r);

    ranges = (0: FAR_M - 1) * delta_r;
    velos = (0: FAR_N - 1) * delta_v;
    delays = 2 * ranges / c;

    t_single = (0: 1/f_s: T_p-1/f_s)';
    t = (0: 1/f_s: FAR_N * T_r-1/f_s)';
    s_T = zeros(length(t), 1);
    T_p_N = floor(T_p * f_s);
    T_r_N = floor(T_r * f_s);

    freqs = zeros(FAR_N, 1);
    for i = 1: FAR_N
        C_n = floor(rand() * FAR_M);
        freqs(i) = f_c + C_n * delta_f;
        % freqs(i) = f_c + i * delta_f;
    end
    
    for i = 1: FAR_N
        N_l = (i - 1) * T_r_N + 1;
        N_r = (i - 1) * T_r_N + T_p_N;
        s_T(N_l:N_r) = exp(1j * 2 * pi * freqs(i) * t_single);
    end

    % echo
    s_R = zeros(length(s_T), FAR_M, FAR_N);
    s_D = zeros(length(s_T), FAR_M, FAR_N);
    for range_idx = 1: FAR_M % R
        r = ranges(range_idx);
        for doppler_idx = 1: FAR_N % D
            v = velos(doppler_idx);
            f_n = freqs(doppler_idx);
            for pulse_idx = 1: FAR_N % N
                nT_r = (pulse_idx - 1) * T_r;
                rr = r + v * nT_r;
                delay = 2 * rr / c;            
                t_l = delay + nT_r;
                N_l = max(1, round(t_l * f_s) + 1);
                N_r = N_l + T_p_N - 1;
    
                s_D(N_l: N_r, range_idx, doppler_idx) = exp(-1j * 2 * pi * f_n * (2 / c) * (r + v * (nT_r + t_single)));
            end
        end
    end

    % Sampling
    A = zeros(FAR_N, FAR_N * FAR_M);
    for range_idx = 1: FAR_M % R
        for doppler_idx = 1: FAR_N % D
            index = (doppler_idx - 1) * FAR_M + range_idx;
            for pulse_idx = 1: FAR_N % N
                t_sample = delays(range_idx) + T_p + (pulse_idx - 1) * T_r;
                A(pulse_idx, index) = s_D(round(t_sample * f_s), range_idx, doppler_idx);
            end
        end
    end
    A = get_Psi(FAR_N, FAR_M, 0);
    AAH = A * A';
    A = A ./ sqrt(AAH(1, 1));

    if 1 == 0
        Lambda = 1:10;
        x_wide = zeros(N_wide, 1);
        x_wide(Lambda) = 1;
        n_wide = randn(length(s_T_FAR), 1) * 0.01;
        y = A_wide * x_wide;
        y_noise = y + n_wide;
        x_hat = recovery(A_wide, y_noise, method);

        figure(1);
        subplot(211); plot(abs(x_wide));
        subplot(212); plot(abs(x_hat));

        figure(2);
        subplot(311); plot(real(A_wide * x_wide));
        subplot(312); plot(real(A_wide * x_hat));
        subplot(313); plot(real(A_wide * (x_wide - x_hat)));
    end