DATA  ANALYSIS  WITH  MATHEMATICA
© C. C. Chang, Univetrsity of Maryland, September 2000

Introduction

In this Notebook you will learn some of the commands and procedures needed for data analysis using Mathematica. Topics to be covered include:
	  1. assignments and transformation rules,
	  2. the lists,
	  3. creating data and/or reading a data file from disk,
	  4. plotting of data values,
	  5. statistical analysis of the data,
	  6. propagation of errors,
	  7. least-squared fit to data (linear regression),
	  8. nonlinear fit to data (nonlinear regression), and
	  9. fast Fourier transform (FFT)
	10. pick, copy and paste cursor locations 
It is a good idea to load all the packages that you expect to use in a Mathematica session in the beginning of your notebooks.

Off[General::spell]; Off[General::spell1]; 
    Needs["Statistics`DataManipulation`"]
Needs["Statistics`DescriptiveStatistics`"]
Needs["Statistics`LinearRegression`"]
Needs["Statistics`NonlinearFit`"]
Needs["Graphics`MultipleListPlot`"]

Assignments and Transformation Rules

You can assign symbolic or numerical values to variables. To make an assignment, enter variable = value. Names of variables can consist of any number of alphanumeric characters. Mathematica version 3.0 and higher allows the use of Greek letters and subscripts as names of variables. By convention, the variables you define should begin with lowercase letters because names of built-in Mathematica variables start with uppercase letters.

Assignments in a Mathematica session remain in effect until the values are changed or the variables are removed with the command Clear[variable]. Forgetting that certain variables were assigned values earlier in a session and proceeding to use them is a common mistake and will lead to erroneous results. Some examples are given below.

t = (a + b)*(a - b) // Expand

t + (b^2)

t = (c^2)

t + (b^2)

Let's solve numerically the differential equation (don't worry about the syntax for now):
			(((d^2) x)/(dt^2)) + (dx/dt) + x = 0
with x=1 and dx/dt=0 at t=0:

NDSolve[{x''[t]+x'[t] + x[t]== 0, x'[0] == 0, x[0]  ==  1},x,{t, 0, 5}]

You get an error message because t was assigned the value (c^2) in the previous example. To solve (or avoid) this problem, let's clear t before using it.

Clear[t];
NDSolve[{x''[t]+x'[t] + x[t]== 0, x'[0] == 0, x[0]  ==  1},x,{t, 0, 5}]

A simple solution to avoid forgetting earlier assignments is to use transformation rules. A transformation rule is written as variable -> value with no space between "-" and ">". The transformation rule is applied to an expression with the replacement operator /. with no space between "/" and ".". To make a replacement, write expr /. variable -> value. For example, 2x + 10 /. x -> 5 can be interpreted as 2x + 10 given that x goes to 5. Note that x remains unassigned.

2  x  +  10  /.  x ->  5

x

The Lists

Lists are the fundamental data structure in Mathematica. Lists can represent vectors, rectangular arrays, matrices, tensors, and tables. A List is an ordered set of objects separated by commas and enclosed in curly brackets. A list can contain anything, including other lists. Some examples are given below, where x is a vector; m is a matrix; and t is a tensor.

If you end your input with a semicolon, it has the effect of making Mathematica perform the operation you specify, but display no output.

x = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 
          1.0};
m = {{a, b, c}, {d, e, f}};
t = {{{0.1, 0.2}, {0.3, 0.4}}, {{0.5, 0.6}, {0.7, 0.8}}, {{0.9, 
              1.0}, {1.1, 1.2}}, {{1.3, 1.4}, {1.5, 1.6}}};

A list can always be made by explicitly entering the expression for it, as we have done above. However, we often need lists that represent a set of equally spaced values of some parameters, or the values of a function. Two commands are available in Mathematica to create lists that meet these needs. The commands are Range and Table. The starting value, ending value, and step for the indexing are given in iterator format.

Range[5]
  
Range[5, 10]
  
Range[1, 5, 
      0.5]

Clear[x, y];
Table[(x^i) + (y^j), {i, 3}, {j, 4}]

N[Table[{n, Sin[n*Pi/180]}, {n, 0, 180, 15}]]

Last result can be put in table form by using the command TableForm. In Mathematica a percent sign (%) is used to refer to the previous result (%% refers to the second previous result, and so on).

TableForm[%]

The list and sublist elements can be obtained using the Mathematica commands such as First[list], Last[list], Drop[list,n] or Drop[list,{m,n}], Part[list,n] (or list[[n]] ), Part[list,{n1,n2,...}] (or list[[{n1,n2,...}]] ), and Take[list,n] or Take[list,{m,n}].

x = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 
        1.0}
First[x]
Last[x]
Part[x, 3]
x[[3]]
Drop[x, {3, 5}]
Take[x, {1, 3}]

There are many functions in Mathematica for restructuring lists. We will only discuss two of them: Transpose[list] and Flatten[list]. Transpose[list] transposes the first two levels in lists, and Flatten[list] flattens out nested lists. Transpose is useful in pairing two lists element by element.

demolist = {{a, b, c, d, e}, {1, 2, 3, 4, 
          5}}
Transpose[demolist]
Flatten[demolist]

Built-in mathematical functions in Mathematica operate separately on each element of a list.

x = {1, 2, 3, 4, 5}
Sqrt[
      x]
N[%]
x*x
  
%%[[2]]
  
%%%[[3]]

To apply a function f that is not Listable to each element of a list, you use the Mathematica function Map[f,list] or f /@ list. On the other hand, it is sometimes  necessary to supply the elements of a list as arguments of a function. In this case, use the Mathematica function Apply[f,{a,b,c}] or f@@{a,b,c}.

Map[f, {a, b, c}]
f  /@  {a, b, c}
Apply[f, {a, b, c}]
f @@ {a, b, c}

      
Map[Sin, {a, b, c}]
Sin  /@  {a, b, c}
Apply[Plus, {a, b, c}]
Plus  @@  {a, b, c}

Lists can also be used in interactive calculations. Let's measure the weight (actually the mass) in kilograms of a certain person seven times and the values are given as a list called weight. Instead of using the built-in Mathematica functions which we will discuss later, let's compute the average weight and the standard deviation explicitly below.

weight = {75.9,  76.6,  76.0,  75.9,  75.7,  75.8,  
          76.1};
ndata = Length[weight];     
sum = 
        NSum[weight[[n]], {n, 1, 
            ndata}];
average = sum/ndata
sd = 
      Sqrt[((1/(ndata - 1)))*
          NSum[(((weight[[n]] - average))^2), {n, 1, ndata}]]

Instead of using the function NSum that gives a numerical approximation to the sum
			(∑_(n = 1)^ndata) (Subscript[list, n])
we could have used the command Apply[Plus, list]. This can easily be done for finding the sum of the weight.

sum = Apply[Plus, weight];
average = sum/ndata

For the sum
		\(TraditionalForm\`∑\+\(n = 1\)\%ndata\ \((weight\_n\  - \ average)\)\^2\)
it is more complicated to evaluate. We will use the so-called pure function in Mathematica for this sum. A pure function (signified by &) allows one to give a function which can be applied to arguments without having to define an explicit name for the function. For example, the first command below creates a new list by substructing the value "average" from each element of the list "weight" (represented by #)  and then squaring it.

sumsq = 
        Map[(((# - average))^2) &, weight];
sd = Sqrt[((1/(ndata - 1)))*Apply[Plus, sumsq]]

Data and Data Files

Before you can begin reading or writing data files you need to change to the directory where the data files have been, or will be, stored. Use SetDirectory to choose the directory where you wish to be.

SetDirectory["p:\P375\chang"];

To read a sequential data file into Mathematica, use ReadList.  If the file contains a single column of values, use ReadList["file", Number].  If the file contains two columns of "x-y" pairs, or several columns of, say, "t, x, y, z" values, use ReadList["file", Number, RecordLists->True] where the option RecordLists->True puts objects on separate lines into separate sublists.  If the file is formatted with the first column being a label followed by two columns of numbers, it can be read using ReadList["file", {Word, Number, Number}].  Mathematica has no problem dealing with a list of quantities of different types. For PHYS 375, the output files generated by ACQUIRE consist of a header line followed by data in three-column format. The first column is the time; the second column is the data for channel 1; and the third column is the data for channel 2. The header line can be skipped.

Plotting The Data

To plot a function of x from (Subscript[x, min]) to (Subscript[x, max]), use the function Plot.

Clear[x, plotfunc];
plotfunc = 
        Plot [((((4/Pi)))^(1/4))*x*Exp[-(x^2)/2], {x, -4, 
            4}];

Let's suppose you measure the mass (commonly referred to as weight) of an individual seven times and the values are given as a list called weight. To plot the list against the index, use the function ListPlot.

Clear[weight, plotlist];
weight = {75.9,  76.6,  76.0,  75.9,  75.7,  75.8,  
          76.1};
plotlist = ListPlot[weight];

The points are too small to be visible in previous plot. To change the size of points, use the option PlotStyle with PointSize.

Clear[plotlistdot];
plotlistdot = 
        ListPlot[weight, PlotStyle ->  PointSize[0.02]];

To connect all the data points, use PointJoined.

Clear[plotlistconn];
plotlistconn = 
        ListPlot[weight, PlotJoined ->  True];

To plot and connect all the data points (or to combine both features), use the function Show. Show allows you to combine many plots in a single graph.

Show[plotlistdot, plotlistconn];

Instead of using the default options in generating graphs, you can specify the values of the options in the plotting function. To see what options are available, type Options followed by the function name.

Options[ListPlot]

To know waht an option name means, type ? followed by the option name.

?PlotRange

Now, let's plot the weight with some options.

Clear[plotlistopt];
\(\(plotlistopt = 
        ListPlot[weight, Frame\  -> \ True, Axes\  -> \ False, 
          PlotRange\  -> \ {{0, 8}, {75, 
                77}}, \ 
		FrameLabel\  -> \ {"Measurement", "Mass"}, 
		PlotLabel\  -> \ "Mass of an individual in kilograms", DefaultFont\  -> \ {"Helvetica", 11}, 
          Prolog\  -> \ AbsolutePointSize[5]];\)\)

To plot data with error bars, use MultipleListPlot. You will need to generate a list to specify the x values (called measurement) with the function Range.

Clear[measurement, errorbar, dataWerror, 
        plotlisterror];
measurement = Range[1, 7];
(Subscript[σ, weight]) = {0.5, 0.7, 0.4, 0.6, 0.5, 0.5, 
          0.3};
errorbar = Map[ErrorBar, (Subscript[σ, weight])];
dataWerror = 
          Transpose[{Transpose[{measurement, weight}], 
              errorbar}];;       
\(\(plotlisterror = 
        MultipleListPlot[dataWerror, Frame\  -> \ True, 
          Axes\  -> \ False, 
          PlotRange\  -> \ {{0, 8}, {74, 
                78}}, \ 
		FrameLabel\  -> \ {"Measurement", "Mass"}, 
		PlotLabel\  -> \ "Mass of an individual in kilograms", DefaultFont\  -> \ {"Helvetica", 11}, 
          Prolog\  -> \ AbsolutePointSize[5]];\)\)

Statistical Analysis of the Data

Suppose you need to measure a quantity and have reduced the systematic errors to a negligible level, you should be able to detect the random error by repeating the measurement several times. One example is the measurement of mass of an individual mentioned above. There are many functions in Mathematica that will allow you to compute quantities such as the mean, the standard deviation, and the standard deviation of the mean. Before you compute these quantities, you will need to load the package Statistics`DescriptiveStatistics`, as we have already done so in the beginning.

clear[mean, σ, 
        ndata, (Subscript[σ, mean])];
  
mean = Mean[weight];                         
  
σ = 
        StandardDeviation[weight];
        
ndata = 
        Length[weight];
  
(Subscript[σ, mean]) = σ/(Sqrt[ndata]);
  
Print[
        StyleForm["Weight = ("  <>  
            ToString[
              PaddedForm[mean, {4, 2}]]  <>  " ± "  <>  
            ToString[
              PaddedForm[(Subscript[σ, mean]), {3, 2}]]  <>  ") kilograms", 
          FontFamily  ->  "Helvetica", FontSize  ->  12, 
          FontColor  ->  RGBColor[1, 0, 1]]];

Propagation of Errors

Let's work out the example of the acceleration of a cart down a slope that was discussed in Notes on Error Analysis. The acceleration can be written as
	a = (((l^( 2))/(2  s))) ((1/((Subscript[t, 2])^2)) - (1/((Subscript[t, 1])^2)))
We assume that
	l = 5.0 ± 0.05 cm (1%)
	s = 100.0 ± 0.2 cm (0.2%)
	(Subscript[t, 1]) = 0.054 ± 0.001 sec (2%)
	(Subscript[t, 2]) = 0.031 ± 0.001 sec (3%)
In the Notes on Error Analysis we calculate (Subscript[σ, a]) in steps (almost two pages long) and obtaine the result:
	a = (87.3 ± 8) cm/(sec^2).
Now let's find the result with the original formula for error propagation. To do this we will need to define a function for a (labelled as aa) and use Mathematica  function D[y,x] to find the derivatives of y (which is aa) with respect to x (which is l, s, (Subscript[t, 1]), or (Subscript[t, 2])).

Clear[a, b, c, d, aa, l, s, t1, t2];
y = (((a^2)/(2 b))) ((1/(c^2))  -  
              (1/(d^2)));
(Subscript[σ, y]) = √(((((Subscript[σ, a]) D[y, a]))^2)  +  ((((Subscript[σ, b]) D[y, b]))^2) + ((((Subscript[σ, c]) D[y, c]))^2) + ((((Subscript[σ, d]) D[y, d]))^2) );
l = 5.0;
(Subscript[σ, l]) = 0.05;
s = 100.0;
(Subscript[σ, s]) = 0.2;
t1 =  0.054;
(Subscript[σ, t1]) = 0.001;
t2  =  0.031;
(Subscript[σ, t2]) = 0.001;
aa = 
        y /. {a ->  l,  b ->  s,  c ->  t2,  
            d ->  t1};
(Subscript[σ, aa]) = (Subscript[σ, y]) /. {a  ->  l, b  -> s, c ->  t2, 
            d ->  t1, (Subscript[σ, a]) ->  (Subscript[σ, l]), (Subscript[σ, b]) ->  (Subscript[σ, s]), (Subscript[σ, c]) ->  (Subscript[σ, t2]), (Subscript[σ, d]) ->  (Subscript[σ, t1])};
\(Print[
      StyleForm["a = ("\  <> \ 
          ToString[PaddedForm[aa, {3, 1}]]\  <> \ " ± "\  <> \ 
          ToString[
            PaddedForm[σ\_aa, {2, 
                1}]] <> \*"") cm/\!\(sec\^2\)""\ , 
        FontFamily\  -> \ "Helvetica", FontSize\  -> \ 12, 
        FontColor\  -> \ RGBColor[1, 0, 1]]]\)

Linear Regression

Suppose you measure two physical variables and they are connected by a linear relation of the form y = a + bx. If you were to measure N different values (Subscript[x, 1]),...,(Subscript[x, N]) and the corresponding values (Subscript[y, 1]),...,(Subscript[y, N]), you can find the straight line that fits the measurements; that is, to find the best estimates for the constants a and b based on the data. The Mathematica functions Fit and Regress enable you to do just that. The built-in function Fit finds a least-squares fit to the data as a linear combination of the specified basis functions. The functions Regress augments Fit by giving a list of commonly required diagnostics such as the coefficient of determination RSquared, the analysis of variance table ANOVATable, and the mean squared error EstimatedVariance. The output of regression functions can be controlled so that only needed information is produced. As an example, let's measure the temperature and pressure of an ideal gas kept at a constant volume. You have T = a + bP, where a is the temperature at which the pressure would drop to zero (if the gas did not condense into a liquid first); it is called the absolute zero of temperature. Let's find the absolute zero of temperature from the data by performing the least-squares fittings.

Let's first generate the input data and then plot the data.

Clear[pressure, temperature, dataset, 
        dataplot];
  
pressure = {65. , 75. , 85. , 95. , 105. };    
temperature = {-20. , 17. , 42. , 94. , 127. };   
dataset = 
        Transpose[{pressure, 
            temperature}];
\(dataplot = 
        ListPlot[dataset, Frame\  -> \ True, Axes\  -> \ False, 
          PlotRange\  -> \ {{60, 110}, {\(-30\), 
                130}}, \ 
		FrameLabel\  -> \ {"Pressure ( mm Hg)", \*""Temperature ( (^o)C)""}, 
		DefaultFont\  -> \ {"Helvetica", 11}, Prolog\  -> \ AbsolutePointSize[5]];\)

You will fit the data with the simple Fit function first. The intercept is fetched with datafit1[[1]] and the slope is fetched with datafit1[[2,1]].

Clear[x, datafit1, fitplot1];
datafit1 = Fit[dataset, {1, x}, x];
fitplot1 = 
        Plot[datafit1, {x, 0, 110}, 
          DisplayFunction -> Identity];
Print[
      StyleForm["T = " <> 
          ToString[
            PaddedForm[
              datafit1[[1]], {4, 
                2}]] <> " + " <> 
          ToString[
            PaddedForm[
              datafit1[[2, 1]], {3, 
                2}]] <> " P", FontFamily -> "Helvetica", 
        FontSize -> 12, FontColor -> RGBColor[1, 0, 1]]]; 
    Show[dataplot, fitplot1, 
      PlotRange  ->  {{0, 110}, {-300, 
            130}}];

Next, let's fit the data with the function Regress. Regress also gives a list of commonly required diagnostics such as fitted parameters, the variance, etc. The 6th command generates three lists, one for the fitted parameters, one for the standard deviations, and one for the confidence parameter intervals.

Clear[x, datafit2, parameter, sigma, ci, a, b, 
        fitplot2];
datafit2 = 
        Regress[dataset, {1, x}, x, 
          RegressionReport ->  {BestFit,  ParameterCITable, 
              FitResiduals, EstimatedVariance}];
residuals = FitResiduals /. datafit2;
(Subscript[σ, T])  = (Sqrt[EstimatedVariance /. datafit2]);
reducedChiSquare = 
        Apply[Plus, (((residuals/(Subscript[σ, T])))^2)]/(Length[residuals] - 
              2);
{parameter, sigma, ci} = 
        Transpose[(ParameterCITable /. 
                datafit2)[[1]]];
a = 
        parameter[[1]];
(Subscript[σ, a]) = 
        sigma[[1]];
b = 
        parameter[[2]];
(Subscript[σ, b]) = 
        sigma[[2]];
fitplot2 = 
        Plot[BestFit /. datafit2, {x, 0, 110}, 
          DisplayFunction -> Identity];
\(Print[
      StyleForm[\*""Reduced \!\(χ\^2\) = ""\  <> \ 
          ToString[PaddedForm[reducedChiSquare, {4, 2}]], 
        FontFamily\  -> \ "Helvetica", FontSize\  -> \ 12, 
        FontColor\  -> \ RGBColor[1, 0, 1]]]; 
    Print[StyleForm["T = ("\  <> \ 
          ToString[PaddedForm[a, {4, 2}]]\  <> \ " ± "\  <> \ 
          ToString[PaddedForm[σ\_a, {4, 2}]] <> ") + ("\  <> \ 
          ToString[PaddedForm[b, {3, 2}]]\  <> " ± "\  <> \ 
          ToString[PaddedForm[σ\_b, {4, 2}]]\ \  <> ") P", 
        FontFamily\  -> \ "Helvetica", FontSize\  -> \ 12, 
        FontColor\  -> \ RGBColor[1, 0, 1]]]; 
    Show[dataplot, fitplot2, 
      PlotRange\  -> \ {{0, 110}, {\(-300\), 130}}];\)

If the errors for dependent variable (T in this example) are known, they can be included in the linear regression by adding the option "Weights -> 1/(σ^2)" in Regress.

Clear[x, datafit2, parameter, sigma, ci, a, b, 
        fitplot2];
(Subscript[σ, T])  = {5.0, 6.0, 5.0, 6.0, 
          6.0};
datafit2 = 
        Regress[dataset, {1, x}, x, Weights ->  1/((Subscript[σ, T])^2), 
          RegressionReport ->  {BestFit,  ParameterCITable, FitResiduals, 
              EstimatedVariance}];
residuals = FitResiduals /. datafit2;
reducedChiSquare = 
        Apply[Plus, (((residuals/(Subscript[σ, T])))^2)]/(Length[residuals] - 
              2);
{parameter, sigmaRaw, ci} = 
        Transpose[(ParameterCITable /. 
                datafit2)[[1]]];
sigma = sigmaRaw/(Sqrt[reducedChiSquare]);
a = 
        parameter[[1]];
(Subscript[σ, a]) = 
        sigma[[1]];
b = 
        parameter[[2]];
(Subscript[σ, b]) = 
        sigma[[2]];
fitplot2 = 
        Plot[BestFit /. datafit2, {x, 0, 110}, 
          DisplayFunction -> Identity];
\(Print[
      StyleForm[\*""Reduced \!\(χ\^2\) = ""\  <> \ 
          ToString[PaddedForm[reducedChiSquare, {4, 2}]], 
        FontFamily\  -> \ "Helvetica", FontSize\  -> \ 12, 
        FontColor\  -> \ RGBColor[1, 0, 1]]]; 
    Print[StyleForm["T = ("\  <> \ 
          ToString[PaddedForm[a, {4, 2}]]\  <> \ " ± "\  <> \ 
          ToString[PaddedForm[σ\_a, {4, 2}]] <> ") + ("\  <> \ 
          ToString[PaddedForm[b, {3, 2}]]\  <> " ± "\  <> \ 
          ToString[PaddedForm[σ\_b, {4, 2}]]\ \  <> ") P", 
        FontFamily\  -> \ "Helvetica", FontSize\  -> \ 12, 
        FontColor\  -> \ RGBColor[1, 0, 1]]]; 
    Show[dataplot, fitplot2, 
      PlotRange\  -> \ {{0, 110}, {\(-300\), 130}}];\)

Nonlinear Regression

Sometimes the data can not be described by a linear function (functions that are linear in their parameters). In these cases you would have to use nonlinear functions such as y = a + b*(sin^2)(ω t + ϕ) where a, b, ω, and ϕ are the parameters, to describe the data. You can use two built-in Mathematica functions NonlinearFit and NonlinearRegress to perform the nonlinear fit. Before you do the nonlinear fit you will need to load the package Statistics`NonlinearFit` , as we have already done so in the beginning.

Let's first generate a data set with the function y = (sin^2)(π t) using the Mathematica function Table, and plot the data. A random noise with amplitude up to ± 0.05 is added to each data point.

Clear[t, oscdata, oscdataplot];
oscdata = 
        Table[Chop[N[{t, Sin[Pi*t]^2 + 0.1*(Random[] - 0.5)}]], {t, 0, 3, 
            0.05}];
\(\(oscdataplot = 
        ListPlot[oscdata, Frame\  -> \ True, Axes\  -> \ False, 
          PlotRange\  -> \ {{0, 3}, {\(-0.2\), 
                1.2}}, \ 
		FrameLabel\  -> \ {"t", "amplitude"}, 
		DefaultFont\  -> \ {"Helvetica", 11}, 
          Prolog\  -> \ 
            AbsolutePointSize[4]];\)\)

Let's fit the data to the nonlinear function y = a + b*(sin^2)(ω t + ϕ) with a, b, ω, and ϕ as parameters. In a nonlinear fit, not only should you provide the fit function (or model), you should also provide initial search values for the parameters. The function Timing is used to monitor how long it takes to do the fit. In some cases, it is useful to have this information. The fitted function is returned in BestFit which can be used to plot the fit result.

Clear[a, b, ω, ϕ, model, startingvalues, 
        fitresult];
model = 
        a + b*Sin[ω*t + ϕ]^2;
startingvalues = {{a, 0. }, {b, 1. }, {ω, Pi}, {ϕ, 
            0. }};
Timing[fitresult = 
          NonlinearRegress[oscdata, model, t, startingvalues, 
            RegressionReport  ->  {BestFit, ParameterCITable, 
                FitResiduals}];]

Next, let's extract the fitted parameters and their errors and compare them to the expected values.

Clear[parameter, sigma, ci];
{parameter, sigma, ci} = 
        Transpose[(ParameterCITable /. 
                fitresult)[[1]]];
afit = 
        parameter[[1]];
(Subscript[σ, afit]) = 
        sigma[[1]];
bfit = 
        parameter[[2]];
(Subscript[σ, bfit]) = 
        sigma[[2]];
ωfit = 
        parameter[[3]];
(Subscript[σ, ω])fit = 
        sigma[[3]];
ϕfit = 
        parameter[[4]];
(Subscript[σ, ϕ])fit = 
        sigma[[4]];
Print[
      StyleForm["a expected = 0;   a from fit = ("  <>  
          ToString[
            PaddedForm[afit, {5, 4}]]  <>  " ± "  <>  
          ToString[PaddedForm[(Subscript[σ, afit]), {5, 4}]] <> ")" , 
        FontFamily  ->  "Helvetica", FontSize  ->  12, 
        FontColor  ->  RGBColor[1, 0, 1]]]
Print[
      StyleForm["b expected = 1;   b from fit = ("  <>  
          ToString[
            PaddedForm[bfit, {5, 4}]]  <>  " ± "  <>  
          ToString[PaddedForm[(Subscript[σ, bfit]), {5, 4}]] <> ")" , 
        FontFamily  ->  "Helvetica", FontSize  ->  12, 
        FontColor  ->  RGBColor[1, 0, 1]]]
\(Print[
      StyleForm["ω expected = π;   ω from fit = ("\  <> \ ToString[
            PaddedForm[ωfit, {5, 
                4}]]\  <> \ " ± "\  <> \ 
          ToString[PaddedForm[σ\_ωfit, {5, 4}]] <> ")"\ , 
        FontFamily\  -> \ "Helvetica", FontSize\  -> \ 12, 
        FontColor\  -> \ RGBColor[1, 0, 1]]]\)
Print[
      StyleForm["ϕ expected = 0;   ϕ from fit = ("  <>  
          ToString[
            PaddedForm[ϕfit, {5, 4}]]  <>  " ± "  <>  
          ToString[PaddedForm[(Subscript[σ, ϕ])fit, {5, 4}]] <> ")" , 
        FontFamily  ->  "Helvetica", FontSize  ->  12, 
        FontColor  ->  RGBColor[1, 0, 1]]]

Finally, let's plot the fit together with the data.

\(\(oscfitplot = 
        Plot[\((BestFit /. fitresult)\), {t, 0, 3}, Frame\  -> \ True, 
          Axes\  -> \ False, 
          PlotRange\  -> \ {{0, 3}, {\(-0.2\), 
                1.2}}, \ 
		FrameLabel\  -> \ {"t", "amplitude"}, 
		DefaultFont\  -> \ {"Helvetica", 11}, 
          PlotStyle -> \ {RGBColor[1, 0, 1]}\ , 
          DisplayFunction -> Identity];\)\)
Show[oscdataplot, oscfitplot];

Let's plot the difference between data and the fit below.

time = Range[0, 3, 0.05];
ListPlot[Transpose[{time, (FitResiduals /. fitresult)}], 
        Frame -> True, Axes -> False, 
        PlotRange -> {{0, 3}, {-0.1, 0.1}}, 
        FrameLabel -> {"t", "Difference"}, 
        DefaultFont -> {"Helvetica", 11}, 
        Prolog ->  AbsolutePointSize[5]];

Fast Fourier Transform

A physical process can be described either in the time domain, by the values of some quantity h as a function of time t, i.e. h(t), or else in the frequency domain, where the process is specified by giving its amplitude H (generally a complex number indicating phase also) as a function of frequency ν, that is H(ν), with -∞ < ν < ∞. For many purposes it is useful to think of h(t) and H(ν) as being two different representations of the same function. One goes back and forth between these two representations by means of the Fourier transform equations,
		H(ν) = (1/(Sqrt[2 π]))\(TraditionalForm\`∫\_\(-∞\)\%∞\( h(
            t)\)\ \(e\^\(\(\ \)\(2  π\ i\ ν\ t\)\)\) dt\)
		  h(t) = (1/(Sqrt[2  π]))\(∫\_\(-∞\)\%∞\( 
                H(ν)\) \(e\^\(\(\ \)\(\(-2\) π\ i\ ν\ t\)\)\) dν\)
	If t is measured in seconds, then ν is in cycles per second, or Hertz. However, the equations work with other units. If h is a function of position (in meters), H will be a function of inverse wavelength (or k in cycles per meter), and so on. 
	The total power in a signal is the same whether we compute it in the time domain or in the frequency domain. This result is known as Parseval's therorem:
		Total Power = \(TraditionalForm\`∫\_\(-∞\)\%∞\(\(|\)\(\ \)\(h(t)\)\( | \^\(\(\ \)\(2\)\)\)\(dt\)\)\  = \ ∫\_\(-∞\)\%∞\(\(|\)\(\ \)\(H(ν)\)\( | \^\(\(\ \)\(2\)\)\)\(dν\)\)\)
Frequently one wants to know "how much power" is contained in the frequency interval between ν and ν + dν. In such circumstances one does not usually distinguish between positive and negative ν, but rather regards ν as varying from 0 (or D.C.) to +∞. In such cases, one defines the one-sided power spectral density (PSD) of the function h as
		\(\(\(\(P\_h\)(ν)\)\(\ \)\(≡\)\)\  | \ 
            H(ν)\( | \^\(\(\ \)\(2\)\)\)\ \(+\ \(\(|\)\(H(\(-ν\))\)\( | \^\(\(\ \)\(2\)\)\)\(\ \ \ \ \ \ \ \ \ \ \ \ \)\(0\  ≤ \ ν\  ≺ \ ∞\)\)\)\)
so that the total power is just the integral of (Subscript[P, h])(ν) from ν = 0 to ν = ∞.  Notice that if the function h(t) is real, H(ν) is still complex. In fact, in this case H(-ν) is the complex conjugate of H(ν); therefore, the two terms above (\(\(|\)\(H(ν)\)\( | \^2\)\) and \(\(|\)\(H(\(-ν\))\)\( | \^2\)\)) are equal.
	In the most common situations, function h(t) is sampled at evenly spaced intervals in time. This is the case for most of the data you will encounter in PHYS375. Let Δ denote the time interval between consecutive samples, so that the sequence of sampled values is
		(Subscript[h, n]) = h (nΔ)        n = ..., -3, -2, -1, 0, 1, 2, 3, ...
The reciprocal of the time interval Δ is called the sampling rate; if Δ is measured in seconds, for example, then the sampling rate is the number of samples recorded in second. For any sampling interval Δ, there is also a special frequency (Subscript[ν, c]), called the Nyquist critical frequency, given by
		(Subscript[ν, c]) ≡ (1/(2 Δ))
	Instead of evaluating the Fourier transform of the original continuous function h(t), we are now going to estimate the Fourier transform of the function from a finite number of its sampled points. Suppose that we have N consecutive sampled values:
		(Subscript[h, k]) ≡ h((Subscript[t, k])),      (Subscript[t, k]) ≡ kΔ,     k = 0, 1, 2, ..., N-1
With N numbers of input, we will evidently be able to produce no more than N independent numbers of output. So, instead of trying to estimate the Fourier transform H(ν) at all values of ν in the range -(Subscript[ν, c]) to (Subscript[ν, c]), let's seek estimates only at the discrete values:
		(Subscript[ν, n]) ≡ (n/(N Δ)),    n = -(N/2), ..., (N/2)
And the Fourier transform of the continuous function h(t) is now to be approximated by a discrete sum:
	H((Subscript[ν, n])) = (1/(Sqrt[2  π]))(∫_(-∞)^∞) h(
            t) (e^( 2 π i (Subscript[ν, n]) 
                t)) dt ≈ (1/(Sqrt[N]))(∑_(k = 0)^(N - 1))(Subscript[h, k])(e^( 2 π i (Subscript[ν, n]) (Subscript[t, k])))Δ 
				                 = (1/(Sqrt[N]))Δ (∑_(k = 0)^(N - 1))(Subscript[h, k])(e^( 2 π i n k /N))
The last summation in equation above is called the discrete Fourier transform of the N points (Subscript[h, k]). Let's denote it by (Subscript[H, n]),
		(Subscript[H, n]) ≡ (1/(Sqrt[N]))(∑_(k = 0)^(N - 1))(Subscript[h, k])(e^( 2 π i n k/N))
and
		H((Subscript[ν, n])) ≈ Δ (Subscript[H, n])
The discrete Fourier transform maps N complex numbers (the (Subscript[h, k])'s) into N complex numbers (the (Subscript[H, n])'s). It does not depend on any dimensional parameter, such as the time scale Δ. 
	The formula for the discrete inverse Fourier transform, which recovers the set of (Subscript[h, k])'s exactly from the (Subscript[H, n])'s is:
		(Subscript[h, k]) = (1/(Sqrt[N]))(∑_(n = 0)^(N - 1))(Subscript[H, n])(e^( -2 π i n k/N))
	Up until the mid-1960, it was thought that a direct calculation of a discrete Fourier transform would mean about (N^2) operations. A new algorithm called the Fast Fourier Transform, or FFT, was proposed and the discrete Fourier transform can then be computed in (N/2)(Subscript[log, 2])N operations. The difference between the two algorithms in terms of the computer CPU time is immense. For example, if N = 1024 (= (2^10)), the CPU time could be reduced by a factor of over 200! In Mathematica, the command Fourier[data] is used to perform the FFT.

Let's first look at a simple example with only eight elements for (Subscript[h, k]).

(Subscript[h, k]) = {4, 3, 7, -9, 1, 0, 0, 0};
TableForm[(Subscript[H, n]) = Fourier[(Subscript[h, k])]]

Notice that although the list (Subscript[h, k]) is real, (Subscript[H, n]) is comples. The first element of the transformed data is the zero frequency term which equal ∑(Subscript[h, k])/(Sqrt[N]) = 6/(Sqrt[8]) = 2.12132. The last three values of (Subscript[H, n]) correspond to negative frequencies and, for the real (Subscript[h, k]), they are complex conjugate of three (2nd - 4th) elements in the first half of (Subscript[H, n]).

Let's look at a more realistic example by generating a list of 200 elements containing two periodic signals (ν = 30 and 40) with random noise added. Notice that (Subscript[t, k]) = k/200, i.e., Δ = 1/200. This corresponds to a sampling rate of 200, or (Subscript[ν, c]) = (1/(2  Δ)) = 100. Notice that time for the last bin in (Subscript[h, k]) ((Subscript[h, 200])) is 1 second and (Subscript[ν, c]) is the frequency for the bin (Subscript[H, n]) at n = N/2 ((Subscript[H, (200/2)])).

data = 
        Table[N[Sin[30 2 π k/200] + 
              Sin[40 2 π k/200] + (Random[] - 1/2)], {k, 
            200}];
ndata = Length[data];

The data looks fairly random and complicated if you plot it directly.

ListPlot[data, PlotJoined -> True];

However, the power spectral density (\(TraditionalForm\`\(\(|\)\(H(ν)\)\( | \^2\)\)\)) shows two strong peaks at 30 and 40, reflecting the frequency components of the original signal at 30 and 40.

fourierTransform = 
        Drop[Abs[Fourier[data]]^2, 1];
fourierTransform = 
        Take[fourierTransform, ndata/2];
samplingPeriod = 1.0/200.0;
(Subscript[ν, c]) = 1. /(2 samplingPeriod);
frequency = 
        Range[1, ndata/2, 
            1] ((Subscript[ν, c])/(ndata/2));
fourierPlot = 
        ListPlot[Transpose[{frequency, fourierTransform}], 
          PlotJoined -> True, PlotRange -> All];

Pick, Copy and Paste Cursor Locations

Occasionally you may want to find the location or the height of a peak. This can be accomplished in Mathematica by using the mouse. To use the cursor, click once anywhere inside  the graph to select the graph. Then position the cursor at the top of the peak and simultaneously hold the CTRL key and left-click (holding it down for at least a few seconds) on the mouse. You can repeat the clicking of the left button if you want to make multiple measurements. The cursor locations can be copied into the paste buffer by either typing CTRL^C or left-clicking "File" and then "Copy". To paste the cursor locations just copied, position the cursor at the location where you want to paste the information and either type CTRL^V or left-click "File" and then "Paste".

Let's replot the graph for the Fourier transform. To determine the location and the height of the second peak, first click once anywhere inside the graph and postition the cursor at the top of the second peak. After copying the cursor locations to paste buffer, position the cursor after "cursorLocations =" and paste.

Show[fourierPlot];

cursorLocations = {{39.7938, 45.5499}, {39.7938, 45.5499}, {40.0581, 
            45.7448}, {39.7938, 45.1211}};

Two lists, one for the location and the other the height, can be generated by transposing cursorLocations. From these two lists, the average location and height of the second peak can then be determined.

{locations, heights} = Transpose[cursorLocations];
averageLocation = Mean[locations]
averageHeight = Mean[heights]