text: # ---
  # jupyter:
  #   jupytext:
  #     formats: ipynb,py
  #     text_representation:
  #       extension: .py
  #       format_name: light
  #       format_version: '1.5'
  #       jupytext_version: 1.13.7
  #   kernelspec:
  #     display_name: Python 3 (ipykernel)
  #     language: python
  #     name: python3
  # ---
  
  # +
  # setup
  from bokeh.plotting import figure, show, output_notebook
  from bokeh.models import (
      ColumnDataSource, Range1d, DataRange1d, DatetimeAxis,
      TickFormatter, DatetimeTickFormatter, FuncTickFormatter,
      Grid, Legend, Plot, BoxAnnotation, Span, CustomJS, Rect, Circle, Line,
      HoverTool, BoxZoomTool, PanTool, WheelZoomTool,
      WMTSTileSource, LabelSet
      )
  from bokeh.models.widgets import DataTable, DateFormatter, TableColumn, Div
  
  import numpy as np
  import sys
  import os
  
  sys.path.append(os.path.join(os.getcwd(), 'plot_app'))
  from pyulog import *
  from pyulog.px4 import *
  from app.plot_app.plotting import *
  from app.plot_app.config import *
  from app.notebook_helper import *
  
  output_notebook()
  
  # +
  file_name = '/home/markagate/px4Analysis/00045-2022-03-10_21-16-58.ulg' # TODO: fill in file name
  ulog = ULog(file_name)
  data = ulog.data_list
  px4_ulog = PX4ULog(ulog)
  px4_ulog.add_roll_pitch_yaw()
  use_downsample = False
  
  print("message names: {:}".format(sorted([d.name for d in data])))
  print_ulog_info(ulog)
  # -
  
  
  
  use_downsample = True # you may want to activate this for large logs (But you will not see all samples when zooming in)
  
  # load a dataset to inspect the field names and types
  sensor_data = [ elem for elem in data if elem.name == 'vehicle_local_position' and elem.multi_id == 0][0]
  types_list = [(f.type_str, f.field_name) for f in sensor_data.field_data]
  for ftype, fname in types_list: print("{:10s} {:}".format(ftype, fname))
  
  plot = plot_map(ulog, plot_config) # gps map
  if plot != None: show(plot)
  
  # plot raw acceleration sensor data
  data_plot = DataPlot(data, plot_config, 'sensor_combined',
          y_axis_label = '[m/s^2]', title = 'Raw Acceleration', plot_height = 'normal')
  data_plot.add_graph(['accelerometer_m_s2[0]', 'accelerometer_m_s2[1]', 'accelerometer_m_s2[2]'],
                      colors3, ['x', 'y', 'z'], use_downsample=use_downsample)
  show(data_plot.finalize())
  
  # +
  # a more complex plot with multiple datasets
  data_plot = DataPlot(data, plot_config, 'vehicle_magnetometer', y_start=0, title = 'Thrust and Magnetic Field',
                       plot_height='normal')
  data_plot.add_graph([lambda data: ('len_mag', 
      np.sqrt(data['magnetometer_ga[0]']**2 + data['magnetometer_ga[1]']**2 + data['magnetometer_ga[2]']**2))],
      colors2[0:1], ['Norm of Magnetic Field'], use_downsample=use_downsample)
  
  data_plot.change_dataset('actuator_controls_0')
  data_plot.add_graph([lambda data: ('thrust', data['control[3]'])], colors2[1:2], ['Thrust'],
                      use_downsample=use_downsample)
  
  show(data_plot.finalize())
  # -
  
  
  
  # +
  # plot low-pass filtered raw acceleration sensor data
  from scipy.signal import butter, lfilter
  cur_data = ulog.get_dataset('sensor_combined').data
  t = cur_data['timestamp']
  
  fs = len(t) / ((t[-1]-t[0])/1e6) # sample rate [Hz]
  cutoff = 10 # cutoff freq [Hz]
  order = 5
  B, A = butter(order, cutoff / (fs / 2), btype='low') # Butterworth low-pass
  filtered_signal = lfilter(B, A, cur_data['accelerometer_m_s2[0]'])
  
  p = figure(plot_width = 800, active_scroll='wheel_zoom')
  p.line(t, cur_data['accelerometer_m_s2[0]'], color='red', alpha = 0.5)
  p.line(t, filtered_signal, color='blue', alpha = 0.8)
  show(p)
  # -
  
  
  
  # get the raw acceleration
  sensor_combined = ulog.get_dataset('sensor_combined').data
  ax = sensor_combined['accelerometer_m_s2[0]']
  ay = sensor_combined['accelerometer_m_s2[1]']
  az = sensor_combined['accelerometer_m_s2[2]']
  t = sensor_combined['timestamp']
  
  # and plot it
  p = figure(plot_width = 800, active_scroll='wheel_zoom')
  p.line(t, ax, color='red', alpha = 0.8, legend="accel x")
  p.line(t, ay, color='green', alpha = 0.8, legend="accel y")
  p.line(t, az, color='blue', alpha = 0.8, legend="accel z")
  show(p)
  
  # select a window
  dt = (t[-1]-t[0]) / len(t) / 1e6 # delta t in seconds
  start_index = int(20 / dt) # select start of window (seconds): make sure the vehicle is hovering at that point
  window_len_s = 3 # window length in seconds
  tw=t[start_index:int(start_index+window_len_s/dt)]
  axw=ax[start_index:int(start_index+window_len_s/dt)]
  ayw=ay[start_index:int(start_index+window_len_s/dt)]
  azw=az[start_index:int(start_index+window_len_s/dt)]
  
  # and plot it
  p = figure(plot_width = 800, active_scroll='wheel_zoom')
  p.line(tw, axw, color='red', alpha = 0.8, legend="accel x")
  p.line(tw, ayw, color='green', alpha = 0.8, legend="accel y")
  p.line(tw, azw, color='blue', alpha = 0.8, legend="accel z")
  show(p)
  
  # +
  # FFT frequency plot
  import scipy
  import scipy.fftpack
  from scipy import pi
  
  FFT_x = abs(scipy.fft(axw))
  FFT_y = abs(scipy.fft(ayw))
  FFT_z = abs(scipy.fft(azw))
  
  freqs = scipy.fftpack.fftfreq(len(axw), dt)
  
  p = figure(plot_width = 800, active_scroll='wheel_zoom')
  p.line(freqs,20*scipy.log10(FFT_x), color='red', alpha = 0.8, legend="x")
  p.line(freqs,20*scipy.log10(FFT_y), color='green', alpha = 0.8, legend="y")
  p.line(freqs,20*scipy.log10(FFT_z), color='blue', alpha = 0.8, legend="z")
  p.legend.click_policy="hide"
  show(p)
  # -
  
  # ### spectogram
  
  # +
  window = 'hann'
  window_length = 256
  noverlap = 128
  
  data_set = {'ax': ax, 'ay': ay, 'az': az}
  legend = ['X','Y''Z']
  
  fs = 250
  
  # calculate the spectogram
  psd = dict()
  for key, val in data_set.items():
      f, t, psd[key] = scipy.signal.spectrogram(
          val,fs=fs, window=window, nperseg=window_length, noverlap=noverlap, scaling='density')
      
  # sum all psd's
  key_it = iter(psd)
  sum_psd = psd[next(key_it)]
  for key in key_it:
      sum_psd += psd[key]
  
  
  # -
  
  # #### plot with bokeh
  
  def plot_spec(t, f, psd, legend):
      color_mapper = LinearColorMapper(palette=viridis(256),low=-80, high=0)
  
      im = [10 * np.log10(psd)]
      p = figure(title='Acceleration Power Spectral Density ' + legend + ' [dB]',
          plot_width=800,x_range=(t[0], t[-1]),y_range=(f[0], f[-1]),
          x_axis_label='Time',y_axis_label='[Hz]',toolbar_location='above')
      p.image(image=im, x=t[0], y=f[0],dw=(t[-1]-t[0]), dh=(f[-1]-f[0]),color_mapper=color_mapper)
      color_bar = ColorBar(color_mapper=color_mapper,
                           major_label_text_font_size="5pt",
                           ticker=BasicTicker(desired_num_ticks=5),
                           formatter=PrintfTickFormatter(format="%f"),
                           label_standoff=6, border_line_color=None, location=(0, 0))
      p.add_layout(color_bar,'right')
      p.add_tools(BoxZoomTool(dimensions="width"))
      return p
  
  
  show(plot_spec(t,f,sum_psd,'X Y Z'))
  
  
  