‎In this paper‎, ‎a numerical method based on cubic hat functions is presented for solving a class of fractional stochastic integro-differential equations driven by multidimensional fractional Brownian motion with Hurst parameter‎ ‎$ H \in \left(\frac{1}{2}‎ , ‎1\right) $‎.‎The fractional derivative is considered in the Caputo sense‎. ‎First‎, ‎an approximation of the Caputo derivative of the unknown function with unknown coefficients is constructed based on the basis functions‎. ‎Then‎, ‎using this approximation and the properties of the Caputo derivative and the Riemann-Liouville integral‎, ‎approximations for the solution and other functions appearing in the equation are introduced‎. ‎By substituting these approximations into the equation and using suitable nodal points along with a backward approximation for simulating the fractional Gaussian noise‎, ‎the original equation is transformed into a system of algebraic equations‎. ‎Due to the properties of the cubic hat functions and their operational matrix of integration‎, ‎the resulting system has low computational complexity‎. ‎Solving this system yields an approximation of the solution to the original problem‎. ‎Furthermore‎, ‎the error analysis of the proposed method is investigated‎. ‎Finally‎, ‎several numerical examples are provided to confirm the efficiency of the proposed method‎.