Consider a N x n non-centered matrix Sigma(n) with a separable variance profile: Sigma(n) = (DnXnDn-1/2)-X-1/2/root n + A(n). Matrices D-n and (D) over tilde (n) are non-negative deterministic diagonal, while matrix A(n) is deterministic, and X-n is a random matrix with complex independent and identically distributed random variables, each with mean zero and variance one. Denote by Q(n) (z) the resolvent associated to Sigma(n) Sigma(n)*, i.e. Q(n) (Z) = (Sigma(n) Sigma(n)* - zI(N))(-1). Given two sequences of deterministic vectors (u(n)) and (upsilon(n)) with bounded Euclidean norms, we study the limiting behavior of the random bilinear form: u(n)* Q(n) (Z)upsilon(n) for all z is an element of C - R+, as the dimensions of matrix En go to infinity at the same pace. Such quantities arise in the study of functionals of Sigma(n) Sigma(n)* which do not only depend on the eigenvalues of Sigma(n) Sigma(n)*, and are pivotal in the study of problems related to non-centered Gram matrices such as central limit theorems, individual entries of the resolvent, and eigenvalue separation.