In this paper, an efficient numerical technique for solving fractional stochastic differential equations driven by fractional Brownian motion is provided. The proposed method is based on the cubic hat functions and their operational matrix of fractional integration. Properties of these basis functions together with a set of grid points help us to reduce the problem to a system of algebraic equations. An example and its numerical results demonstrate the effectiveness of the proposed method.