Correct method for frequency shift and processing

Multiplication of your discrete time signal by a time varying complex exponential is the correct way to effect a circular shift of the digital spectrum (i.e. rotate the spectrum about the origin of the z-plane):

$$ B_s[n] = B[n] \cdot e^{j\pi \frac{f_B}{F_s/2} n }$$

I will speculate that the artifacts you are seeing are due to numerical problems.

With frequency rotators, numerical problems usually show up in one of two ways:

1. When using the complex exponential directly, most numerical implementations of $\sin()$ and $\cos()$ will start to behave poorly when their arguments get too far from the interval $[-2\pi, 2\pi]$
2. When using a phase increment complex number that is applied repeatedly, repeated multiplications can cause the resultant number to have a modulus that moves away from 1.

Here is some GNU Octave/Matlab code, that I wrote years ago, that generates a vector of values to effect frequency rotation that avoids these numerical problems. The output of vector $W$ of this function should be multiplied point by point with your signal samples vector.

function [W, ZN] = freq_rotator(F, FS, N, MODE, A)
% [W, ZN] = freq_rotator(F, FS, N, 'arith')
% [W, ZN] = freq_rotator(F, FS, N, 'arith', Z0)
% W = freq_rotator(F, FS, N, 'exp')
% W = freq_rotator(F, FS, N, 'exp', PHI_C)
%
% *** This is a looping implementation to assist in porting to C/C++/CUDA
% Generate a vector of complex values W, that when multiplied term by
% term with a time domain sequence, effects a cyclic DFT spectrum shift
% in the DFT frequency domain.
%
%     F  the freqeuncy shift desired, with values in the range [-FS/2, FS/2]
%        encompassing the complete output range of possible rotations.
%
%     FS the sample rate in samples/second.
%
%     N  the number of output values to generate, i.e. the length of W.
%
%     MODE is the algorithm to use 'arith' for arithmetic or 'exp' for
%        complex exponential.  The arithmetic mode is useful for streaming,
%        but is serial.  The exponential mode uses complex exponentials
%        to compute the terms.  This function does it serially, but it could
%        be parallelized easily.
%
%     Z0 the initial complex value to use as the start of the rotation sequence.
%        It will be normalized to unit modulus.  It allows successive calls to
%        this function to maintain continuous phase of the rotator.  If not
%        specified, it defaults to 1+0j.
%
%     ZN the complex value to use in a subsequent call to this function to
%        maintain continuous phase, if continuing the rotation sequence.
%
%     PHI_C is the desired phase angle of the center term in radians.  If not
%        specified, it defaults to 0 radians.

if (nargin < 4)
    MODE = 'arith';
end


W = zeros(1,N);
radian_phase_inc = pi/(FS/2).*F;

if (strcmp (MODE, 'arith'))
    % Arithmetic mode algorithm. Good for streaming, but strictly serial,
    % because it's recursive:
    %     z[n] = z_phase_inc * z[n-1]

    if (nargin < 5)
        z0 = 1+0j; % exp(1j*0) because why not.
    else
        z0 = A/abs(A); % normalize to ensure it is on the unit circle.
    end

    z_phase_inc = exp(1j*radian_phase_inc);

    % N.B. all values in W should be unit modulus, i.e. on the unit circle
    W(1) = z0;
    j = 0;
    for i = 2:N
        W(i) = W(i-1)*z_phase_inc; % mult. of unit vectors adds phase angles.

        % renormalize once in a while, so accumlated error doesn't move us too
        % far away from the unit circle.
        j = j+1;
        if j == 256
            W(i) = W(i)/abs(W(i));
            j = 0;
        end
    end

    ZN = W(N)*z_phase_inc; % mult. of complex unit vectors adds phase angles.
    ZN = ZN/abs(ZN); % normalize to ensure it is on the unit circle.
else
    % Exponential mode algorithm - can be parallelized, but calls sin() & cos().
    % Conceptually this algorithm does this (ignoring index shifts) and can
    % be parallelized on the index n:
    %
    %   z[n] = exp(j*n*theta)*exp(j*phi_c) = 
    %             (cos(n*theta) + j*sin(n*theta))*(cos(phi_c) + j * sin(phi_c))
    %
    % There is a symmetry around the center value, so only half the terms
    % actually need to be computed by calling cos() & sin().  The other half
    % of the terms are obtained by reflection and conjugation.

    if (nargin < 5)
        w_phase_shift_vec = 1+0j; % i.e. no additional phase shift.
    else
        % Ensure within [0,2*pi).
        % When porting to C/C++/CUDA, take care with how mod() handles
        % negative signs. See the Octave or MatLab documentation.
        w_phase_shift_vec = exp(1j*mod(A, 2*pi));
    end

    if (mod(N, 2) == 0)
        center = N/2+1;
        neg_count = N/2;
        pos_count = N/2-1;
    else
        center = (N-1)/2+1;
        neg_count = (N-1)/2;
        pos_count = (N-1)/2;
    end

    % Center value.
    W(center) = 1+0j;
    % Add user requested phase shift.
    if (nargin >= 5)
        W(center) = W(center) * w_phase_shift_vec;
    end

    for i = 1:neg_count
        % Compute phase angle for negative offset from center.
        % Ensure within [0,2*pi).
        % When porting to C/C++/CUDA, take care with how mod() handles
        % negative signs. See the Octave or MatLab documentation.
        phase_angle = mod(-i*radian_phase_inc, 2*pi);

        % Compute required complex exponential values from phase angle

        % Negative offset from center.
        W(center-i) = exp(1j*phase_angle);

        % Positive offset from center.
        if (i <= pos_count)
            % We're symmetrical around the center point, so just conjugate and
            % copy for the positive offset from center.
            W(center+i) = conj(W(center-i));
        end

        % Add user requested phase shift.
        if (nargin >= 5)
            % Negative offset from center.
            W(center-i) = W(center-i) * w_phase_shift_vec;
            % Positive offset from center.
            if (i <= pos_count)
                W(center+i) = W(center+i) * w_phase_shift_vec;
            end
        end
    end
    ZN = 0; % meaningless output
end

The shortest shorthand notation for frequency shifting a real signal $x(t)$ by an amount $f_B$ is:

$$ y(t) = \Re e \Big\{ \big(x(t) + j \hat{x}(t) \big) e^{j 2 \pi f_B t} \Big\} $$

where

$$\begin{align} \hat{x}(t) & \triangleq \mathscr{H} \Big\{ x(t) \Big\} \\ \\ & = \int\limits_{-\infty}^{\infty} \frac{1}{\pi u} \ x(t-u) \ \mathrm{d}u \end{align}$$

And $\mathscr{H} \{\cdot\}$ is the Hilbert transform.

Now, to do this digitally (in discrete time) means you need to make a discrete-time Hilbert transformer. That means turning an infinite impulse response into a finite impulse response (suggest using the Kaiser window) and sampling it. This FIR (which is centered at $t=0$ and therefore acausal) must be delayed by a known amount to make it causal so what you get out of the Hilbert transformer is really:

$$ \hat{x}(t- T_D) = \int\limits_{-T_D}^{+T_D} \frac{1}{\pi u} w(u) \ x(t - T_D - u) \ \mathrm{d}u $$

$w(t)$ is the Kaiser window centered at $t=0$. Since the Hilbert output, $\hat{x}(t-T_D)$ is delayed by $T_D$, so also should $x(t)$ be delayed by the same amount for this to work. So you get instead:

$$ y(t) = \Re e \Big\{ \big(x(t-T_D) + j \hat{x}(t-T_D) \big) e^{j 2 \pi f_B t} \Big\} $$

In this answer, I go through what you need to do to make a discrete Hilbert transformer.

First of all I'm assuming you are only interested in a band-limited section of signal B around DC of bandwidth $B_i$.

If $|f_B| + B_i/2 < F_s/2$ then you can just mix and BP filter, otherwise you will have the situation where your signal of interest occupies sections of both the positive and negative digital bandwidths (i.e., it gets split across the image.

In that case you would first need to resample the input signal B to a higher sample rate $F_{s2}$ so that $|f_B| + B_i/2 < F_{s2}/2$.

Notes:

1. $f_B$ can be either positive or negative.
2. the BP filter will be complex
3. the output will be complex
4. I'm unsure of the relevence of signal A: you define it but its use in your question is not clear.