<jupyter_start><jupyter_text>In this demo we will analyze the mechanical vibration of a quadrotor during flight using publicly available logs.We obtain the logs from review.px4.io, https://review.px4.io/plot_app?log=8aa57955-db73-46d5-8800-17903bbc3575  Open the file<jupyter_code>%pylab inline

import pyulog

with open('data/8aa57955-db73-46d5-8800-17903bbc3575.ulg', 'rb') as f:
    log = pyulog.ULog(f)<jupyter_output>Populating the interactive namespace from numpy and matplotlib<jupyter_text>List the Data Topics<jupyter_code>for data in log.data_list:
    print(data.name)<jupyter_output>vehicle_local_position_setpoint
vehicle_magnetometer
vehicle_air_data
ekf2_timestamps
vehicle_status_flags
vehicle_status
vehicle_rates_setpoint
vehicle_local_position
vehicle_land_detected
vehicle_gps_position
vehicle_global_position
vehicle_attitude_setpoint
vehicle_attitude
telemetry_status
system_power
sensor_preflight
sensor_combined
safety
rate_ctrl_status
manual_control_setpoint
input_rc
home_position
estimator_status
ekf2_innovations
distance_sensor
cpuload
battery_status
actuator_outputs
actuator_outputs
actuator_controls_0<jupyter_text>Plot the raw x acceleration, about the aircraft roll axis<jupyter_code>sensor = log.get_dataset('sensor_combined')
accx = sensor.data['accelerometer_m_s2[0]']
t = sensor.data['timestamp']/1.0e6
plt.figure(figsize=(10, 5))
plt.plot(t, accx)
plt.gca().set_ylim(-10, 10)
plt.grid()
plt.xlabel('t, sec')
plt.ylabel('accel m/s^2')<jupyter_output><empty_output><jupyter_text>Compute the Sample Time<jupyter_code>plt.plot(np.diff(t))
Ts = np.mean(np.diff(t))
Ts, 1/Ts<jupyter_output><empty_output><jupyter_text>Find the Vibrational Frequencies<jupyter_code>def plot_fft(t, y):
    n = len(y)
    Y = (2/n)*np.fft.fft(y)
    n_half = int(np.ceil(n/2))
    freq = np.fft.fftfreq(n, Ts)
    plt.plot(freq[:n_half], np.abs(Y[:n_half]))
    plt.xlabel('Hz')
    plt.grid()
plot_fft(t, accx)
test_sig = np.sin(2*np.pi*20*t) + np.sin(2*np.pi*10*t)
plot_fft(t, test_sig)<jupyter_output><empty_output><jupyter_text>Attitude Rates Control Response<jupyter_code>def plot_rate_response(log, name):
    attitude = log.get_dataset('vehicle_attitude')
    rates_setpoint = log.get_dataset('vehicle_rates_setpoint')
    t = attitude.data['timestamp']/1e6
    y = attitude.data[name + 'speed']
    t_u = rates_setpoint.data['timestamp']/1e6
    u = rates_setpoint.data[name]
    plt.plot(t, np.rad2deg(y))
    plt.plot(t_u, np.rad2deg(u))
    plt.xlabel('t, sec')
    plt.ylabel(name + 'speed deg/s')
    plt.grid()

plt.figure(figsize=(10, 5))
plot_rate_response(log, 'pitch')
plt.gca().set_xlim([244, 253]);
plt.gca().set_ylim([-90, 60])
import numpy as np
import pylab as pl
import sysid

def find_plant(log, name, t0, t1, dt):
    attitude = log.get_dataset('vehicle_attitude')
    rates_setpoint = log.get_dataset('vehicle_rates_setpoint')
    t = attitude.data['timestamp']/1e6
    y = attitude.data[name + 'speed']
    t_u = rates_setpoint.data['timestamp']/1e6
    u = rates_setpoint.data[name]
    y_r = []
    u_r = []
    for ti in np.arange(t0, t1, dt):
        y_r.append(np.interp(ti, t, y))
        u_r.append(np.interp(ti, t_u, u))
    y_r = np.array(y_r)
    u_r = np.array(u_r)
    
    sys = sysid.subspace_det_algo1(y_r, u_r, 6, 3, 0.001, 0.01)
    return sys
a = np.zeros(100)
a[50] = 1
a[65] = 1
np.where(a  > 0.3)
sys = find_plant(log, 'roll', 245, 252, 0.1)
sys
import control
sys = control.ss(sys.A, sys.B, sys.C, sys.D, sys.dt)
sys
G_d = control.ss2tf(sys)
G_d<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1319581802769144, 'fraction_numerical': 0.03729867194122634, 'mean_word_length': 5.310160427807487, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 24, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '350026', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1549, 'n_tokens_neox': 1445, 'n_words': 256, 'path': '03.DroneSystemID.ipynb', 'repo_name': 'jgoppert/aae364_notebook', 'size': '3539'}