## 1 Introduction  
The rotor-bearing system is the main component of large-scale, rotatable equipment such as, e.g., aircraft engines and gas turbines, etc., and most of the rotor-bearing systems are significantly affected by nonlinearities [1]. In order to investigate the properties of nonlinear, rotor-bearing systems, mathematical models are important [2]. The mathematical models can generally be divided into two categories: the numerical model and the physical model [3]. For nonlinear, rotor-bearing systems, the systems are usually simplified into physical models based on mechanical or electrical-related theories [4, 5]. For example, [PERSON] et al. [6] established a 5 degree of freedom (5 DOF) physical model about the aircraft engine spindle dual-rotor system. [PERSON] and [PERSON] [7] introduced a dual rotor-ball and bearing-stator coupling dynamic system. However, those physical models are relatively simple to describe complex, nonlinear, rotor-bearing systems such as, e.g., large centrifugal compressors and steam turbines, etc., where strong nonlinearities have to be taken into account [8]. Moreover, in practice, it is also difficult to establish the physical model of a system due to its complexity and the lack of physical knowledge [9]. Fortunately, it is possible to build a data-driven model by using only input and output signals, which is widely applied in system analysis, design, and control for the advantages of efficiency and the capacity of tracking changes of the system [10]. In 1980s, the NARX (Nonlinear Auto-Regressive with Exogenous Inputs) model was introduced by [PERSON] as a new representation for a wide class of discrete, nonlinear systems [11]. The Volterra series model [12], the block-structured model [13], and many neural network architectures [14] can all be considered as subsets of the NARX model [15]. In practice, many systems have been investigated by using the NARX model [16; 17; 18]. For example, [PERSON] et al. [19] established a NARX model of aluminum plate with structure damage and then detected the location of damage by frequency domain analysis, which provided a theoretical support for structural damage detection in practical engineering. Besides, the NARX model can also be used in other industrial scenarios such as modeling the large horizontal axis wind turbine [20] and large horizontal lathes [21]. Considering the nonlinear, rotor-bearing system, [PERSON] et al. [22] established a Volterra series model by using date sets of the vibration response under different positions; based on this model, the fault modes were identified. [PERSON] et al. [23] developed the Volterra series model to detect the crack of the two discs in the rotor-bearing system. In general, a random signal is selected as the system input in the traditional NARX modeling process. This is because the random signal contains different frequency and amplitude characteristics, and by using this, different properties of the system can be tested when the knowledge about system parameters and model structures is absent [24]. However, in practice, it is impossible to produce random input excitations when identifying rotor-bearing systems [25]. To address the aforementioned issue, in this study, a multi-harmonic signal generated by a speed-up process is applied as the system input to establish the NARX model. Moreover, a new approach to determine the NARX model by using the speed-up harmonic signal is introduced. Many researchers have analyzed the rotor-bearing system based on the speed-up process. For example, [PERSON] et al. [26] and [PERSON] et al. [27] extracted fault features of the rotor-bearing system by using the speed-up process, based on which the fault modes were identified. In this paper, the approach of identifying the NARX model of the nonlinear, rotor-bearing system is proposed by using the speed-up harmonic signal.