Spectrogram Analysis : An Approach to Identify The Quadrotor Thrust Model

This paper deals with a non-contact method to identify the aerodynamic propeller parameters of the Parrot AR.Drone quadrotor. The experimental set consists in a camera recording the vehicle flights, the audio signal is extracted and is used a spectrogram analysis to estimates the propeller velocity. First, the aerial vehicle takes off and starts a hovering maneuver. The experiment is repeated with different additional masses attached to its rigid body. If the weight over the UAV increases/decreases, then the propeller must rotates faster/slower to produce a higher/lower thrust, and consequently, the sound frequency increases/decreases. Finally, this proposition is validated experimentally, and the estimated velocity is used to identify the quadrotor thrust parameters.


Introduction
Recent advances in technology have generated a wide range of applications of Unmanned Aerial Vehicles (UAV).Most of them are related to the the increasing computational power of embedded systems and the improvement of sensors measurement.Rotorcraft navigation is still a challenge, which motivates research groups worldwide to develop autonomous systems.
Nowadays, rotorcraft UAVs have been a great testbed, due to its three-dimensional mobility when compared with ground or other aerial vehicles (airplanes and balloons, for instance).In contrast, such UAVs are inherently unstable, non-linear and multi-variable systems, with complex and highly coupled dynamics.Therefore, control them is another challenge, which starts from a suitable dynamic model [1].In our previous work [2], the quadrotor mathematical model, obtained by Euler-Lagrange equations, has been represented as a cascade connection of four interconnected subsystems, which compose the low-and high-level models.The low-level one represents the actuator model, which relates the input signals, the internal PD control loop, the Electronic Speed Controller, the motor dynamics, its angular velocity and the generated thrust.While the high-level one relates the rigid body dynamics, which includes the thrust applied on the UAV body and its pose response as displacement on 3-D space.
An autonomous flight requires reliable sensor information.For AR.Drone Parrot quadrotor, the UAV used in this work, these information are available from a inertial measurement unity (IMU), an ultrasound sensor and a built-in camera.Then, a filtering algorithm embedded on its firmware is responsible for determining its altitude, attitude, linear and angular velocities.Finally, this UAV can perform several flight tasks, including path planning and collision avoidance [3][4][5], but even with good sensors quality, it is necessary a reliable model, with tuned parameters, once the nonlinear controller depends on it.
The whole parameter identification of a quadrotor requires input signal information, rotor-speed values and 3-D pose during a non-constraint flight.For a rotor-speed estimation, one can use the mathematical model of the propeller system, and then validate it with a tachometer sensor.In [6] an identification technique determines the brushless DC-motor parameters of RC helicopters and quadrotors.Rotary encoder and hall sensors inform respectively the rotor-speed and the required current.In such a way, the propeller velocity is indirectly estimated.
Spectral analysis, or spectrogram, is another way to estimate an angular velocity.Such method avoids mechanical contact or the proximity needed by optical sensor and a specific equipment for measurements, only a microphone.According to [7], there is a proportional relationship between the audio frequency and the angular frequency.In such a context, the present work aims to estimate the propeller velocity of a rotary-wing 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 a room captures the acoustic signal, and with spectrogram analysis estimates the rotor speed.Thus, the propellers constants of the UAV model can be determined.This method avoids errors caused by the test platform, when the vehicle is attached to it.In such a case, the measurements are done by force sensors.
Figure 1 illustrates the top view of the room used in the experiments.The light color in the figure shows the useful flight area for the quadrotor, the hatched line shows the forbidden zone (due to room partial occupation).The computer is located in a table near the wall and the microphone is placed 2 meters above ground, in the position illustrated, to capture all the sound in the room.It important to stress that, even a standard microphone can measure the frequency and the intensity, due to the high intensity sound produced by UAV.This paper is hereinafter split as follows: First, a brief state of the art on spectral analysis is presented.In the sequel, the quadrotor dynamics model is introduced, emphasizing its rotary wing dynamics (also called thrust model).The experimental setup is described, and the results are shown and discussed.Finally, the main conclusions are highlighted and some suggestions for future works are also presented.

Spectrogram approach: State of the Art
The main advantage of spectral analysis to determine angular velocity is its non-contact approach.Moreover, this method does not require proximity to the equipment, as usual in optical or magnetic sensor applications.In contrast, the only restriction is the capture and interpretation of the sound signal of the source of interest.
The authors in [8] emphasize that a typical technique to measure the fundamental frequency of rotary machines requires proximity probes or photo-sensors which can be hard task.In this work, mechanical vibration and/or acoustic signals are used to exact the fundamental frequency of an automotive passenger car.A Bayesian approach combines the fundamental frequency and its harmonics and provides a better velocity estimation.
In the aerial robotic perspective, spectrogram and sound processing can be used to identify acoustic signature and classify moving UAVs [9].In such a case, the micro-Doppler signature classifier is capable to distinguish planes, quadrocopter, helicopters and stationary rotors.The patterns of each vehicle are extract form their spectral signature, and the classifier accuracy is around 95%. Sonogram is applied in [10], to obtain the relationship between distance sound level is correlated, and the sound profile in the propeller plane, for distinct rotor velocities.In [11], a collision avoidance technique is performed by acoustic signals, similar to conventional radars.In such a case, a fixed-wing UAV is equipped with a set of two microphones that detects and avoids other nearby UAVs.
The spectrogram analysis is also applied in [12] to determine the altitude, velocity and/or revolutions per minute of an aircraft, based on a static beacon signal recorder (ground-based microphone).In turn, the identification of UAV at long distance is the proposal presented in [13].Doppler effect estimates the number of revolutions per minute of the drone rotors, and the Plotted Image Machine Learning (PIL) and K-Nearest Neighbors (KNN) classify them, based on average distance similarity of a preset signature data-base.The identification accuracy are 83% and 61%, respectively.In [14], a coaxial UAV called Flyper helicopter, is mounted in a platform, where a set of microphones records its acoustic signal.Then, the rotor speeds are estimated without any prior knowledge about its acoustic properties.The proposed method based on neural network and genetic algorithm could identify each rotor velocity, without being affected (or corrupted) by the other.
Another perspective is to use a microphone array in the quadrotor to use in missions as "search and locate", or, sound source localization as it is called [15,16].In [17], the the potential of using the acoustic signatures for precision moving target detection and tracking.Four types of airplanes perform a set of predefined maneuvers over the microphone array, such as approaching, departing, turning left, and turning right.Then, spectrogram is performed to classify the airplanes and estimates its height, speed and fundamental frequency.
Acoustic technologies have also become an important tool for solving problems of UAV swarn navigation.In [18], the noise emitted by the propellers is used to detect and identify the UAV positions, minimizing the mid-air collision risk.In a distinct way, when the UAV carries the microphone array, another set of missions can be performed, such as "find and rescue" of humans in a disaster situation, using a quickly speech recognition [19,20].The UAV ego-noise is suppressed to emphasize, recognize and localize the voice signal (which generally occupy different time-frequency bins).

The Quadrotor Model
The complete model of a quadrotor can be represented by four interconnected subsystems [21], as shown in Figure 2 and the reference systems and the propeller forces acting on it are shown in Figure 3.The actuator dynamics is responsible for transforming the servo inputs into rotor-speed.The rotary wing dynamic 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 dynamic describes the aircraft displacement in the Cartesian space.
A quadrotor is an aircraft comprising four rotors symmetrically distributed around a rigid body.Its displacement is generated by changing the velocities of all motors (which are directly coupled to Actuator Dynamics

Rotary Wing Dynamics
Torque and Forces Generation Rigid Body Dynamics  the propellers).The profile of the propeller, associated to its rotational speed, generates a normal propulsion.As seen in Figure 3, the thrust points towards b z at the body reference frame, and is always positive.
For the AR.Drone Parrot quadrotor, the low-level system is responsible for its stabilization.In other words, an inner controller executes a hovering maneuver, whenever none control signal is sent to it.Thus, the first subsystem is the embedded low-level controller, composed by an inner control loop, normally a PD controller.Such block receives the joystick commands in the following order: pitch angle θ d , roll angle φ d , yaw rate ψd and vertical rate żd .Its output is the voltage increase applied to the motors, given by where φ = φ d − φ is the roll error.Similarly, we have the pitch error θ, yaw rate error ψ and vertical rate error z˜.Notice that k pi and k di are the proportional and the derivative positive gains of the embedded low level controller.
Remark 1.The commands sent by the joystick, u i , are normalized in ±1.Thus, the desired value of roll, for instance, is determined by φ d = u φ φ max , where φ max defines the maximum scale value (the same applies for the other low-level control signals).
Remark 2. The AR.Drone motors are not aligned with b x and b y axes.Actually, they are in a direction that is displaced 45 • of such axes.Thereby, to perform any lateral or longitudinal maneuver it is required a joint action of all engines, unlike other available quadrotors.
For the second subsystem, assuming that each AR.Drone brushless motor is modeled like as a brush DC one, one can describe from [22,23] the relationship among the input voltage v, the load torque τ l and the rotor speed ω as where R represents the rotor resistance, k m and k b are magnetic field flow constants, r is the gear ratio, and J m and B m are the moment of inertia and a dissipative terms of the motor.According to [23][24][25][26], the rotary wing dynamics (also called thrust model) can be written as where f i and τ li are the thrust and load torque generated by the i−th motor.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, the air density, among other phenomena.These parameters are considered as constants for restricted applications, such as non-aggressive maneuvers and low altitude flights, for instance.
The focus of this work is identify the C f coefficient using spectrogram analysis.In this current proposal, C τ estimation is not considered, because its computation requires torque measurements, which are not available in our lab.As comment the two first subsystems in Figure 2 are commonly labeled Low-level Model.The other two subsystems are the High-level Model, which describes the aircraft movement at the 3-D space.Details about the complete model can be found in [27].

The Experimental Setup
A set of dynamic equations can describe the whole AR.Drone Parrot model (see [27]), i.e., a representation through 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 low-level model, a force sensor could be used to relate straightforward the inputs u i and the thrusts f i .But, such sensor is commonly dedicated to specific application of robotic manipulators (force feedback control for welding tasks, for instance), besides being not available for usage in many research laboratories.Hence, the objective is estimate the rotor speed, and then determine indirectly the propeller thrust.In our approach, spectrogram analysis estimates the rotor speed, where the spectrogram is obtained by where x[n] is the discrete signal, w[n] is the discrete window function and N is the number of samples in the analysis window.The sound signature is measured during flight, and then a relationship between the UAV plus additional load (the hoods that came within the box of the ArDrone Parrot).In such a way, C f can be identified from (3) during a hovering flight.
According to [7], the audio signature is represented by where N is the revolution per minute, n is the number of blades, and ζ is the audio frequency.In our case, n = 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, also, there was some background noise, as people talking outside and running equipment,but it does not affect the experiments, since the analyzed frequency range is distinct from most of the environment noise.Furthermore, the technique used is robust against these noises.The proposed method considers that even a regular camera can measure the desired frequencies and intensities to estimate the aerodynamic constants.So all the experiments was filmed with a Nikon d7000 digital camera, and then the audio was extracted.The only caution needed is to know the frequency range of the equipment to better performance.In our proposal, the rotor speed varies according with the weigh.Maintaining the same altitude, the controller increases the control signals to compensate the extra weight.Thus the propellers rotate faster, and, consequently the thrust along the z-axis increases, keeping the vehicle in the same altitude.Three different configurations shown in Figure 4 are used:

Preprints
Without any hull: total weight of 413g, with 2000mA battery; Small hull: with additional 31g that yields 442g, i.e. 7.5% of its total mass; Big hull: with additional 59g that yields 462g, i.e. 14.3% of its total mass.
Although not being designed for load transportation, these was the items design by the manufacturer that has in any of the product boxes.
Remark 4. According to the manufacturer, the main characteristics of the AR.Drone propellers are: 4 brushless inrunner motors; 14.5 watt and 28,500 rpm in hovering mode; self-lubricating bronze bearings; low noise Nylatron gears for 8.625 propeller reduction; and specific high propelled drag for great maneuverability.
In order to achieve a good amount of data for estimation, three experiments was run for altitude control (hovering task), where the vehicle should keep at an altitude of 1m height, with each hull configuration.Each experiment was repeated five times, during 10s.Each experiment repetition is called a trial.The take off and the landing maneuver were removed from the original audio, once the part of interest is the hovering task.For each trial, a new set of smaller experiments is obtained by a rectangular window with 40% of overlap, totaling 25 samples of approximately 5s long.
Two methods are proposed to estimate the mean velocity of the rotors.Both of them require the fundamental frequency, its second and third harmonics to estimate the propeller speed.So, three bandpass filters are applied, creating three subsets.There is no need to use a narrow band or a high order filter, it only must contain the desired frequency.Figure 5 illustrates the database structure, highlighting the 75 available samples per experiment.
Maximum Trace and Gaussian approximation are proposed to estimate the C f constant, as describe in the following subsection.

Maximum Trace Method
For each one of the samples, Maximum Trace method consists in looking for the highest module of the FFT along each i-th time instant, and identify its correspondent frequency.Then, apply (6) to obtain the velocity.If only the fundamental frequency is considered, one can apply to get the frequency average ζ ¯, where ζ i is i-th frequency, for each i-th time instant, and n is the number of observations in the sample.If the fundamental and the second and third harmonics are taken into account, the standard deviation is also required for each interest frequency.
The weighted average consists into improve the velocity estimation using the fundamental frequency and its harmonics.In theory, the first harmonic has the double frequency of the fundamental, hence the obtained value should be divided to match the fundamental frequency.In our proposal, the standard deviation weights the estimation and also compensates any uncertainty or inaccuracies of the camera microphone.Thus, the speed estimation becomes where ζ ¯ω is the weighted frequency mean, H is the number of the harmonic (0 for the fundamental frequency).Applying ( 6) and (3) with ζ ¯ω, one has C f .

Gaussian Approximation Method
The second method consists in not only consider the highest values but a percentage of this value.For example, the 80% of the upper values are collected of each one of the three filters.Then, a Gaussian approximation given by represents each i-th time instant, where ζ ¯ is the frequency average and Gaussian centroid.Once more, applying ( 6) and ( 3), one can estimate C f .

Experimental Results and Discussion
Three different experiments were run five times, each one for a distinct hull configuration, during a hovering task.The AR.Drone quadrotor flies freely in the 3D Cartesian space, with its altitude controlled.The inner PID controller of the UAV guides its longitudinal and lateral movements.Figure 6 illustrates the specrogram of each experiment.Notice that all five trials were jointed to compute the spectral graph and to show it as a figure.Applying the proposed Maximum Trace approach, one can identify the highest values of each time window, as shown in Figure 8 for the fundamental frequency, for no-hull case.Table 1 presents the C f average value, as well the other parameter used to estimate it.Notice that C f = 28.5180× 10 −4 Ns 2 is the average value for all samples of all trials.For Gaussian approximation method, after finding the amplitude peak in a time instant, the frequency and its neighbor values with at least 80% are taken to determine the Gaussian parameters.Figure 9 illustrates the trace of the fundamental frequency and its harmonics, as well the fitted Gaussian curve.
The estimated C f with other other interesting values can be seen in Table 2.The final value of the aerodynamic constant, using a regular mean, is C f = 28.5290× 10 −4 Ns 2 .Remark 5. From 50% to 90% of the highest frequency (peak in amplitude), C f has negligible variations.This work presents a non-contact method to identify the aerodynamic propeller constants of the Parrot AR.Drone quadrotor.First, the experiments are filmed, and then the audio is extracted.In the sequel, the spectrogram is performed, and the frequency profile is analyzed.The experiments demonstrate that additional masses on the UAV require more thrust to keep the altitude on the preset value.Thus, one can conclude that there is a direct relationship between the rotorspeed and the sound frequency.Such assumption is used to determines indirectly the propeller aerodynamic coefficient.
The aerodynamic constant C f 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, the air density, among other phenomena.So, it is important a simple method of its identification for the development of the quadrotor dynamic model.In this paper, Maximum Trace and the Gaussian approximation method get aerodynamic constant with the same order of magnitude, C f = 28.5 × 10 −4 Ns 2 .And is worth remember that this constant will be different from the one used in these experiments, unless the blades are the same as well as the atmospheric properties.
Next steps of this research involve the whole identification of the AR.Drone low-level model.As shown in Figure 2, the rotor speed can be used to identify the parameters of the motor controller and its dynamic model separately.Another option is a grey box model that describes both set of parameters, and relates the joystick inputs with the rotor speeds.However the challenge is to discriminate, from the spectrogram, which acoustic signature correspond to which motor.

Figure 1 .
Figure 1.Top view of room used in the experiments.

Figure 3 .
Figure 3. 6-DOF CAD model of a quadrotor, including the reference frames and abstract control inputs (forces f i ) associated to it.The inertial, the spatial and the body reference frames are referred to as e , s and b , respectively.

Remark 3 .
The relationship between v and ω represented in (2) is applied to each motor of the aircraft, having v = v o + ∆v, where v o is the voltage value correspondent to a hovering flight.In other words, v o represents the voltage required to compensate the own AR.Drone weight through the propellers.

Figure 6 .
Figure 6.Altitude Flight with no hull, small and big one, respectively.

Table 1 .
Measured velocity values and the calculated C f with Maximum Trace Method.

Table 2 .
Measured velocity values and the calculated C f with Gaussian Method.