This paper presents an original algorithm to simulate quasi-static loading of fluid-saturated porous material for numerical upscaling of the complex fluid-saturated fractured-porous media. The resulting model is anisotropic viscoelastic which means that it is defined by a complex-valued frequency-dependent stiffness tensor that can be recomputed into the frequency-dependent phase velocities and attenuation of seismic waves. Numerical upscaling requires a solution of parabolic approximation of the Biot poroelastic equation for anisotropic media in frequency space for a series of frequencies and multiple right-hand sides. The approach is based on the finite-difference approximation of the quasi-static formulation of the Biot equation with the use of a direct solver to resolve the obtained system of linear algebraic equations. The direct solver allows for efficiently treating multiple right-hand sides which is an essential statement of the upscaling problem. The presented realization of the algorithm allows solving the problems of the size of up to 2000 points in each spatial direction using a single machine, which allows dealing with representative fractured-porous models. To illustrate the applicability of the algorithm we provide several series of experiments illustrating the effects of fracture connectivity and anisotropy of fracture-filling material on seismic waves dispersion and attenuation if propagating in complex fluid-saturated fractured-porous media.