fixed PSD averaging bug here too...

88b655e1 · Mr. SQUID · 7a2044f2 · 88b655e1 · 88b655e1
Commit 88b655e1 authored 4 years ago by Mr. SQUID
--- a/SQUID measurement software.ipynb
+++ b/SQUID measurement software.ipynb
@@ -1164,7 +1164,7 @@
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
-    "hide_input": true,
+    "hide_input": false,
    "scrolled": false
   },
   "outputs": [
@@ -1294,8 +1294,8 @@
    "            I_psd = I_psd_trace\n",
    "            V_psd = V_psd_trace\n",
    "        elif n_filled < nmax:\n",
-    "            I_psd = I_psd * n_filled / (n_filled+1) + I_psd_trace\n",
-    "            V_psd = V_psd * n_filled / (n_filled+1) + V_psd_trace\n",
+    "            I_psd = (I_psd * n_filled + I_psd_trace)  / (n_filled+1) \n",
+    "            V_psd = V_psd * n_filled + V_psd_trace) / (n_filled+1) \n",
    "        else:\n",
    "            # this should work because of numpy array index wrapping\n",
    "            I_psd += I_psd_trace / nmax - I_psd_array[ind_next-nmax]/nmax\n",
 %% Cell type:markdown id: tags:
  
 # Load libraries + setup functions
  
 %% Cell type:code id: tags:
  
 ``` python
 import warnings
 warnings.filterwarnings('ignore')
  
 import sys
 import os
 import shutil
 from pyrpl import Pyrpl
 import numpy as np
 import time
 import matplotlib.pyplot as plts
  
 import IPython
 import ipywidgets as widgets
  
 from bokeh.plotting import figure, show
 from bokeh.io import output_notebook, push_notebook
 from bokeh.layouts import column, row
 from bokeh.models import ColumnDataSource, Toggle, Range1d
  
 from scipy.signal import welch
  
 from datetime import datetime
  
 #Arduino Serial connection
 import serial
  
 # Red Pitaya Stuff
  
 HOSTNAME = "rp-f06897.local"
 p = Pyrpl(hostname=HOSTNAME, config='test',gui=False)
 r = p.rp
 s = r.scope
  
 # Some good default settings
 s.input1 = 'in1'
 s.input2 = 'in2'
 s.decimation = 1024 # or s.duration =0.01
 s.average = True
 s.trigger_source = 'immediately'
  
 # For the gui updates during the live view loops
 ipython = IPython.get_ipython()
  
 # For Bokeh
 output_notebook()
  
 # The full scale ranges
 I_full_scale = 520 #uA
 V_full_scale = 800 #uV
  
 # An autoscaling routine
 def auto_scale(r, d):
    d_max = np.max(d)
    d_min = np.min(d)
    d_range = d_max - d_min
    if d_range < 0.01:
        d_range = 0.01
    r.start = d_min - d_range*0.1
    r.end = d_max + d_range*0.1
  
 def get_IV_t(decimation = 1024):
    s.decimation = decimation
    s._start_acquisition()
    time.sleep(s.duration)
    I,V = s._get_curve()
    I *= 20 # we have RP in HV mode
    V *= 20 # we have RP in HV mode
    I /= 0.91 # scaling correction for my RP
    V /= 0.91 # scaling correction for my RP
    I /= 1e4 # Conversion factor Mr SQUID 10000 A/V
    V /= 1e4 # Conversion factor Mr SQuID 10000 V/V
    t = s.times
    return I*1e6,V*1e6,t
 ```
  
 %% Output
  
    INFO:pyrpl:All your PyRPL settings will be saved to the config file
        C:\Users\LocalAdmin\pyrpl_user_dir\config\test.yml
    If you would like to restart PyRPL with these settings, type "pyrpl.exe test" in a windows terminal or
        from pyrpl import Pyrpl
        p = Pyrpl('test')
    in a python terminal.
    INFO:pyrpl.redpitaya:Successfully connected to Redpitaya with hostname rp-f06897.local.
  
  
  
 %% Cell type:markdown id: tags:
  
 # Your information
  
 Please update the practicum round
  
 %% Cell type:code id: tags:
  
 ``` python
 academic_year = "2020_2021"
 practicum_round = 0
 ```
  
 %% Cell type:markdown id: tags:
  
 Run this code to create the folder where your data will be saved and make a subfolder with a copy of the latest analysis scripts.
  
 %% Cell type:code id: tags:
  
 ``` python
 # Make a folder for this group based on the "round number"
 folder = academic_year + "_round_%02d/" % practicum_round
 if not os.path.exists(folder):
    os.makedirs(folder)
  
 gitfolder = "../../Documents/GitHub/rp-squid-python-code/"
  
 # Copy the latest copy of the example analysis code to a subfolder
 if not os.path.exists(folder+"analysis_code/"):
    shutil.copytree(gitfolder + "/analysis_code/", folder+"analysis_code/")
  
 # And put these in the top level folder
 for file in "autoplot.py", "Plot all my data.ipynb":
    if not os.path.exists(folder+file):
        shutil.copy2(gitfolder + "/analysis_code/" + file , folder+file)
 ```
  
 %% Cell type:markdown id: tags:
  
 # IV vs Time Live Display
  
 %% Cell type:code id: tags:
  
 ``` python
 I = []
 V = []
 t = []
  
 def update_display():
    global I,V,t
    I,V,t = get_IV_t()
    I_avg = np.average(I)
    I_ptp = np.max(I) - np.min(I)
    p1.title.text = "Current          Average = %.2f uA     Peak-to-Peak %.2f uA" % (I_avg, I_ptp)
    V_avg = np.average(V)
    V_ptp = np.max(V) - np.min(V)
    p2.title.text = "Voltage          Average = %.2f uV     Peak-to-Peak %.2f uV" % (V_avg, V_ptp)
    R_ptp.value = "%.2f" % (V_ptp / I_ptp)
    R_avg.value = "%.2f" % (V_avg / I_avg)
    source1.data = dict(x=t, y=I)
    source2.data = dict(x=t, y=V)
    if scale_mode.value == "Autoscale":# and False:
        auto_scale(p1.y_range,I)
        auto_scale(p2.y_range,V)
    elif scale_mode.value == "Full range":
        p1.y_range.start = -I_full_scale
        p1.y_range.end = I_full_scale
        p2.y_range.start = -V_full_scale
        p2.y_range.end = V_full_scale
  
 range1 = Range1d()
 range2 = Range1d()
 source1 = ColumnDataSource()
 source2 = ColumnDataSource()
  
 p1 = figure(title="Current", plot_height=300, plot_width=900,toolbar_location=None, y_range = range1)
 #p1 = figure(title="Current", plot_height=300, plot_width=900,toolbar_location=None)
 p1.yaxis.axis_label = 'Curret (uA)'
  
 p2 = figure(title="Voltage", plot_height=300, plot_width=900,toolbar_location=None, y_range = range2)
 #p2 = figure(title="Voltage", plot_height=300, plot_width=900,toolbar_location=None)
  
 p2.xaxis.axis_label = 'Time (s)'
 p2.yaxis.axis_label = 'Voltage (uV)'
  
 p1.line('x', 'y', source=source1)
 p2.line('x', 'y', source=source2)
  
 style = {'description_width': 'initial'}
  
 stop_button = widgets.ToggleButton(description='Stop')
 pause_button = widgets.ToggleButton(description='Pause')
  
 def save_data(w):
    fmt = f"mrsquid_{mrsquid_mode.value}_%Y-%m-%d-%H_%M_%S.dat"
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([t,I,V]).T)
    filename.value = outname
  
 save_button = widgets.Button(description='Save data')
 save_button.on_click(save_data)
 filename = widgets.HTML(description="Last filename: ",style=style)
 filename.value = "(none)"
  
 scale_mode = widgets.RadioButtons(options=['Autoscale', 'Full range', 'Freeze range'])
 mrsquid_mode = widgets.RadioButtons(options=['V-I', 'V-Phi'])
  
 R_ptp = widgets.Text(description="V/I ptp (ohms)",style=style)
 R_avg = widgets.Text(description="V/I avg (ohms)",style=style)
 R_ptp.layout.width = '150px'
 R_avg.layout.width = '150px'
  
 update_display()
 target = show(column(p1,p2), notebook_handle=True)
 rows = []
 rows.append(widgets.HBox([stop_button,pause_button,scale_mode, mrsquid_mode]))
 rows.append(widgets.HBox([save_button, filename]))
 rows.append(widgets.HBox([R_ptp, R_avg]))
 display(widgets.VBox(rows))
  
 while True:
    ipython.kernel.do_one_iteration()
    if not pause_button.value:
        update_display()
    ipython.kernel.do_one_iteration()
    if stop_button.value:
        break
    push_notebook(handle=target)
  
 print("Live view done")
 ```
  
 %% Output
  
  
  
  
    Live view done
  
 %% Cell type:markdown id: tags:
  
 # IV Live Display
  
 %% Cell type:code id: tags:
  
 ``` python
 I = []
 V = []
 t = []
  
 def update_display():
    global I,V,t
    ipython.kernel.do_one_iteration()
    I,V,t = get_IV_t()
    ipython.kernel.do_one_iteration()
    source1.data = dict(x=I, y=V)
  
    # Update widgets
    I_avg = np.average(I)
    I_ptp = np.max(I) - np.min(I)
    V_avg = np.average(V)
    V_ptp = np.max(V) - np.min(V)
    widgets = (I_avg_w, V_avg_w,  R_avg_w, I_ptp_w, V_ptp_w, R_ptp_w)
    vals = (I_avg, V_avg, V_avg / I_avg, I_ptp, V_ptp, V_ptp / I_ptp)
    for w,v in zip(widgets,vals):
        w.value = "%.2f" % v
  
    # Update scaling
    if scale_mode.value == "Autoscale":
        if plot_ref_button.value:
            auto_scale(p1.x_range,np.concatenate(I,np.array(source3.data['x'])))
            auto_scale(p1.y_range,np.concatenate(V,np.array(source3.data['y']))
        else:
            auto_scale(p1.x_range,I)
            auto_scale(p1.y_range,V)
    elif scale_mode.value == "Full range":
        p1.x_range.start = -x_full_scale
        p1.x_range.end = x_full_scale
        p1.y_range.start = -y_full_scale
        p1.y_range.end = y_full_scale
  
    x = np.array([np.min(I), np.max(I)])
    R = line_slope.value
    y = np.array([V_avg - I_ptp*R/2, V_avg+I_ptp*R/2])
    source2.data = dict(x=x, y=y)
  
    if plot_line_button.value:
        slope_line.visible = True
    else:
        slope_line.visible = False
  
    if plot_ref_button.value:
        ref_trace.visible = True
    else:
        ref_trace.visible = False
  
 range1 = Range1d()
 range2 = Range1d()
 source1 = ColumnDataSource()
 source2 = ColumnDataSource()
 source3 = ColumnDataSource()
  
 p1 = figure(title="IV", plot_height=500, plot_width=600,toolbar_location=None, x_range = range1, y_range = range2)
 p1.yaxis.axis_label = 'Measured Voltage (uV)'
 p1.xaxis.axis_label = 'Applied Current (uA)'
 ref_trace = p1.line('x','y',source=source3, color='magenta')
 ref_trace.visible = False
 p1.line('x', 'y', source=source1)
 slope_line = p1.line('x','y',source=source2, color='red')
 slope_line.visible = False
  
  
 x_full_scale = I_full_scale #uA
 y_full_scale = V_full_scale # uV
  
 stop_button = widgets.ToggleButton(description='Stop')
 pause_button = widgets.ToggleButton(description='Pause')
  
 scale_mode = widgets.RadioButtons(options=['Autoscale', 'Full range', 'Freeze range'])
 mrsquid_mode = widgets.RadioButtons(options=['V-I', 'V-Phi'])
  
 style = {'description_width': 'initial'}
  
 I_avg_w = widgets.Text(description="I_avg (uA)",style=style)
 V_avg_w = widgets.Text(description="V_avg (uV)",style=style)
 I_ptp_w = widgets.Text(description="I_ptp (uA)",style=style)
 V_ptp_w = widgets.Text(description="V_ptp (uV)",style=style)
 R_ptp_w = widgets.Text(description="V/I ptp (ohms)",style=style)
 R_avg_w  = widgets.Text(description="V/I avg (ohms)",style=style)
  
 plot_line_button = widgets.ToggleButton(description="Plot line")
 line_slope = widgets.FloatText(description="R_line (ohms)")
 line_slope.value = 50
 fit_slope_button = widgets.Button(description="Fit")
  
 def fit_slope(w):
    coef = np.polyfit(I,V,1)
    line_slope.value = coef[0]
 fit_slope_button.on_click(fit_slope)
  
 plot_ref_button = widgets.ToggleButton(description="Plot ref")
 take_ref_button = widgets.Button(description="Take as ref")
  
 def take_ref(w):
    global source3, I, V
    source3.data = dict(x=I, y=V)
 take_ref_button.on_click(take_ref)
  
  
 for w in (I_avg_w, V_avg_w, I_ptp_w, V_ptp_w, R_ptp_w, R_avg_w):
    w.layout.width = '150px'
  
 def save_data(w):
    fmt = f"mrsquid_{mrsquid_mode.value}_%Y-%m-%d-%H_%M_%S.dat"
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([t,I,V]).T)
    filename.value = outname
  
 save_button = widgets.Button(description='Save data')
 save_button.on_click(save_data)
 filename = widgets.HTML(description="Last filename: ",style=style)
 filename.value = "(none)"
  
 take_ref([])
 update_display()
 target = show(p1, notebook_handle=True)
 rows = []
  
 rows.append(widgets.HBox([I_avg_w, V_avg_w, R_avg_w]))
 rows.append(widgets.HBox([I_ptp_w, V_ptp_w, R_ptp_w]))
 rows.append(widgets.HBox([plot_line_button, line_slope, fit_slope_button]))
 rows.append(widgets.HBox([plot_ref_button, take_ref_button]))
 rows.append(widgets.HBox([save_button, filename]))
 rows.append(widgets.HBox([stop_button,pause_button,scale_mode,mrsquid_mode]))
 controls = widgets.VBox(rows)
 display(controls)
  
 while True:
    ipython.kernel.do_one_iteration()
    if not pause_button.value:
        update_display()
    if stop_button.value:
        break
    push_notebook(handle=target)
  
 print("Live view stopped")
 ```
  
 %% Output
  
  
  
  
    ---------------------------------------------------------------------------
    TypeError                                 Traceback (most recent call last)
    <ipython-input-45-67091be2a7a5> in <module>
        133     ipython.kernel.do_one_iteration()
        134     if not pause_button.value:
    --> 135         update_display()
        136     if stop_button.value:
        137         break
    <ipython-input-45-67091be2a7a5> in update_display()
         23     if scale_mode.value == "Autoscale":
         24         if plot_ref_button.value:
    ---> 25             auto_scale(p1.x_range,np.concatenate(I,source3.data['x']))
         26             auto_scale(p1.y_range,np.concatenate(V,source3.data['y']))
         27         else:
    <__array_function__ internals> in concatenate(*args, **kwargs)
    TypeError: 'list' object cannot be interpreted as an integer
  
 %% Cell type:markdown id: tags:
  
 # Temperature trace recorder
  
 %% Cell type:code id: tags:
  
 ``` python
 # Global variables
 R = []
 T = []
 tm = []
 Iavg = []
 Iptp = []
 Vavg = []
 Vptp = []
 t0 = time.time()
 timer_t0 = time.time()
  
 # The functions
  
 def update_iv():
    I,V,_ = get_IV_t()
    source1.data = dict(x=I, y=V)
    p1.title.text = "IV    Iptp = %.2f uA" % (np.max(I) - np.min(I))
  
  
 def get_temperature(serial_port, baudrate, timeout):
    # serial_port='COM7', baudrate=9600, timeout=1
  
    #variables (could be adjusted to inputs?)
    alpha = 4.1139e-3; R0 = 1000; R1=1081; Vin=3.272        #should be according to literature: alpha = 3.85e-3 & R0 = 1000
  
    #open Serial port to Arduino + flush what's still on it
    ser = serial.Serial(serial_port,baudrate,timeout=0.5)
    ser.flush()
  
    #sent command to Arduino so we can read it out:
    time.sleep(0.5)                       #needs to sleep, otherwise we get a: 'can not convert float... at the Vout-line'
    ser.reset_input_buffer()
    ser.write(b'g')
  
    #Serial read + Temperature calculations:
    Vout = float(str(ser.readline().decode("utf-8")))*(Vin/1023.0)          #calculates the voltage from the voltage divider (over the PT1000)
    Rpt = (Vout*R1)/(Vin-Vout)                                              #calculates the resistance of the PT1000
    Temp_C = ((Rpt/R0)-1)/alpha                                             #with the above calculated resistance and 'known' variable we calculate the temperature in Celsius
    Temp_K = Temp_C + 273.15                                                #Temperature in Kelvin
  
    return Temp_K
  
  
 def take_measurement():
    global timer_t0
    timer_t0 = time.time()
    T.append(get_temperature('COM7',9600,1))
    I,V,_ = get_IV_t()
    R.append(np.polyfit(I,V,1)[0])
    Iavg.append(np.average(I))
    Iptp.append(np.max(I)-np.min(I))
    Vavg.append(np.average(V))
    Vptp.append(np.max(V)-np.min(V))
    tm.append(time.time()-t0)
    source2.data = dict(x=tm,y=R)
    source3.data = dict(x=tm,y=T)
    source4.data = dict(x=T,y=R)
    p2.title.text = "R vs time   R = %.2f Ohms" % R[-1]
    p3.title.text = "T vs time   T = %.2f K" % T[-1]
  
 # The plots
  
 p1 = figure(title="IV", plot_height=300, plot_width=300,toolbar_location=None)
 p1.yaxis.axis_label = 'Measured Voltage (uV)'
 p1.xaxis.axis_label = 'Applied Current (uA)'
 p2 = figure(title="R vs time", plot_height=300, plot_width=300,toolbar_location=None)
 p2.yaxis.axis_label = 'Resistance (Ohm)'
 p2.xaxis.axis_label = 'Time (s)'
 p3 = figure(title="T vs time", plot_height=300, plot_width=300,toolbar_location=None)
 p3.yaxis.axis_label = 'Temperature (K)'
 p3.xaxis.axis_label = 'Time (s)'
 p4 = figure(title="R vs T", plot_height=300, plot_width=900,toolbar_location=None)
 p4.yaxis.axis_label = 'Resistance (Ohm)'
 p4.xaxis.axis_label = 'Temperature (K)'
  
 init_data = dict(x=[0,1],y=[0,1])
  
 source1 = ColumnDataSource()
 source1.data = init_data
 p1.line('x', 'y', source=source1)
 source2 = ColumnDataSource()
 source2.data = init_data
 p2.line('x', 'y', source=source2)
 p2.circle('x', 'y', source=source2, fill_color="white", size=8)
 source3 = ColumnDataSource()
 source3.data = init_data
 p3.line('x', 'y', source=source3)
 p3.circle('x', 'y', source=source3, fill_color="white", size=8)
 source4 = ColumnDataSource()
 source4.data = init_data
 p4.circle('x', 'y', source=source4, fill_color="white", size=8)
  
 take_measurement()
  
 #target = show(column(row(p1,p2,p3),p4), notebook_handle=True)
 target = show(row(p1,p2,p3), notebook_handle=True)
  
 # The widgets
  
 style = {'description_width': 'initial'}
  
 save_button = widgets.Button(description='Save data')
 def save_data(w):
    fmt = f"mrsquid_R_vs_T_%Y-%m-%d-%H_%M_%S.dat"
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([R,T,tm,Iavg,Iptp,Vavg,Vptp]).T)
    filename.value = outname
 save_button.on_click(save_data)
  
 filename = widgets.HTML(description="Last filename: ",style=style)
 filename.value = "(none)"
  
 reset_button = widgets.Button(description='Restart data')
 def reset(w):
    global R,T,tm,Iavg,Iptp,Vavg,Vptp
    R = []
    T = []
    tm = []
    Iavg = []
    Iptp = []
    Vavg = []
    Vptp = []
    t0 = time.time()
    take_measurement()
 reset_button.on_click(reset)
  
 meas_button = widgets.Button(description='Measure now')
 def meas(w):
    take_measurement()
 meas_button.on_click(meas)
  
 stop_button = widgets.ToggleButton(description='Stop')
 pause_button = widgets.ToggleButton(description='Pause')
 interval = widgets.FloatText(description="Measurement delay (s)",style={'description_width': 'initial'})
 interval.value = 0
  
 rows = []
 rows.append(widgets.HBox([interval,meas_button]))
 rows.append(widgets.HBox([pause_button, save_button, reset_button, stop_button]))
 rows.append(widgets.HBox([filename]))
 controls = widgets.VBox(rows)
 display(controls)
  
 while True:
    ipython.kernel.do_one_iteration()
    update_iv()
    if (time.time()-timer_t0) > interval.value and not pause_button.value:
        take_measurement()
    if stop_button.value:
        break
    push_notebook(handle=target)
    if (time.time() - t0) > 3600:
        save_data([])
        break
  
 print("Live view stopped")
 ```
  
 %% Output
  
  
  
  
    Live view stopped
  
 %% Cell type:markdown id: tags:
  
 # Power Spectral Density Live Display
  
 %% Cell type:code id: tags:
  
 ``` python
 # The traces we download
 I = []
 V = []
 t = []
  
 # The calculated spectral densities
 I_psd = []
 V_psd = []
 f = []
  
 I_avg = 0
 I_ptp = 0
 V_avg = 0
 V_ptp = 0
  
 # Each PSD trace is 8193 points long
 n_psd = 8193
 nmax  =  100
 I_psd_array = np.zeros([nmax,n_psd])
 V_psd_array = np.zeros([nmax,n_psd])
  
 # For keeping track of how many we've filled and where the next goes
 ind_next = 0
 n_filled = 0
  
 def update_display():
    global I,V,t,f,I_psd,V_psd,last_trace_time,ind_next,n_filled
    global I_avg, I_ptp, V_avg, V_ptp
    I,V,t = get_IV_t(decimation=decimation.value)
    I_avg = np.average(I)
    I_ptp = np.max(I) - np.min(I)
    p1.title.text = "Current          Average = %.2f uA     Peak-to-Peak %.2f uA" % (I_avg, I_ptp)
    V_avg = np.average(V)
    V_ptp = np.max(V) - np.min(V)
    p2.title.text = "Voltage          Average = %.2f uV     Peak-to-Peak %.2f uV" % (V_avg, V_ptp)
    f, I_psd_trace = welch(I,1/t[1],nperseg=len(I))
    f, V_psd_trace = welch(V,1/t[1],nperseg=len(V))
  
    if enable_average.value:
        # ind_next: the running index
        # n_filled: the number of traces taken since the last reset
  
        I_psd_array[ind_next,:] = I_psd_trace
        V_psd_array[ind_next,:] = V_psd_trace
  
        if n_filled  == 0:
            I_psd = I_psd_trace
            V_psd = V_psd_trace
        elif n_filled < nmax:
-            I_psd = I_psd * n_filled / (n_filled+1) + I_psd_trace
-            V_psd = V_psd * n_filled / (n_filled+1) + V_psd_trace
+            I_psd = (I_psd * n_filled + I_psd_trace)  / (n_filled+1)
+            V_psd = V_psd * n_filled + V_psd_trace) / (n_filled+1)
        else:
            # this should work because of numpy array index wrapping
            I_psd += I_psd_trace / nmax - I_psd_array[ind_next-nmax]/nmax
            V_psd += V_psd_trace / nmax - V_psd_array[ind_next-nmax]/nmax
        ind_next += 1
        if ind_next == 100:
            ind_next = 0
        if n_filled < 100:
            n_filled += 1
            averaged_num.value = "%d" % n_filled
  
    else:
        I_psd = I_psd_trace
        V_psd = V_psd_trace
        n_filled = 0
        ind_next = 0
        if averaged_num.value != "1":
            averaged_num.value = "1"
  
    source1.data = dict(x=f, y=I_psd)
    source2.data = dict(x=f, y=V_psd)
    if manual_freq.value:
        p1.x_range.start = freq_start.value
        p1.x_range.end = freq_stop.value
        p2.x_range.start = freq_start.value
        p2.x_range.end = freq_stop.value
    else:
        p1.x_range.start = f[0]
        p1.x_range.end = f[-1]
        p2.x_range.start = f[0]
        p2.x_range.end = f[-1]
  
 range1 = Range1d()
 range2 = Range1d()
 source1 = ColumnDataSource()
 source2 = ColumnDataSource()
  
 p1 = figure(title="Current", plot_height=300, y_axis_type="log",
            plot_width=900, toolbar_location=None, x_range = range1)
 p1.yaxis.axis_label = 'Curret PSD dB(A^2/Hz)'
  
 p2 = figure(title="Voltage", plot_height=300, y_axis_type="log",
            plot_width=900, toolbar_location=None, x_range = range2)
 p2.xaxis.axis_label = 'Frequency (Hz)'
 p2.yaxis.axis_label = 'Voltage PSD dB(V^2/Hz)'
  
 p1.line('x', 'y', source=source1)
 p2.line('x', 'y', source=source2)
  
 stop_button = widgets.ToggleButton(description='Stop')
 pause_button = widgets.ToggleButton(description='Pause')
  
 def reset_average(w):
    global n_filled, ind_next
    n_filled = 0
    ind_next = 0
  
 # For the decimation GUI
 options = []
 for d,t in zip(s.decimation_options[9:], s.duration_options[9:]):
    if t < 1e-3:
        t*= 1e6
        u = "us"
    elif t<1:
        t*= 1e3
        u = "ms"
    else:
        u = "s"
    options.append((str(d) + ", " "{0:.3g}".format(t) + " " + u,d))
 options
  
 def update_freq(w):
    s.decimation = decimation.value
    trace_info.value = "Frequency resolution: %.2f Hz &nbsp; &nbsp; &nbsp; &nbsp;" % (1/s.duration)
    trace_info.value += "Max frequency: %.2f Hz" % (1/s.duration*n_psd)
  
 decimation = widgets.Dropdown(description="Decimation:",
                              options=options)
 decimation.observe(reset_average)
 decimation.observe(update_freq)
 trace_info = widgets.HTML()
  
 enable_average = widgets.Checkbox(description="Enable averaging (max 100 traces)", style=style)
 averaged_num = widgets.HTML(value="1", description="Number averaged:", style=style)
 reset_button = widgets.Button(description="Reset averaging")
 reset_button.on_click(reset_average)
  
 manual_freq = widgets.Checkbox(description="Manual frequency range")
 freq_start = widgets.FloatText(value=0,description="Start (Hz)")
 freq_stop = widgets.FloatText(value=1e5,description="Stop (Hz)")
  
 def save_data(w):
    fmt = f"mrsquid_PSD_%Y-%m-%d-%H_%M_%S.dat"
    header = "I_avg %e I_ptp %e V_avg %e V_ptp %e" % (I_avg, I_ptp, V_avg, V_ptp )
    outname = folder + "/" + datetime.now().strftime(fmt)
    np.savetxt(outname, np.array([f,I_psd,V_psd]).T, header=header)
    filename.value = outname
  
 save_button = widgets.Button(description='Save data')
 save_button.on_click(save_data)
 filename = widgets.HTML(description="Last filename: ",style=style)
 filename.value = "(none)"
  
 update_display()
 update_freq(0)
 target = show(column(p1,p2), notebook_handle=True)
 rows = []
 rows.append(widgets.HBox([decimation,trace_info]))
 rows.append(widgets.HBox([enable_average, averaged_num, reset_button]))
 rows.append(widgets.HBox([manual_freq,freq_start,freq_stop]))
 rows.append(widgets.HBox([save_button, filename]))
 rows.append(widgets.HBox([stop_button,pause_button]))
 display(widgets.VBox(rows))
  
 while True:
    ipython.kernel.do_one_iteration()
    if not pause_button.value:
        update_display()
    ipython.kernel.do_one_iteration()
    if stop_button.value:
        break
    push_notebook(handle=target)
    time.sleep(0.1)
  
 print("Live view done")
 ```
  
 %% Output
  
  
  
  
    Live view done
--- a/SQUID measurement software.md
+++ b/SQUID measurement software.md
@@ -527,7 +527,7 @@ print("Live view stopped")
 # Power Spectral Density Live Display
 <!-- #endregion -->

-```python hide_input=true
+```python hide_input=false
 # The traces we download
 I = []
 V = []
@@ -577,8 +577,8 @@ def update_display():
            I_psd = I_psd_trace
            V_psd = V_psd_trace
        elif n_filled < nmax:
-            I_psd = I_psd * n_filled / (n_filled+1) + I_psd_trace
-            V_psd = V_psd * n_filled / (n_filled+1) + V_psd_trace
+            I_psd = (I_psd * n_filled + I_psd_trace)  / (n_filled+1) 
+            V_psd = V_psd * n_filled + V_psd_trace) / (n_filled+1) 
        else:
            # this should work because of numpy array index wrapping
            I_psd += I_psd_trace / nmax - I_psd_array[ind_next-nmax]/nmax