function test_image = test_image_frequency_response_bw_1(nx, ny, f, n)
% Compute ...
% TEST_IMAGE_FREQUENCY_RESPONSE_BW_1(NX, NY, F, N) returns ...

    % Set pattern characteristics.
    fmax  = 0.5*(min(nx, ny) - 1);
    fmin  = 0.0*(min(nx, ny) - 1);
    amin  = 0.10;
    amax  = 0.90;

    % Compute the coefficients
    % of the ideal 2D low-pass filter.
    [fx, fy] = freqspace(n*f, 'meshgrid');
    H = zeros(n*f, n*f);
    H(sqrt(fx.^2 + fy.^2) <= 1/f) = 1;
    h = fwind1(H, hamming(n*f));

    % Precompute
    % often used constants.
    coriginx = (1 + nx)/2;
    coriginy = (1 + ny)/2;
    cxscale  = 1/(min(nx, ny) - 1);
    cyscale  = 1/(min(nx, ny) - 1);
    foriginx = (1 + n*f)/2;
    foriginy = (1 + n*f)/2;
    fxscale  = cxscale/(f - 1);
    fyscale  = cyscale/(f - 1);
    rayonmin = 0.0;
    rayonmax = 0.5;

    % Coarse move along the X-axis.
    for cxindex = 1:1:nx

        % Compute the coarse position
        % along the X-axis.
        cxposition = cxscale*(cxindex - coriginx);

        % Coarse move along the Y-axis.
        for cyindex = 1:1:ny

            % Compute the coarse position
            % along the Y-axis.
            cyposition = cyscale*(cyindex - coriginy);

            % Fine move along the X-axis.
            for fxindex = 1:1:n*f

                % Compute the fine position
                % along the X-axis.
                fxposition = cxposition + fxscale*(fxindex - foriginx);

                % Fine move along the Y-axis.
                for fyindex = 1:1:n*f

                    % Compute the fine position
                    % along the Y-axis.
                    fyposition = cyposition + fyscale*(fyindex - foriginy);

                    % Compute the pattern.
                    if (abs(fxposition) > abs(fyposition))
                         rayon = abs(fxposition);
                    else rayon = abs(fyposition);
                    end
                    phase = pi*fmax*(rayon - rayonmin)^2/(rayonmax - rayonmin) + pi*fmin*(rayon - rayonmax)^2/(rayonmin - rayonmax);
                    fvalue(fyindex, fxindex) = (1 + sin(phase))/2;

                end

            end

            % Compute the pattern.
            cvalue = sum(sum(fvalue.*h));         % Low-pass filter, and sub-sample.
            cvalue = amin + (amax - amin)*cvalue; % Re-scale amplitude.
            cvalue = min(1, max(0, cvalue));      % Saturate amplitude.
            test_image(cyindex, cxindex) = cvalue;

        end

    end

end