An efficient Galerkin averaging-incremental harmonic balance (EGA-IHB) method is developed based on the fast Fourier transform (FFT) and tensor contraction to increase efficiency and robustness of the IHB method when calculating periodic responses of complex nonlinear systems with non-polynomial nonlinearities. As a semi-analytical method, derivation of formulae and programming are significantly simplified in the EGA-IHB method. The residual vector and Jacobian matrix corresponding to nonlinear terms in the EGA-IHB method are expressed using truncated Fourier series. After calculating Fourier coefficient vectors using the FFT, tensor contraction is used to calculate the Jacobian matrix, which can significantly improve numerical efficiency. Since inaccurate results may be obtained from discrete Fourier transform-based methods when aliasing occurs, the minimal non-aliasing sampling rate is determined for the EGA-IHB method. Performances of the EGA-IHB method are analyzed using several benchmark examples; its accuracy, efficiency, convergence, and robustness are analyzed and compared with several widely used semi-analytical methods. The EGA-IHB method has high efficiency and good robustness for both polynomial and non-polynomial nonlinearities, and it has considerable advantages over the other methods.