function testimage = test_image_impulse_response_bw_1(nx, ny, f, n)
% Compute ...
% TEST_IMAGE_IMPULSE_RESPONSE_BW_1(NX, NY, F, N) returns ...

    % Set pattern characteristics.
    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) <= sqrt(2)/n) = 1;
    h = h/sum(sum(h));

    % 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
                    if (rayon <= 0.20*rayonmax - max(cxscale, cyscale)/2)
                         tvalue(fyindex, fxindex) = 0.0;
                    else if (rayon <= 0.20*rayonmax + max(cxscale, cyscale)/2)
                              tvalue(fyindex, fxindex) = 1.0;
                         else if (rayon <= 0.40*rayonmax)
                                   tvalue(fyindex, fxindex) = 0.0;
                              else if (rayon <= 0.60*rayonmax - max(cxscale, cyscale)/2)
                                        tvalue(fyindex, fxindex) = 1.0;
                                   else if (rayon <= 0.60*rayonmax + max(cxscale, cyscale)/2)
                                             tvalue(fyindex, fxindex) = 0.0;
                                        else if (rayon <= 0.80*rayonmax)
                                                  tvalue(fyindex, fxindex) = 1.0;
                                             else
                                                       tvalue(fyindex, fxindex) = 0.0;
                                             end
                                        end
                                   end
                              end
                         end
                    end

                end

            end

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

        end

    end

end