UAV Thrust Model Identiﬁcation Using Spectrogram Analysis

: This paper deals with a non-contact method to identify the aerodynamic propeller constants of the Parrot AR.Drone quadrotor. The experimental setup consists of a microphone installed in the ﬂight arena to record audio data. In terms of methodology, a spectrogram analysis is adopted to estimate the propeller velocity based on the ﬁltered sound signal. It is known that, in a hovering maneuver, when the UAV mass increases, the propellers rotate faster to produce the necessary thrust increment. In this work, the rotorcraft takes off with its factory settings, ﬁrst with no hull, corresponding to a mass of 413 g, and after with a small hull, corresponding to a mass of 444 g, and a bigger hull, corresponding to a mass of 462 g. In the sequence, the velocity of the propellers are estimated for each of these three cases using spectrograms of audio recorded by a microphone, corresponding to the sound generated by the four rotors. Finally, the estimated velocity is used to identify the aerodynamic parameters, thus validating the proposal. and 25 M.Sc. students. His research interests are nonlinear control, mobile robot navigation, multi robot systems, coordinated control of mobile robots, unmanned aerial vehicles and coordinated control of ground and aerial robots.


Introduction
Recent advances in technology have generated high-performance Unmanned Aerial Vehicles (UAVs), which have been used in a wide range of applications. Most such advances are related to the increasing computational power of embedded systems and the improvement in sensor accuracy. However, UAV navigation is still a challenge, motivating research groups worldwide to develop UAV navigation systems.
Rotorcraft UAVs, in particular, have become a great testbed, due to its three-dimensional mobility when compared with ground vehicles or even other aerial ones (airplanes and balloons, for instance). However, such an advantage has a price: they are inherently unstable nonlinear multi-variable systems with complex and highly coupled dynamics [1].
The design of a suitable controller for such UAVs starts from a dynamic model representing them. However, precise modeling of the UAV model is also a challenge. In [2], the mathematical model of a quadrotor, obtained using Euler-Lagrange equations, is represented as a cascade connection of four interconnected subsystems, which compose the lowand high-level models. The low-level part represents the actuator model (which involves the input signals, the internal PD control loop, the electronic speed controller, the dynamics, and the angular velocities of the motors and the generated thrust), whereas the high-level one corresponds to the rigid body dynamics (which considers the thrust applied to the UAV body and its 3-D displacement. Nonetheless, to automatically control the flight of a UAV using a controller designed considering its dynamic model, reliable sensory information is necessary. For the Parrot AR.Drone quadrotor, the UAV used in this work, this information is provided by a inertial measurement unity (IMU), an ultrasound sensor, and a built-in camera. Then, a Kalman filter algorithm embedded on its firmware is responsible for determining its altitude, attitude, and linear and angular velocities by performing a data fusion procedure. As a result, such a UAV can accomplish several flight tasks, including collision avoidance, as described in [3][4][5].
The whole parameter identification for a quadrotor requires input signal information, rotor-speed values, and a 3-D pose during a non-constrained flight. In the case of rotor-speed estimation, one can use the mathematical model of the propeller system and then validate it with a tachometer sensor, for instance. This is illustrated in [6], where an identification technique allows for estimating the brushless motor parameters of RC helicopters and quadrotors. A rotary encoder and a hall sensor inform about the rotor-speed and the electric current supplied to reach such speed, respectively, thus making it possible to indirectly estimate the velocity of the propellers.
Using spectrograms is another way to estimate rotational velocities, with the advantage of not demanding the mechanical contact or proximity necessary when using optical sensors or other measuring equipment. According to [7], there is a proportional relationship between the audio frequency and the angular frequency. An example is discussed in [8], where acoustic signals and mechanical vibration produced by a motor are recorded and then analyzed. In the case of aerial robotics, spectrogram and sound processing can be used to identify the acoustic signature and to classify UAVs [9]; to determine the altitude, velocity, and revolutions per minute of an aircraft [10]; and to estimate the UAV location based on beacon signals [11]. In this last case, the strategy adopted allowed us to identify the UAV at long distances and to determine the number of revolutions per minute of the rotors as well, using the Doppler effect. Moreover, in [12], a microphone array was used to track a light-weight aircraft and to measure its variables.
Another perspective is to use a microphone array in the quadrotor for use in missions involving "search and locate" or sound source localization. In such applications, however, one certainly deals with a noisy environment, since the microphone array is mounted in the quadrotor. The same happened in the cases reported in [13,14], which used a microphone array to detect humans in disaster situations in missions of "search and rescue". Reference [15] proposed a compass made of microphones that are used for sound source localization.
Acoustic technologies have also become an important tool for solving problems related to the navigation of UAV swarms. In [16], the noise emitted by the propellers was used to detect and identify the position of the UAVs, minimizing the risk of mid-air collision (it is well known that there is an imminent risk of collision when UAVs share the airspace).
The sound signature of a UAV was studied in [17]. The Flyper autonomous helicopter was mounted on a platform, where a set of microphones recorded its acoustic signal. Then, the rotor speed information was extracted and evaluated without any prior knowledge about its acoustic properties. In such a case, the genetic algorithm identifies frequency bands while the neural network estimates the rotor speed.
In a similar way, the present work aims to estimate the propeller velocity of a rotarywing vehicle, the AR.Drone Parrot quadrotor. However, in our approach, the UAV flies without constraints in an indoor environment. A single stationary microphone installed in the flight arena captures its acoustic signal, and spectrogram analysis is used to estimate the rotation speeds of its rotors. Thus, the constants related to the propellers, associated with the low-level model of the quadrotor, can be determined. Notice that this is a noncontact measurement, which avoids risk to the user when trying, for instance, to manually measure the rotor speed, and allows us to obtain the value of such parameters even when dealing with vehicles that do not provide such information through the embedded sensors (close vehicles).
To deal with this subject, this paper is hereinafter split in a few sections, starting with Section 2, which briefly presents the quadrotor dynamic model, emphasizing its rotary wing dynamics (thrust model). Subsequently, Section 3 describes the experimental setup adopted, whereas Section 4 shows and discusses the results obtained. The analysis of the acoustic spectrogram is also presented. Finally, Section 5 highlights the main conclusions and gives suggestions for the continuity of the research. Figure 1 illustrates the reference frames and the propeller forces acting on the quadrotor. Its complete model can be represented by four interconnected subsystems [18], as shown in Figure 2. The actuator dynamics is responsible for transforming the servo inputs into rotor speed. The rotary wing dynamics relates the aerodynamic parameters to the thrust/propulsion. The force/torques generation decomposes the thrust into the 3-D forces/torques actuating on the UAV rigid body. Finally, the rigid body dynamics describes the aircraft displacement in the Cartesian space.

The Dynamic Model of the Quadrotor
The two first subsystems in Figure 2 are commonly labeled the Low-level Model. They receive the control signals and transform them into real forces applied to the aircraft, with all of them in the normal direction (see Figure 1). The other two subsystems are the High-level Model, which describes the aircraft movement in 3-D space. More details about the complete model can be found in [19]. Six-DOF CAD model of a quadrotor, including the reference frames and abstract control inputs (forces f i , i = 1, · · · , 4) associated with it. The inertial, the spatial, and the body reference frames are referred to as e , s , and b , respectively. The rotary wing dynamics (also called thrust model) is the focus of this work. According to [20][21][22], these dynamics can be approximated and written as where f i and τ li are the thrust and load torque generated by the ith motor, respectively.
Moreover, C f and C τ are the aerodynamic constants, which are dependent on the number, the width and the shape of the rotor blades, the inner and outer radius of the airflow through the rotor, and the air density, among other phenomena. These parameters are considered constants for restricted applications, such as non-aggressive maneuvers and low altitude flights.
In this work, however, the estimation of C τ is not considered because its computation requires torque measurements, which are not available. The constant C f depends only on the force exerted on the UAV, which can be easily obtained considering the UAV hovering and knowing the mass of the vehicle, as described in a later section.

The Experimental Setup
A set of dynamic equations can describe the whole AR.Drone Parrot model (see [19]), i.e., the vehicle can be represented using a white box model. However, the parameter identification of each subsystem shown in Figure 2 requires the input signals, the rotor speeds, the thrust produced by each propeller, and the 3-D pose of the UAV. For the lowlevel model, a force sensor could be used to straightforwardly relate the inputs u i and the thrusts f i . However, such a sensor is commonly dedicated for specific application of the robotic manipulators (force feedback control in welding tasks, for instance), besides not being available for use in many research laboratories. Hence, the objective is to estimate the rotor speed and, then, to indirectly determine the thrust of the propellers.
In our approach, the spectrogram analysis estimates the rotor speed. Then, the first step is to obtain the spectrogram, which is obtained as where x[n] is the discrete signal, w[n] is the discrete window function adopted, and N is the number of samples in the window. The sound signature is measured during flight, and then, a relationship between the mass of the UAV plus additional load (the hoods inside the box of the Parrot Ar.Drone quadrotor) and the thrust is obtained, allowing us to write (1), from which C f can be estimated.
According to [7], the audio signature is represented by where RPM is the number of revolutions per minute, b is the number of blades, and ζ is the audio frequency. In our case, b = 8, once each propeller has two blades, and all of them contributes constructively to the audio signature. For the AR.Drone Parrot, the propeller gear reduction r is 8.625. The experiments took place in a regular laboratory with no concern about the acoustic performance of the room. There was also some background noise, such as people talking outside and the sound of the equipment running. The proposed method considers that even a regular camera can measure the desired frequencies and intensities to estimate the aerodynamic constants. Therefore, all of the experiments were recorded with a digital camera and the audio was extracted.
In our proposal, the rotor speed varies according to the weight, increasing somewhat to maintain the same altitude, with the addition of extra mass. Then, the controller increases the control signals to compensate for the extra weight. Thus, the propellers rotate faster, and consequently, the thrust along the z-axis increases, keeping the vehicle hovering at the same altitude.
Three different configurations were adopted, regarding the mass of the quadrotor, as illustrated in Figure 3.
In the first case, no hull is added to the quadrotor, as Figure 3a shows, resulting in a total mass of 413 g, with a battery of 2000 mA included. In the second case, illustrated in Figure 3b, the small hull designed for the quadrotor is added, in which the mass is 31 g (7.5% of the total mass of the vehicle), resulting in a total mass of 444 g. For the third case, the bigger hull designed for the quadrotor (see Figure 3c) is installed, adding 59 g to the 413 g of the vehicle (a mass increment of 14.3%), resulting in a total mass of 462 g.  As additional information, the AR.Drone data sheet gives its main characteristics. It has 4 brush-less in-runner motors, consuming 14.5 Watts and working at 28,500 rpm in hovering mode; self-lubricating bronze bearings; low noise Nylatron gears with 8.625:1 speed reduction; and specific, high-propelled drag for great maneuverability.
In order to achieve a good amount of data to perform a reliable estimation, three experiments of altitude control were run, in which the vehicle was controlled by an external computer to maintain a hovering altitude of 1 m for no hull, the small hull, and the big hull installed in the quadrotor (see Figure 3). Each experiment was repeated five times, with a standard duration of approximately 10 s, thus generating five trials per experiment or 15 trials in total. All of the trials were recorded with a digital camera, and the correspondending audio was extracted. To guarantee regular values, the parts of the audio correspondending to take off and landing were removed, resulting in about 8 s of useful audio per trial. Then, the audio correspondending to each trial was split in five segments, each one containing 60% of the original trial, which corresponds to about 5 s, with 80 % of superposition with the previous segment, resulting in 25 sets of audio signal for each experiment. Figure 4 shows the spectrogram corresponding to one such audio sets, whereas Figure 5 shows the structure of the resulting database.
Two methods are proposed to estimate the mean velocity of the rotors from the audio recorded, both using the fundamental frequency of the propellers, and its second and third harmonics. Then, three band-pass filters are applied, creating three subsets, one for each filter. There is no need to use a narrow band or a high-order filter, it only must contain the desired frequency. Therefore, each one of the 25 audio sets generated three, leading to 75 audio subsets per experiment.    In this work, two simple methods are proposed to estimate the C f constant, the Maximum Trace, and the Gaussian approximation, which is made to allow for a comparison, since using one of them is enough to estimate the C f parameter.

The Maximum Trace Method
For each one of the audio subsets, the Maximum Trace method consists of finding the highest magnitudes of the FFT along each time window and of identifying in which frequency it occurs. Then, using the frequency to velocity transformation in (4), a velocity vector is obtained. If only the fundamental frequency is considered, one can usē to find the frequency meanζ, where ζ i is the ith entry of the frequency vector. The case with the fundamental, second, and third harmonics are considered; it is used a weighted mean, with the frequencies weighted by their respective standard deviation (STD), this last one given by where n is the number of observations in the sample. In the weighted mean approach, harmonics are used to improve the velocity estimation. As the first harmonic has twice the frequency of the fundamental one, it must be divided to match the same frequency as the fundamental one. Here, the standard deviation is used as a weight to compensate for any variation in the camera microphone. Thus, the rotor speed becomesζ whereζ ω is the weighted frequency mean and H is the number of the harmonic, with H = 1 for the fundamental one.
Having the mean frequency of the propellers, one has to just solve (1) and (4) for C f .

Gaussian Approximation Method
The second method adopted consists in considering not only the highest magnitudes of the signal FFT but also a percentage of this value. In this work, the 80% upper values are collected for the output of each filter. Then, the Gaussian approximation is adopted. The valueζ in (8) is the frequency mean and centroid of the Gaussian, so that using (1) and (4), one can obtain C f .

Experimental Results and Discussion
As aforementioned, three experiments were undertaken. For all of them, an altitude controller was used to keep the quadrotor hovering at the desired altitude. Nonetheless, longitudinal and lateral oscillations occur since there is no external active controller for such variables, so that only the embedded PID controller of the vehicle takes care of them, with a zero value as the setup.
The spectrograms for the three experiments can be seen in Figure 6. The duration is a little more than 40 s because all five trials of each experiment were added to calculate the spectrogram. Furthermore, an attentive reader can see the slight difference between the trials, more easily seen at higher frequencies, showing constant frequency values for the same propeller force. Additionally, comparing the spectrograms, one can see an increase in the frequency with the mass increase, corresponding to an increased propeller speed. Indeed, increasing the total mass of the vehicle, the natural response of the controller is to increase the thrust to compensate for the additional mass, i.e., the velocities of the propellers increase evenly, to keep the vehicle hovering at the same altitude.  The resulting spectrogram of an audio set is shown in Figure 7a, with the result of the implemented filtering in Figure 7b.

The Maximum Trace Method
Using the first approach, the Maximum Trace, one has to identify the highest values of the FFT magnitude for each time window (see Figure 8). The average value of the estimated C f together with other interesting values for each experiment is shown in Table 1. The final value of the aerodynamic constant, using a regular mean, is C f t = 0.028518 Ns 2 .

Gaussian Approximation Method
The Gaussian approximation method was used taking the magnitude values above 80% of the highest value for each time window. The result for a whole set is shown in Figure 9. Each subplot shows the result of the three subsets, one for each band-pass frequency. The frequencies were merged, considering that the velocity found with the first harmonic should be divided by two, whereas for the second harmonic, it should be divided by three. |Power/frequency| (dB/Hz) Figure 9. Gaussian approximation method. Each subplot shows the result of the application of the three bandpass filters to an audio set obtained as illustrated in Figure 5, for the case of Figure 3a.
The blue, yellow, and green lines correspond to the readings for the fundamental, first, and second harmonic frequencies, respectively. The red line corresponds to the Gaussian estimation of such frequencies.
The estimated C f together with other values of interest are shown in Table 2. The final value of the aerodynamic constant using a regular mean is C f g = 0.028529 Ns 2 . An interesting remark is that, even when using 50% or 90% of the highest magnitude values, instead of the 80% used to estimate the C f value, the variation is negligible. The difference between the results obtained with the two methods, in what concerns C f , is |C f t − C f g | = 1.110 −5 , which is absolutely negligible. This means that any of the two methods could be used without any special concern.
Another worthy remark is that, although a big amount of data was used to estimate the parameter C f , using only one random subset, one obtains a value of C f = 0.028378 Ns 2 , which is quite a good approximation.

Concluding Remarks
This work presents a non-contact method to identify the aerodynamic propeller constants of the quadrotor Parrot AR.Drone. The procedure involves two great steps, namely recordings of the sound generated by the rotors of a vehicle hovering at a fixed controlled altitude, and the computation and analysis of the frequency profile. The results obtained demonstrated that masses added to the UAV body demand more thrust from the propellers to maintain the preset altitude. Therefore, it is verified that there is a direct relationship between the rotor speed and the sound frequency, which is used to determine the aerodynamic coefficient for the propellers.
Two approaches are proposed to analyze the frequency profile obtained, namely the Maximum Trace method and the Gaussian approximation, and it is verified that both yield the same result. It is worth emphasizing that a mathematical determination of the aerodynamic constant C f is not easy, since it depends on the number, the width and the shape of the rotor blades, the inner and outer radius of the airflow through the rotor, and the air density, among other phenomena. Thus, a simple method to estimate it, as the one proposed here, is useful to obtain a reliable model for the UAV, mainly when the vehicle does not provide information on the rotor speed, which is common regarding on-the-shelf quadrotors.
To continue this research, the authors intend to perform whole identification of the AR.Drone low-level model, which includes estimation of the aerodynamic constant C τ . As shown in this text, the rotor speed can be used to identify the parameters of the motor controller and its dynamic model separately. Another option is a gray box model that describes both sets of parameters and relates the joystick inputs to the rotor speeds. However the challenge is to discriminate, from the spectrogram, which acoustic signature corresponds to which motor.