A Short Guide to Using Python For Data Analysis In Experimental Physics (V3)

A Short Guide to Using Python For Data Analysis In Experimental Physics (V3) | Authorea try { document.documentElement.classList.add('js'); } catch (e) { } var _gaq = _gaq || []; _gaq.push(['_setAccount', 'G-8VDV14Y67G']); _gaq.push(['_trackPageview']); (function() { var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true; ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); })(); Skip to main content Preprints Collections Wiley Open Research IET Open Research Ecological Society of Japan All Collections About About Authorea FAQs Contact Us Quick Search anywhere Search for preprint articles, keywords, etc. Search Search ADVANCED SEARCH SCROLL This is a preprint and has not been peer reviewed. Data may be preliminary. 1 April 2026 V3 Latest version Share on A Short Guide to Using Python For Data Analysis In Experimental Physics (V3) Authors : Nathanael A. Fortune 0000-0002-2704-9969 [email protected] and Rebecca Webster Authors Info & Affiliations https://doi.org/10.22541/au.152961025.52816543/v3 12991 views 486 downloads Contents Abstract Importing and exporting data loading data from a CSV file importing the data file saving data to a CSV file Other file handling methods Plotting data using error bars Curve fitting model fit: V_0 laser intensity model fit: polarizer angle error model fit: errors in calculated quantities putting it all together: estimate of uncertainty fitting the model to data graphing the results Evaluating the goodness of fit chi-square test Data Smoothing and Differentiation Scipy.Signal.Savgol_filter Sample Python code Butterworth filters Differentiation Interpolation and Peak Finding InterpolatedUnivariateSpline peak finding Working with units and uncertainties Calculations with units Calculations with uncertainties Installing Python Installing Anaconda Navigator Building your own distribution Using command line installers conda install pip install Information & Authors Metrics & Citations View Options References Figures Tables Media Share Abstract Common signal processing tasks in the numerical handling of experimental data include interpolation, smoothing, and propagation of uncertainty. A comparison of experimental results to a theoretical model further requires curve fitting, the plotting of functions and data, and a determination of the goodness of fit. These tasks often typically require an interactive, exploratory approach to the data, yet for the results to be reliable, the original data needs to be freely available and resulting analysis readily reproducible. In this article, we provide examples of how to use the Numerical Python (Numpy) and Scientific Python (SciPy) packages and interactive Jupyter Notebooks to accomplish these goals for data stored in a common plain text spreadsheet format. Sample Jupyter notebooks containing the Python code used to carry out these tasks are included and can be used as templates for the analysis of new data. On your own computer The Authorea server will work for all of the examples except the (aptly named) Python packages Pint and Uncertainties used for numerical calculations using units and/or uncertainties, as those packages are not currently installed on that server (but maybe someday?). That said, you will want to write, run, and store your own data and programs on your own system. If you don't have access to a local Jupyter webserver with the packages you need, you will want to install and run Python and Jupyter on your own computer. The good news is that everything needed is available for free (except the computer)! See Section \ref{662059}: Installing Python for instructions on how to do this. Importing and exporting data There are many ways to import, export, and represent data in files. Here we provide just enough to get you started but it might very well be all you need. In these examples, we assume you have first entered the data into a spreadsheet program, then exported that data in 'CSV' (comma separated variable) file format. We use the CSV file format here because it is ubiquitous: all spreadsheets have the ability to export and import data in the CSV file format. Data in other common plain text file formats (such as tab delimited) can also be imported by making a few small modifications to the examples provided below. CSV spreadsheet file format Let's look at a particular CSV data file titled Calibration_650nm.csv . The file consists of a single header row of text which we need to skip over when loading the numerical data, followed by three columns of data, one for each measured variable. Here's what the first few lines look like when opened in a text editor (notice the commas separating the values within each row): angle, V_pd, V_pd_error0.000000000000000000e+00,3.236249999999999716e+01,2.492398249033692462e-021.500000000000000000e+01,3.008079999999999998e+01,1.536232648449061544e-013.000000000000000000e+01,2.402850000000000108e+01,2.527732747539992997e-01 loading data from a CSV file We are now going to use the Numerical Python (numpy) command loadtxt to load data from this CSV text file. Each column of data will become a Numerical Python ( numpy ) array. The process consists of two steps: 1. importing the loadtxt command from Numerical Python 2. using the loadtxt command to transfer the data from the file importing commands loadtxt is not part of the core Python language but is instead part of the Numerical Python package ( numpy ). We must therefore "import" it into our program before using. There are two common methods of doing this. In the first method , each individual function or module (such as loadtxt from numpy ) is imported as needed . We call this the direct numpy import style . This method is used, for example, in the excellent Python-based textbook Computational Physics by Mark Newman \cite{mark2013} . Here is another example: with result In the second method , we first import all of the numpy package ( and, optionally, provide an abbreviation for numpy such as np .) We then add numpy. (or the abbreviation np. ) as a prefix when using a function from the numpy package. We call this the traditional import style . We call this the traditional method because this is what you will find in the examples included in the official numpy user guide and quickstart tutorials . The key to using the traditional method is to be sure to specify that you are using a function or module from the numpy package by preceding it either with numpy or the abbreviation of your choice (such as np ). This is why the loadtxt commands used in the examples below are written np.loadtxt instead of loadtxt . Here is another example: with the same result as before: One one hand, the direct import method has the advantage of being lean, clean, and easy to read. Computer scientists love it. On the other hand, it does not work well if we need to use multiple packages containing identically named functions. For example, the math package math , the complex mathematics package cmath , the numerical python package numpy , and the numerical calculation of uncertainties package uncertainties all have trigonometric functions called sin and cos . To keep clear which function we are using from which package, and when, we will usually use the traditional import style for packages. importing the data file We now show how to use loadtxt to import the csv format spreadsheet file Calibration_650nm.csv . Unless specified otherwise, the data file is assumed to be in the same file folder as the Python program. The loadtxt command: What these commands do : data_file = file_folder + file_name tells Python what the name of the file is and where to find it! Tip: If you have trouble determining how to specify the file_folder location for your data, an easy workaround is to first put the data in the same folder as your Python program, then either (a) set file_folder = '' as in the example above or (b) replace np.loadtxt(file_folder + file_name, ) with np.loadtxt(file_name, ) . delimiter = ',' tells loadtxt that your data is in comma separated variable format (CSV). The character used to separate one column of data from another is called the 'delimiter.' See the numpy.loadtxt manual page for details. skiprows = 1 tells loadtxt to skip over one row of text in the CSV file before looking for data to load into arrays. This is because the first row of text contains names for each column of data instead of data values (as shown in Table \ref{296817}). unpack = True tells loadtxt that the data is in a 'packed' data format (in which each variable corresponds to a different column instead of to a different row )and therefore needs to be 'unpacked' when loaded. This is the typical arrangement for data in spreadsheets. Use unpack = False if the data is in an 'unpacked' data format (in which each variable corresponds to a different row instead of a different column ). usecols = (0, 1, 2) says the data you are looking for is in the first 3 columns, which are numbered 0, 1, and 2 (because Python always starts from zero when counting). In our case, since there are only 3 columns of data and we want to use all three, this command is unnecessary. You could leave it out and everything would work just fine for this particular data file. If, however, you wanted to load data from a file with many columns but only needed data from column number 0, column number 3, and column number 4, you would need to include usecols = (0,3,4) within the loadtxt command. angle, V_pd, V_pd_error = np.loadtxt() tells Numerical Python to create an array called angle and fill it with values from the first column of data in the CSV spreadsheet file, then create an array called V_pd and fill it with values from the second column of data, and finally create an array called V_pd_error and fill it with values from the third column of data. The result is three shiny new numpy data arrays we can use in our calculations. What if your data is in a different text file format (such as tab delimited)? In that case you can still use loadtxt to import your data as long as you modify the delimiter command to match the file format. See the numpy.loadtxt manual page for details. saving data to a CSV file For completeness, we now provide an example of how to use the savetxt command from numpy to save data in csv format. There are many ways to save files in spreadsheet format but this is the simplest method we've found so far using numpy . Direct numpy import style: Traditional numpy import style: Why do we need the line data = array([angle, V_pd]).T ? We need it because ordinarily savetxt would save the data in what Python calls 'unpacked' format, a format in which each variable corresponds to a different row instead of to a different column . This is often convenient but is not what we wanted in this particular case. We therefore did the following clever trick before saving the data to a file: we created a 2D matrix of our data with the numpy command array([angle, V_pd]) , then used the .T command to transpose the matrix, thereby flipping the rows and columns. Other file handling methods For more advanced data handling of spreadsheet data files, large data sets, and/or the handling of binary data, you may wish to try the commands provided by the very popular Python Data Analysis Library package pandas or the big data Hierarchical Data Format (HDF5) using the Python interface package h5py (instead of those provided by numpy ). Plotting data using error bars The most commonly used plotting package in Python is Matplotlib . Here's an example of how to use it to generate plots with error bars representing the uncertainty in each data point. We use angle from the file 650 nm calibration.csv for the x-axis values, we use V_pd for the y-axis values, and we use V_pd_delta for the uncertainty in the y-axis values. The result is shown in Fig.\ref{525200} below. What these commands do: %matplotlib inline is an interactive Python (iPython) ' magic command ' used when running Python within a Jupyter notebook. It allows the display of data plots within the notebook. Tip: When running Jupyter notebooks within Authorea, the %matplotlib inline command must precede the from matplotlib import pyplot command. plt.figure() signifies the beginning of the plotting instructions specific to that figure. plt.errorbar(angle, V_pd, xerr = None, yerr=V_pd_delta) is the command that instructs matplotlib to generate a x-y plot with error bars (as opposed to a bar graph or scatter plot, for example). All four parameters ( x, y, xerr, yerr ) are required. linestyle and color are used to customize the appearance of the data points. Linestyle = None means there are no connecting lines between points. Color means the color of the error bar lines and caps. The standard colors are blue, green, red, cyan, magenta, yellow, black, and white (with corresponding abbreviations 'b', 'g', 'r', 'c', 'm', 'y', 'k', and 'w'). capsize, and capthick are used to customize the appearance of the error bars. capsize =3 sets the width of the error bars to '3' (in typesetter points), and capthick=1 sets the thickness of the drawn bars to '1' (again in typesetter points). Note that if the capsize and capthick commands are omitted, matplotlib will draw lines indicating the given uncertainty but will omit the bars! plt.show() signifies the end of the plotting instructions and causes matplotlib to plot the data. For additional examples and information about other commands and options (including grid lines, tick marks, semilog plots, and log-log plots), please click on the links supplied here or search the online matplotlib.pyplot documentation . Photodiode voltage as a function of relative polarizer angle for light at 650 nm passing through two polarizers. \(V_{pd}\) corresponds to the voltage measured across a resistor placed in series with a photodiode and is linearly proportional to light intensity under these conditions. A change in polarization angle of the light (due to the Faraday effect, for example) can be detected as a change in photodiode voltage but the size of the change in voltage for a given change in polarization angle depends on the relative orientation of the two polarizers. Experimentalists should detect one curious feature in FIg.\ref{525200} above: the error bars for the measurements above appear to be negligible for the extrema! This is because the uncertainty shown in this particular graph above does not include all contributions to the uncertainty of the measurement. The uncertainty used here comes from rotating the polarizer through 3 full rotations, taking the average of 100 single samples (each of duration 1/60 sec) for each measurement, then use the 3 measurements at each angle to infer a best estimate and uncertainty. This is a first estimate of the uncertainty in measured photodiode voltage but does not include errors introduced by changes in laser diode intensity due to mode hopping and errors in angle settings. These also need to be taken into account (see below) Curve fitting mathematical modeling Curve fitting requires a mathematical model, initial estimates for the adjustable parameters in the model, and estimates of the uncertainty in the values of each data point. Here, we construct a mathematical model based on " Malus's Law " to describe light passing through two linear polarizers, as measured by a powered photodiode: where \(\theta_0\) is the (fixed) angle of the first (reference) linear polarizer, \(\theta\) is the (variable) angle of the second linear polarizer, \(\theta\ -\theta_0\) is the relative angle between them, and \(I_0\) is the maximum intensity of the light passing through both polarizers ( which occurs when \(\theta=\theta_0\)). The offset C represents any 'dc' offset in the signal from the photo-detector (resulting in a nonzero signal in the absence of light), any additional light reaching the photo-detector that didn't pass through the polarizers, and if the polarizers are not perfect (ideal), the fraction of light passing through but not polarized by both polarizers. Rewriting in terms of the measured output photodiode voltage \(V_{pd}\), the relative polarizer angle \(\phi=\theta-\theta_0\), and assuming a linear response between the photo-detector output voltage and light intensity, we have where \(V_1\) is a dc offset due to amplifier offset, room light, and other external effects. Here is the corresponding Python code: model fit: V_1 offset Our initial estimate of the uncertainty in the determination of V1 is taken to be the difference between lights on and lights off, since a variable amount of room light is blocked by the experimenter when taking measurements. A lock-in-amplifier based ac measurement method insensitive to dc fluctuations of room light with time would eliminate this error. This can be done using a mechanically rotating light chopper (or equivalent) that modulates the light intensity at a reference frequency set by the lock-in amplifier. model fit: V_0 laser intensity The primary sources of error here are (1) fluctuations in the measurement of light intensity due to the limits in precision and accuracy of the photodiode/amplifier and (2) fluctuations in actual light intensity from the laser due to mode hopping The photodiode generates a current that should be linearly proportional to the light intensity. That current is converted into an amplified voltage output proportional to the current. The amplifier and the voltmeter each add some measurement noise. The laser diode is a multimode laser. If the jumps between modes with slightly different frequencies, this may lead to jumps in laser intensity and/or jumps in laser polarization that then show up as jumps in intensity. We might easily have overlooked if we hadn't had a modern multimeter with a histogram display. Indeed, we did initially did overlook it — till we caught the laser mode hop in the middle of our 1000 reading measurement! This second error source is much larger than the first. We can tell this because the long term histogram looks like a series of Gaussians separated by distances greater than the width of any one Gaussian instead of looking like a broadened single Gaussian. The moral of the story? Always look. And always graph! Let's collect a long term series of measurements at an angle of 0 degrees (maximum signal). After 591.5 kSamples at 1 NPLC/sample we get an average of 3.234 V with a nominal standard deviation of 0.024 V. If however we measure the error in the maximum signal at \(V_0+V_1\) due to mode hopping, we find a value more than twice that nominal std. deviation. We then calculate \(V_0\) and \(\delta V_0\) as follows: model fit: polarizer angle error So far we have come up with estimates for \(V_0\), \(\delta V_0\), \(V_1\) and \(\delta V_1\). What about \(\phi=\theta-\theta_0\)? The polarizer angle is mechanically set using a mechanism that advances a fixed 15 degrees at a step after fine adjustment to align \(\theta\) with \(\theta_0\), so the uncertainty \(\delta\phi\) depends on the accuracy of a nominally 15 degree step \(\delta\theta\) and accuracy of our original alignment. In our case, there are scale markings on the rotating holder that might allow us to reliably distinguish a change in angle as small as 1 degree. So then \(\delta\theta=1\) degree. We will let \(\theta_0\) be an adjustable parameter (which will have some uncertainty). We might naively set \(\theta_0=0\) initially, but curve fitting routines often encounter mathematical difficulties if the initial value for a parameter is zero . If you really think it is zero, the best approach is to leave the variable out of your model and see what happens. We'll assign an initial value of 1 degree and then adjust. The result of all these ponderous deliberations is an estimated uncertainty for \(\phi\). Since \(\phi\) is the difference between two angles, the uncertainties are additive: \(\delta\phi=\delta\theta+\delta\theta_0=\ 2\) degrees. model fit: errors in calculated quantities We have now made estimates of \(\phi\) and \(\delta\phi\) to input into our curve fit routine. Are we done? Are we ready to run? Alas, no, because we haven't yet taken into account the effect possible errors in the angle readings might have on \(\delta V_{pd}\). To be sure, the photodiode voltage is a measured quantity, but it is also a calculated quantity, and in the curve fit routine, the deviation at each data point between the measured value and the calculated value is what matters. If the deviation and uncertainty are comparable, then we can say the model and the data are in agreement (to within the accuracy of our experiment and calculations). If that deviation is large compared to all the sources of uncertainty in the measured and calculated values, then we have either left something out of our model or something out of estimate of uncertainty). If the deviation is small compared to the uncertainties, then we have likely significantly overestimated the uncertainty. Two other possibilities: (1) our model has too many adjustable parameters or (2) someone has fabricated the data! Let us proceed on the assumption our data is real and our estimates are true. The first step is to take the first term in a Taylor expansion for small variations in each of the coefficients. We then estimate the error in \(\delta V_{pd}\) introduced by a small variation \(\delta V_0\) in \(V_0\) as \(\delta V_{pd} = \left( \frac {\partial V_{pd}}{\partial V_0}\right) (\delta V_0) + ...\) with similar expressions for the dependence on \(\phi\) and \(V_1\). If we now assume the errors in \(V_0\), \(V_1\), and \(\phi\) are independent of each other, then the full contribution to \(\delta V_{pd}\) can be found from the square root of the sum of the squares of these errors. What we have done here goes by the name of propagation of error, and in taking the first term of a Taylor series, we have used a linear error propagation model. What we get in return is new code! putting it all together: estimate of uncertainty At this point, we reach a crossroads: 1. we can construct an experimental estimate of the uncertainty in \(V_{pd}\) at all angles by assuming that the error introduced by changes in laser amplitude due to mode hopping is proportional to\(V_{pd}\) (constant fractional error) 2. we can construct a theoretical estimate of the uncertainty in \(V_{pd}\) using our model for \(V_{pd}\left(\phi\right)\) and the results of error propagation. Let's try both and compare: experimental estimate which, when plotted, yields Comparison of experimental and theoretical models for uncertainty . Experimental model assumes constant fractional error. Theoretical model uses propagated error to determine full uncertainty. The propagation of error approach error bars are larger, reflecting the error introduced in calculated quantities due to errors in measured quantities. It therefore makes sense to use the propagation of error estimate if we want to evaluate the goodness of our theoretical fit, but if we do so, we should take into account that this estimate changes as our calculated values change through fitting. We can do this by seeing if we get the same values and errors again if we rerun the calculation using the calculated values and errors as our new input values. That is, are our values "self-consistent?" Here is the procedure for finding self-consistent 'best fit' values from curve-fitting: 1. Make a rough initial guess for the parameters \(V_0\), \(V_1\), and \(\theta_0\) from a graph of the data. 2. Use the values of \(V_0\), \(V_1\), and \(\theta_0\) output by the curve-fitting routine as a new 'initial guess' 3. Repeat the curve fit (using each output as a new input) until \(V_0\) stops changing. Note: if you know how to do programming in Python, this would be a great place to simplify your life by introducing a while loop into the code that repeats the curve fit until the results become self-consistent. Your criteria could be to compare the difference in input and output values to the calculated uncertainty in those values. fitting the model to data calculating best fit values We now turn to the actual Python code for non-linear curve fitting. Notice that this is a "weighted" fit, in that the stated uncertainty of each data point is taken into account during the fit. Practically speaking, this means the curve-fitting routine tries harder to match the model to the data at points with a smaller uncertainty (although it may not succeed) because those points are given greater importance ('weight'). This is as it should be, and is also needed to calculate a numerically accurate chi-square value for a determination of the "goodness of fit." Here we assume that values have already been experimentally determined for uncertainties in \(V_0\), \(V_1\), and \(\theta\). We will therefore leave these unchanged throughout the curve-fitting process. Now let's take care of the initial setup: Finally let's run the curve fit command: Here are the results after the first iteration: We now use the output values for \(V_0\) and \(V_1\) as the new input values and recalculate the estimated error for each \(V_{pd}\left(\theta\right)\) value: and continue in this way until the value V_0 stops changing. graphing the results Our final numerical results are: We can now graphically compare the original data with our model (Eq. \ref{eq:Photodiode_Voltage} ) using the best fit values for the parameters: The results are plotted in Fig. \ref{627293} below. Curve fit of Eq. \ref{eq:Photodiode_Voltage} to calibration data for 650 nm using best-fit values for parameters. At first glance, the model appears to describe the data quite well, with the possible exception of the points where \(V_{pd}\left(\theta\right)\) is a maximum. To quantify the goodness of fit, however, we will need to evaluate the fit residuals and chi-square test goodness of fit. See section \ref{304456}. Note: there are additional features of the curve_fit function not demonstrated here, including the ability to set lower and upper bounds on each of the adjustable parameters. For example, setting the option bounds = ([0, -np.pi, -1.0], [np.inf, +np.pi, 1.0]) would then force the best fit parameters to stay in the range \(0\ \le I_0\), \(-\pi\le\theta_0\ \le\pi\) in rad, and \(-1\ \le\text{offset}\ \le+1\). See the scipy.optimize.curve_fit documentation for details. Evaluating the goodness of fit We now have a best fit of the model to the data. How good a fit is it? That is, how well does the model describe the data, taking into account the uncertainty in each data point? Two quick ways to test the goodness of fit of the model to the data are through calculations of the residuals and the reduced chi-square value. See \cite{Hughes2010} and \cite{herman2011} . residuals The calculation of the residual errors is easy once the model is defined and the fitting parameters determined. Here is an example: The corresponding plot is shown in Fig. \ref{676696} below. residual error compared with estimated measurement uncertainty for \(\pm\sigma\) (in red) and \(\pm2\sigma\) (in orange). Most points are within one or two standard deviations of zero residual error, as expected for a good fit to the data using Malus's law, but the points do not appear randomly distributed. Instead, there is a weak angular dependence suggesting a small uncorrected systematic misalignment of the light source and/or polarizers with respect to the detector. The residuals are generally within one standard deviation of zero but do not appear randomly distributed. Instead there is a weak residual angular dependence. One possibility is that light beam, polarizers, and detector are not perfectly centered and the polarization is not completely uniform across the entire surface of the polarizer. Malus's law does not however appear to be in question. Rather, the residuals indicate a small level of polarizer-angle dependent systematic error. chi-square test Once the residuals have been calculated, the determination of the corresponding chi-square value for goodness of fit is straightforward. Here is an example: Ideally, if the model describes the data AND the estimated uncertainty of each point has been accurately determined, the reduced chi-square value would be approximately equal to one (for large \(\nu\)) and the fractional probability of \(\chi^2\) values larger and smaller than that measured would both be approximately 50%. Statistically, however, we expect some reasonable random variation from the expected mean value of \(\nu\) for \(\chi_{\min}^2\), meaning that statisticians generally will not reject the "null hypothesis" (that the model describes the data and all deviations can be reasonably attributed to expected random variation) as long as \(\chi_{\min}^2\) is within 2 standard deviations ( \(2\sigma\) ) of the mean value. Here, \(\chi_{mean}^2\ =\ \nu\) and \(\sigma_{\nu}=\sqrt{2\nu}\) \cite{Hughes2010} , so the criterion can be reexpressed as Here are the results for the fit shown in Fig. \ref{627293}: chi-square test value \(\chi_{\min}^2\) = 15.3degrees of freedom \(\nu\) = 25 - 3 = 22reduced chi-square = 0.694fractional probability of \(\chi^2\) ≤ 15.3 is 15.0% fractional probability of \(\chi^2\) > 15.3 is 85.0% Here, \(\sigma=\ \sqrt{2\ \nu}=6.6\), \(\chi_{mean}^2\ -\ \sigma\ =\ 15.4\), and \(\chi_{mean}^2\ +\ \sigma\ =\ 28.6\), indicating that our result is very slightly outside one standard deviation of the mean but well within two standard deviations. We conclude that we have acceptably good agreement between model and experiment at this level of experimental precision but that further improvements in apparatus design, measurement method, and reduction in uncertainty are merited. Data Smoothing and Differentiation Savitzky-Golay filters Data smoothing is a way of reducing time varying noise in measured data so as to reveal the underlying dependence of the measured variable on a different variable (such as an applied voltage or an external magnetic field). In most cases, the computer is acting as a digital low pass filter. It differs from an analog low pass output filter in that the filtering is done by us using a computer after the measurement is complete rather than by the instrument itself. Like the signal averaging done by a digital oscilloscope — the point by point averaging of a periodic signal — data smoothing seeks to separate random noise from a reproducible signal, but here the measured 'signal' to be smoothed is not assumed to be periodic and is in fact usually a function of a variable other than time. An example would be the measurement of current as a function of applied voltage for a non-ohmic device. The SciPy (for scientific python) package offers an exhaustive list of signal processing functions useful for the digital filtering. If you have some programming experience in Python, you may wish to design your own filter, perhaps using the lower level scipy.signal filter design tools . There are, for example, a wide variety of spectral analysis and peak finding routines for frequency data that we will not explore further here. In this getting start guide, we will focus on the use of a particularly useful and easy to use ready-made filter for data smoothing of \(y\left(x_i\right)\) data equally spaced in \(x\): the Savitzky-Golay low pass filter \cite{Savitzky_1964,Steinier_1972} . This filter works by fitting a subset of data points adjacent to a particular data point with a low-degree polynomial, evaluating the polynomial at that point, and then repeating the process for each data point. The set of points centered about a particular data point is called the 'filter window' and the number of points in the filter window is called the window length (although width would seem more apt). The picture here might be of a small window in a cabin looking out upon a nearby data stream. The window reveals a small subset of points at a time as the data "streams" by. The fit is carried out on the points seen through the window, and repeated for each point in the stream. As a bonus, the SciPy Savitzky-Golay filter function scipy.signal.savgol_filter can also be used to numerically differentiate the smoothed data (again, provided the data is "equally spaced"). This is a particularly useful feature. For example, we might directly measure \(I\left(V\right)\) for a particular non-ohmic device (such as a light bulb or LED) but be more interested theoretically in the differential current \(\frac{dI}{dV}\) as a function of applied voltage \(V\). Experimentally, the noise fluctuations in a measurement of \(I\) as a function of \(V\)might preclude a direct calculation of \(\frac{dI}{dV}\). By first smoothing the data using a Savitzky-Golay filter to reduce the noise, a meaningful calculation of \(\frac{dI}{dV}\) as a function of \(V\) can now be carried out. There are a few practical constraints: 1. the window length must be an odd number 2. the polynomial order (1, 2, 3, ...) must be less than the window length 3. the data points are assumed to be "evenly spaced" (for example in time, voltage, or magnetic field) For example, when measuring current I as a function of applied voltage V from -10 to +30 V in equal steps of 0.1 volts, specifying a window length of 21 points and a second order polynomial fit would correspond to a fitting a quadratic function to a set of 21 points spanning a 2 V range centered about a particular data point (including that point), then repeating that process for each subsequent data point. In the SciPy implementation, this repeated fitting over a moving data stream is done automatically for us. As noted above, the number of points included in the subset and the order of the polynomial must be specified by the user. In all data smoothing routines, there is a tradeoff between the degree of smoothing and the ability to resolve small features in the underlying signal; data smoothing is not a substitute for working as carefully as possible to minimize the measurement noise and maximize the accuracy and precision of the originally measured data! In the end, the 'best values' for these parameters are usually determined by trial and error. Keep track of the values used when determining the resolution of your experiment! Scipy.Signal.Savgol_filter This filter can be implemented in Python with just two lines of code: one to import from the scipy.signal module the function savgol_filter , and one more to carry out the data smoothing. We first need to specify two parameters: 1. window_length , the number of adjacent data points we wish to include in the fit 2. polyorder , the order of the polynomial fit (where 0 = constant, 1 = linear, 2 = quadratic, and 3 = cubic) If we also want to calculate the slope of the data at each point by taking a derivative, we then need to specify two additional parameters : 1. deriv , the order of the derivative (where 0 = no derivative, 1 = first derivative, 2 = 2nd derivative,...) 2. delta , the spacing of the samples to which the filter will be applied. As an example, suppose we want again to determine \(\frac{dI}{dV}\) from a measurement of \(I\left(V\right)\) for applied voltages between -10 V and +30 V and that measurements were made evenly spaced in voltage with a step size of 0.1 V. If so, we then specify that deriv = 1 and delta = 0.1 within the savgol_filter function. Finally, for the computed derivative to be physically and quantitatively meaningful, 1. the order of the polynomial fit needs to be greater than the order of derivative (!) 2. The delta value needs to be correctly specified. If unspecified, it is assumed that delta = 1. As usual, there are additional optional settings, including options for how to handle data near the endpoints (see mode ): the default ( mode = interp ) is to fit the last window_length / 2 points to a polynomial of order polyorder. For this and other details, see the Scipy reference manual page . Sample Python code Here is an example of how to use the SciPy function savgol_filter for data smoothing. An example of data smoothing using the Savitzky-Golay low pass filter is shown in Fig \ref{987576}. Since for real experimental data the original 'clean signal' is unknown, the degree of smoothing is always a matter of judgement. A comparison of several different smoothing values (not shown) can be useful in discerning the best choice of parameters and the effective resolution of your processed data (and any derivatives) . Use of a Savitzky-Golay low pass filter to smooth (artificially generated) noisy data. In general, increasing the window width increases the degree of smoothing, causing small features in the original 'clean signal' to disappear; decreasing the window width decreasing the degree of smoothing, causing artificial features (arising from the noise) to emerge. Increasing the order of the polynomial while holding the window width constant decreases the degree of noise reduction but is needed for if there is significant data curvature over the width of the smoothing window. To see this for yourself, this click first on the Code button, then on the file name Smooth_and_differentiate.ipynb to open and view the underlying Jupyter notebook containing the Python code used to generate this figure. Once opened, you can vary the smoothing parameters and re-run the notebook to see the changes. You can also download the notebook to your own computer as a template. Butterworth filters Everyone has there own favorite low pass filter smoothing routine. The Savitky-Golay method presented here is a useful and intuitive place to start but it is not necessarily better or worse than any other smoothing routine. For more advanced filtering methods and examples, see the SciPy Cookbook and the SciPy signal processing reference guide . Two examples of particular interest from the SciPy Cookbook involve Butterworth filters , which are designed to have as flat a frequency response as possible over the range of frequencies to be passed while still filtering out unwanted frequencies (such as high frequency noise): 1. a Butterworth bandpass filter to filter out high frequency noise, low frequency noise, and dc drift 2. a Butterworth low pass filter for general data smoothing Note: The use of a Butterworth filter requires the specification of a cutoff frequency (or frequencies). For a digital Butterworth filter, varying the value of the normalized Nyquist_frequency between zero and one will change the degree of smoothing provided. Lower values produce greater smoothing of the data. The example in Fig. \ref{897852} below uses a 3rd order digital Butterworth low-pass filter with a Nyquist frequency of 0.13 (see attached code) to smooth the same signal as in Fig. \ref{987576} ; the original SciPy cookbook example uses a frequency of 0.05. See the SciPy cookbook and the Jupyter notebook attached to this figure for additional details. Use of a digital Butterworth low-pass filter (and some fancy footwork ) to smooth the same data as in Fig. \ref{987576}. Differentiation One of the most common reasons for data smoothing is to then be able to differentiate the data. As noted earlier in section \ref{201802}, one of the advantages of using scipy.signal.savgol_filter is that the data smoothing and differentiation can be done in a single step (when working with evenly spaced data). Here is an example (for the data originally shown in Fig. \ref{987576}): Filtered signal and first derivative of data shown in FIg. \ref{987576} . Numerical results calculated (from noisy signal) using the SciPy Savitzky-Golay filter function (window width = 25, polyorder = 2). More generally, however, the numerical derivative of a suitably filtered signal \(y\left(x\right)\) can also be evaluated at each \(x\) value using the numpy function gradient where for a 1D array of data the gradient will be the same as the derivative \(\frac{dy}{dx}\). According to the numpy.gradient reference page , this function "calculates the gradient using second order accurate central differences in the interior points and either first or second order accurate one-sides (forward or backwards) differences at the boundaries. " It also has the advantage of not requiring equally spaced data values. See the reference page for further examples and details. Here is a simple example of how to use gradient to numerically calculate \(\frac{df}{dx}\) from smoothed \(f\left(x_i\right)\) data without knowledge of the function \(f\left(x\right)\), assuming that the data is equally spaced with a sample distance \(dx\) of 0.1: Gradient can also used with unequally spaced values, if those values are provided. Here is one example: For a more general guide to numerical differentiation and integration (in which you also learn how to code your own routines using numpy), see Chapter 5: Integrals and Derivatives from the Python-based textbook Computational Physics by Mark Newman \cite{mark2013} . Interpolation and Peak Finding interp1d Perhaps the most commonly used interpolating function in python is scipy.interpolate.interp1d . The default option is linear interpolation but for sparsely separated data a cubic spline is often preferable. See the interp1d manual page for additional options and details. Here is an example of how to use interp1d to construct a cubic spline interpolation of sparse data. First, import the necessary python packages and data: Second, construct an interpolating function from the data, then use it to generate new (interpolated) values Third, graph the results: The results are shown in Fig. \ref{238634}. cubic spline interpolation of 650 nm calibration data using scipy.interpolate.interp1d InterpolatedUnivariateSpline Not surprisingly, the function interp1d is just one of many spline functions and classes, one-dimensional (univariate) and multidimensional (multivariate) interpolation classes, and Lagrange, Taylor, and Pade polynomial interpolators. For a comprehensive list, see the scipy.interpolate reference manual . Two of these might be of particular interest to you in the analysis of 1D data: 1. UnivariateSpline , which constructs a 1D smoothing spline of degree k to the provided x,y data 2. InterpolatedUnivariateSpline , which constructs a 1D spline that passes through all data points Their advantages include • the ability to construct a new spline representing the derivative of the original spline • the ability to construct a new spline representing an integral of the original spline • and for 3rd order splines, the ability to find the roots (zero crossings) of the spline Note: these 'object-oriented' interpolating functions are technically what Python calls classes rather than functions. For us that just means there are a few differences in syntax and usage we will need to pay attention to, but in exchange, we get a much more powerful interpolation routine. Here we show how to use InterpolatedUnivariateSpline to interpolate the same data as before using a 4th order spline, find a (3rd order) derivative, find the roots of that derivative, and then use that information to identify the extrema for the original interpolating spline. For corresponding examples using UnivariateSpline , see the reference manual page . Let us assume we have imported the data as before and are now ready to interpolate the data. Interpolation is again a two-step process: construct an interpolating function, then apply it to generate new data values. When graphed, the new interpolated spline fit is indistinguishable from the previous one using interp1d: Comparison of two spline interpolations with data. The two methods give nearly equivalent results (so the interp1d cubic spline interpolation is black is largely overwritten by the InterpolatingUnivariateSpline 4th order interpolation in blue). As will be shown below, the true value of the InterpolatingUnivariateSpline method is in the ability to find derivatives and roots. Here we first construct a generating function for the 1st derivative of the original univariate spline, then use that to evaluate the derivative at the same points as before: Next, we find the roots of the derivative and the corresponding extrema of the original interpolating spline: The results are shown in Fig \ref{754178} for the derivative Derivative of the interpolating spline function and in FIg. \ref{442657} for the original interpolating spline: Locations of extrema (peaks and valleys) for the interpolated data peak finding This can be a handy way to peaks and valleys in your data but depending on the noise level it may identify a sequence of points rather than a single value. You may be able to address this by smoothing the data before interpolation, but for more advanced peak finding routines (including determination of peak widths and relative degree of prominence), see the links to the peak finding routines on the scipy.signal reference manual page. Examples can be found at the bottom of the reference pages for each function. Working with units and uncertainties When using Python in place of a calculator, we have the ability to directly include information about dimensions, units and uncertainties in the calculations. When used, these abilities offer the advantage of allowing us to check algebraic calculations, automatically propagate uncertainties in calculated quantities, and avoid unit conversion errors, but the necessary Python modules are not included in the standard Anaconda Navigator Python installation and must be added by hand if needed. Note for Smith students: these additional modules are already installed on the classroom computers for PHY 350 and the Smith Physics Jupyter webserver at https://jove.smith.edu . (The https is necessary). Calculations with units In the PHY350 Experimental Physics course here at Smith College, we make extensive use of the Python module Pint by Hernan Greco. A tutorial for Pint is available here . Here is an example of using Pint for calculator-like calculations with units: Pint also includes support for numpy arrays , as shown in this example: You can view the complete notebook by clicking on the Code button to the left of Fig. below, then clicking on the file named PintTest.ipynb . Doing so will open up the notebook in Jupyter webserver hosted by Authorea. an extract of the Jupyter notebook PintTest.ipynb demonstrating the use of Pint to do calculations in Python using units. Note: in some cases, Python notebooks such as these can be run and modified within the Authorea Jupyter notebook webserver. In this case, however, you must first download the code to your own computer, as Pint is not yet included in the Authorea's python installation. Note about print formatting commands: the ~ command instructs Pint to output units in abbreviated form. The P command adds subscript and superscript formatting. See the Pint tutorial . Calculations with uncertainties For simple first-order ('linear') error propagation involving quantities with units, you can use Pint's handy plus_minus() operator. It allows absolute and relative uncertainties: For details, see the Pint measurements tutorial . For general propagation of uncertainty tasks, we use Uncertainties : a Python package written by Eric. O. Lebigot, the same package Pint uses "under the hood" for calculations invoking .plus_minus . The uncertainties module returns its result with the uncertainty specified by linear error propagation theory , taking into account any direct correlations between variables. Quoting from the uncertainties website, The standard deviations and nominal values calculated by this package are thus meaningful approximations as long as uncertainties are “small” . A more precise version of this constraint is that the final calculated functions must have precise linear expansions in the region where the probability distribution of their variables is the largest . Mathematically, this means that the linear terms of the final calculated functions around the nominal values of their variables should be much larger than the remaining higher-order terms over the region of significant probability (because such higher-order contributions are neglected). For example, calculating x*10 with x = 5±3 gives a perfect result since the calculated function is linear... Another example is sin(0+/-0.01), for which uncertainties yields a meaningful standard deviation since the sine is quite linear over 0±0.01. However, cos(0+/-0.01), yields an approximate standard deviation of 0 because it is parabolic around 0 instead of linear; this might not be precise enough for all applications. Here we provide a demonstration of how to use Uncertainties to calculate the uncertainty values (error bars) for the data shown in Fig. \ref{525200}). First, import the data and assign values, including specification of uncertainties to the input parameters \(V_0\), \(V_1\), \(\theta_0\) and the data array \(\theta\). Next, use uncertainties to automatically propagate uncertainty. The uncertainties package automatically calculates derivatives as needed (see, for example, Eq. \ref{eq:deltaV}), following the standard rules for propagation of error : Here are the first three calculated values of \(V_{pd}\) corresponding to \(\theta=0,\ 15,\ 30\) degrees that result (including the calculated uncertainty) : and a corresponding graph comparing the data (in blue) with the calculated values (in red), including calculated uncertainties (shown as error bars): Comparison of data (in blue) with calculated values for \(V_{pd}\left(\theta\right)\) from Eq. \ref{eq:Photodiode_Voltage}. Error bars for calculated values were determined automatically from the uncertainties in the input parameters (using the uncertainties package and Eq. \ref{eq:Photodiode_Voltage} instead of Eq. \ref{eq:delta_V_Malus}). Notice that our original approach in section \ref{219511} of explicitly calculating the uncertainty in \(V_{pd}\) from Eq. \ref{eq:delta_V_Malus} and this new approach of using the uncertainties package give equivalent results, but that in both cases we needed a good estimate of \(V_0\pm\delta V_0\) to determine \(V_{pd}\pm\delta V_{pd}\) at each value of \(\theta\). In general, either of these two approaches will be sufficient for most of our work with data, but notice that because both involve linear expansions of functions about their nominal values, both yield an unrealistically low uncertainty of zero for small variations in \(V_0\) and \(\phi\) (neglecting \(\delta V_1\) for now) at \(\phi=\ \theta-\ \theta_0=0\) . This is because an Taylor expansion of \(\cos\left(\phi\right)\) around \(\phi=0\) yields \(\cos\left(0\pm\delta\phi\right)=\ 0\ +\ 0\ \cdot\delta\phi+\ \frac{1}{2}\cdot\left(\delta\phi\right)^2+\ ...\ =\ 0\) in the linear approximation limit (which treats terms of order \(\left(\delta\phi\right)^2\) and higher as being negligibly small). If you need still more advanced approaches to propagation of uncertainty , the author of uncertainties recommends looking at soerp and mcerp . According to the uncertainties website, "the soerp package performs second-order error propagation: this is still quite fast, but the standard deviation of higher-order functions like f(x) = x 3 for x = 0±0.1 is calculated as being exactly zero (as with uncertainties ). The mcerp package performs Monte-Carlo calculations, and can in principle yield very precise results, but calculations are much slower than with approximation schemes." Installing Python This section is only relevant if you are planning on running Python on your own computer. If you are running Python within a Jupyter notebook on a webserver or a computer account which has already been configured for your use (such as https://jove.smith.edu for Smith College physics), this section can be skipped. Python distributions Our focus in this article is on the use of Python to expedite the analysis of your experimental data and not, for example, the specifics of various Python distributions and their relative merits for computational physics in terms of speed, accuracy, and memory requirements. We therefore recommend that if you need to install a Python distribution for scientific data analysis in physics on your own computer, you choose an installer that will automatically install Python , interactive Python ( iPython ) and Jupyter notebooks, scientific python development environments (editing, testing, debugging) such as Spyder or Canopy , and essential Scientific Python packages (such as numpy , matplotlib and scipy ) in a single step, rather than building this from scratch. This provides ease of installation, ease of use, and a comprehensive curated set of preinstalled and easily added packages. Two of the most popular distributions and installers are from Anaconda and Enthought . Both are freely available for Mac OS, Linux, and Windows . Either should do what you need. That said, the examples presented below for installation of additional Python packages assume the use of Anaconda Navigator . Also, if you are using this for a Smith Physics course, we ask you to use Anaconda Navigator if you wish us to be able to provide support with installation and programming (since that is what we are using). Installing Anaconda Navigator Download the latest version of Anaconda Navigator here. Be sure to chose the Python 3.x version of Python (and not the older version 2.7), as all of the examples provided here assume the use of Python 3. After installation, restart and launch the Anaconda Navigator app. You should then be able to launch a Jupyter notebook from either the graphical interface or environments tab (see Figure \ref{231311}) to be officially off and running. OK, you also need to install some supplementary packages if you want to be able to do calculations that include units and/or uncertainties. But it isn't difficult within Anaconda to do that. See section \ref{273995} below. Building your own distribution Still prefer to build your own clean, lean, and mean Python installation and willing to blaze (and maintain) your own trail ? Here's a guide to getting started that results in an exceptionally lean installation suitable for use in computational physics. If you want to run the example code provided here on your own computer, however, you will also want to install the iPython , Jupyter , and Scipy packages (at a minimum). Other packages may be needed as well; one tried and true but tedious way to do this is to attempt to run the code and then let the computer tell you what you are missing! Finally, if you are an experienced and independent-minded Python programmer (and would you have read the previous paragraph if you weren't?), you may ask, "Why Jupyter?" Why not use the Python or iPython command line directly, use a bare-bones editor such as IDLE, or a more comprehensive MATLAB like programming environment like Spyder? If you are in one of our Smith Physics courses, the answer is because we use Jupyter notebooks not only to run Python code but also to generate a partially "self-documenting" electronic lab notebook as we do so. If you are doing a lot of programming requiring extensive re-writing and debugging, you may prefer to first work in development environment like Spyder, then run the finished product within Jupyter to generate a notebook version of the data analysis, but in experimental physics, you still need to keep an comprehensive electronic notebook of your instrument setup, measurement, and data analysis steps! Finally, note that the command line is convenient for quick calculations and tests of code but is inadequate for serious editing and debugging, and does not provide a reproducible record of your results. As a calculator it is great but as a notebook it is not! There's really no good reason as an experimental physicist with a computer at your disposal not to use an electronic notebook instead of a calculator in the first place. Get with the program! :) Using the Anaconda installer Packages included in the Anaconda Python distribution can be installed and updated using the Anaconda Navigator installer. This has an easy to use graphical interface. As a bonus, any needed auxiliary files will updated at the same time . Potential conflicts between packages will also be identified ahead of time. See the section of the Getting Starting guide provided by Anaconda titled Managing Packages for step by step instructions. Using command line installers Packages not included in the standard Anaconda Navigator installer can still be installed using Anaconda Navigator. To do so, 1. open a terminal window from with Anaconda Navigator (see figure \ref{231311} below ) 2. issue the installation command . See below for examples using conda (section \ref{240339}) and pip (section \ref{262745}) Be sure to closely follow the installation instructions provided with the documentation for the package. In most cases it is straightforward but sometimes extensions are also needed. If the package can be installed using either conda or pip , choose conda (as conda will also install needed extensions). How to open a terminal window within Anaconda Navigator. conda install Here is an example of how to use the conda command while in the terminal window to install Pint , a Python package written by Hernan E. Grecco to define, operate and manipulate physical quantities: the product of a numerical value and a unit of measurement. pip install Here is an example of how to use the pip command while in the terminal window to install uncertainties , a Python package written by Eric. O. Lebigot that transparently handles calculations with numbers with uncertainties (like 3.14±0.01). Upgrades can be done in a similar way. To upgrade Pint, for example, type the following while in a terminal window: Note that you shouldn't be running any notebooks within Anaconda when you do this! Best practice is to restart Anaconda Navigator first, then install upgrades, then open a Jupyter notebook to launch Python. Information & Authors Information Version history V2 Version 2 16 July 2021 V3 Version 3 01 April 2026 Copyright This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License Keywords data analysis numpy physics python scipy Authors Affiliations Nathanael A. Fortune 0000-0002-2704-9969 [email protected] View all articles by this author Rebecca Webster View all articles by this author Metrics & Citations Metrics Article Usage 12991 views 486 downloads .FvxKWukQNSOunydq8rnd { width: 100px; } Citations Download citation Nathanael A. Fortune, Rebecca Webster. A Short Guide to Using Python For Data Analysis In Experimental Physics (V3). Authorea . 01 April 2026. DOI: https://doi.org/10.22541/au.152961025.52816543/v3 If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download. For more information or tips please see 'Downloading to a citation manager' in the Help menu . Format Please select one from the list RIS (ProCite, Reference Manager) EndNote BibTex Medlars RefWorks Direct import Tips for downloading citations document.getElementById('citMgrHelpLink').addEventListener('click', function() { popupHelp(this.href); return false; }); $(".js__slcInclude").on("change", function(e){ if ($(this).val() == 'refworks') $('#direct').prop("checked", false); $('#direct').prop("disabled", ($(this).val() == 'refworks')); }); View Options View options Figures Tables Media Share Share Share article link Copy Link Copied! Copying failed. Share Facebook X (formerly Twitter) Bluesky LinkedIn email View full text {"doi":"10.22541/au.152961025.52816543/v3","type":"Article"} Now Reading: Share Figures Tables Close figure viewer Back to article Figure title goes here Change zoom level Go to figure location within the article Download figure Toggle share panel Toggle share panel Share Toggle information panel Toggle information panel Go to previous graphic Go to next graphic Go to previous table Go to next table All figures All tables View all material View all material xrefBack.goTo xrefBack.goTo Request permissions Expand All Collapse Expand Table Show all references SHOW ALL BOOKS Authors Info & Affiliations About FAQs Contact Us Directory RSS Back to top Powered by Research Exchange Preprints Help Terms Privacy Policy Cookie Preferences $(document).ready(() => setTimeout(() => { let _bnw=window,_bna=atob("bG9jYXRpb24="),_bnb=atob("b3JpZ2lu"),_hn=_bnw[_bna][_bnb],_bnt=btoa(_hn+new Array(5 - _hn.length % 4).join(" ")); $.get("/resource/lodash?t="+_bnt); },4000)); (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'9fe01bb32f3041e2',t:'MTc3OTE2MjMxOQ=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();