Quadrature Michelson Interferometer Data Analyzed in MATLAB

Below is an HTML version of the MATLAB Livescript version on how to analyze the data. The original (compressed) CSV data file is attached as the tek0016ALL.zip file as well as the modified MatLab data for Part 2 through 4 which is given as an attachment below in the file tek0016AllV1a.mat.

Click here for a PDF version of the MATLAB Livescript file. The actual script, QuadDataMoessbauer1a.mlx, can be uploaded below.

MATLAB Live Script HTML Version

KW 12/06/17

Overview:

With the new cage assembly from Thorlabs, we get such excellent data that we don't need the more sophisticated algorithm from Teachspin anymore.

The X vs. Y voltage data is so stable and circular, that we only have to subtract its mean (or median) and then use the atan2 and Matlab's phase unwrap function and we are done!

This script is broken into 4 sections.

The first section reads the times and the X and Y voltages from the scope's CSV file and then saves it as a MatLab data (mat) file. You should run this section only once to convert the data.

The second section loads the X and Y voltages and the time into MatLab.

The third section converts the MatLab data from voltage into pathlength by calculating and unwrapping the phase between the X and Y data.

The fourth section averages the pathlength data to reduce and smooth it and then calculates the corresponding velocity.

Part 1: Read the Data from Tektronix Scope's CSV File into Matlab Format and Save it Again

%NOTE: This for the Tektronix Scope CSV File Format; adjust for your own format

filename = 'tek0016AllV1a.csv';

delimiterIn = ',';

headerlinesIn = 15;

A = importdata(filename,delimiterIn,headerlinesIn);

sampleinterval = A.textdata(6) %display the sampling time (interval)

X = A.data(:,2);

Y = A.data(:,3);

dt = 2E-7; %sampling time interval: MAKE SURE THIS IS CORRECT and agrees with value displayed above

t1 = (1:size(X))*dt; %generate time vector

t = t1';

plot(t1, X, t1, Y)

grid on

title('Voltage for X and Y vs. Time')

xlabel('Time (seconds)')

ylabel('Voltage')

legend('X', 'Y')

plot(X, Y)

grid on

title('Vx vs. Vy')

xlabel('X')

ylabel('Y')

filename_mat = replace(filename,'csv','mat');

save( filename_mat, 'X', 'Y','t')

clear all

Part 2: Load the MatLab Formated X and Y Voltages and Times Into MatLab

filename = 'tek0016AllV1a.csv'; %use the original CSV file name; it will load .mat file

filename_mat = replace(filename,'csv','mat');

load( filename_mat);

plot(t, X, t, Y)

grid on

title('Voltage for X and Y vs. Time')

xlabel('Time (seconds)')

ylabel('Voltage')

legend('X', 'Y')

plot(t, X, t, Y)

grid on

title('Voltage for X and Y vs. Time Expanded')

xlabel('Time (seconds)')

ylabel('Voltage')

legend('X', 'Y')

xlim([0.0 0.001]) %Expanded Data range

Part 3: Convert the X Y Voltages into Distance

% calculate the mean or median and subtract it.

% probably could user better algorithm and we don't check for ellipticity

%tosc = (t>0.45E-4) & (t<1.25E-4);

tosc = (t>0.0); %starting time for mean or median calculations

%Xmean = mean(X(tosc))

Xmedian = median(X(tosc))

%Ymean = mean(Y(tosc))

Ymedian = median(Y(tosc))

%subtract the median

Xc = X-Xmedian;

Yc = Y-Ymedian;

plot(Xc, Yc)

grid on

title('Vx vs. Vy with Median Subtracted')

xlabel('X')

ylabel('Y')

% calculate phase relationship between X and Y and unwrap the phase

phiraw =atan2(Yc, Xc);

phirawunw = unwrap(phiraw);

%convert unwrapped phase to pathlength

lam = 633E-9; %HeNe laser wavelength

d = phirawunw*lam/(4*pi);

plot(t, d)

grid on

title('Distance vs. Time')

xlabel('Time (sec)')

ylabel('Distance (m)')

Part 4: Reduce and Smooth the Data Before Differentiating to Get Velocity

%Reduce the matrix size by averaging the data by 'comb' data points

comb = 100; %datapoints to aveage over

tcom = mean(reshape(t, comb, length(t)/comb));

dcom = mean(reshape(d, comb, length(t)/comb));

dt = tcom(2)-tcom(1); %calculate the time interval between the averaged data

%below is an alternative way to smooth the data; use at your own risk

%[dscom window] = smoothdata(dcom,'sgolay',100);

% window

dscom = dcom; % use the data obtained above by averaging data points

vcom=diff(dscom)/dt; %calculate the velocity

plot(tcom(1:length(vcom)), vcom)

title('Velocity vs. Time')

xlabel('Time (sec)')

ylabel('Velocity (m/sec)')

grid on