diff --git a/dev/matlab_files/b_norm.m b/dev/matlab_files/b_norm.m
new file mode 100644
index 00000000000..39e9f88f957
--- /dev/null
+++ b/dev/matlab_files/b_norm.m
@@ -0,0 +1,186 @@
+% b_norm 
+% normalize continuous data from BOXY (unparsed) - function version
+%
+% © 2015 University of Illinois Board of Trustees, All Rights Reserved.
+%
+%
+% elm 03/10/04 cleaned up print messages, supressed polyfit warnings
+%               & chgd "filt' to 'norm'
+% elm 2/24/11 added phase outlier cleaning
+% elm 7/7/11 saving boxy_hdr, good_norm_mrac_chs & good_mrdc_chs in output file
+
+%cycle	exmux	A-AC	A-DC	A-Ph	B-AC	B-DC	B-Ph	time	group	mark	flag
+
+function b_norm(log_file)
+global hdr file_name i_block dc ac ph adc aac aph msg mod_freq id_tags qual
+global boxy_hdr
+
+disp('b_norm 6.3');
+
+n_bad_points=12;    % number of anomalous initial data points to zap
+
+in_path=[hdr.data_path hdr.exp_name 'opt\'];
+
+[status,msg] = mkdir([hdr.data_path hdr.exp_name], ['norm00-00']);
+
+out_path=[hdr.data_path hdr.exp_name 'norm00-00\'];
+
+in_file_flg=fopen([in_path file_name]);
+out_file_flg=fopen([out_path file_name]);
+fclose('all');
+
+% Open the log file if it exists
+%
+if ~strcmp(log_file,'')
+    logid = fopen(log_file,'a');
+else logid = 1; end
+
+if ((in_file_flg>0) & (out_file_flg<0))
+    
+    % read a boxy file
+    %
+    [boxy_hdr,dc,ac,ph,aux]=read_boxy_hdr([in_path file_name]);
+    
+    % Fix bad first points (Dennis' problem :)
+    %
+    for index = 1:n_bad_points
+        dc(index,:)=dc(n_bad_points+1,:);
+        ac(index,:)=ac(n_bad_points+1,:);
+        ph(index,:)=ph(n_bad_points+1,:);
+    end
+    disp('opened file');
+	% fix phase wrap
+    % 1st estimate mean phase of point 1 to 50
+	% if a pont id more than 90 degrees different from the mean,
+    %  add or subtract 360 degrees
+    %
+	fprintf(logid,'Fixing phase wrap\n');
+    
+	for i_chan=1:boxy_hdr.n_chans
+        if mean(ph(1:50,i_chan)) < 180
+            wrapped_pts=find(ph(:,i_chan)>270);
+            ph(wrapped_pts,i_chan)=ph(wrapped_pts,i_chan)-360;
+        else
+            wrapped_pts=find(ph(:,i_chan)<90);
+            ph(wrapped_pts,i_chan)=ph(wrapped_pts,i_chan)+360;
+        end
+	end
+    
+    % Detrend (polynomial) phase data (fix Emily's problem ;-)
+    % (also sets mean = 0 !)
+    %
+    fprintf(logid,'Detrending\n');
+    y=1:boxy_hdr.n_points;
+    x=y';
+    warning off MATLAB:polyfit:RepeatedPointsOrRescale
+    
+    for i_ch=boxy_hdr.n_chans
+        poly_coeffs = polyfit(x,ph(:,i_ch),3);   % 3 => 3rd order
+        tmp_ph = ph(:,i_ch)-polyval(poly_coeffs,x);
+        ph(:,i_ch)=tmp_ph;
+    end
+    
+    % remove mean here so STD threshold makes sense
+    %
+    mrph=mean(ph);
+    for i_chan=1:boxy_hdr.n_chans
+        ph(:,i_chan)=(ph(:,i_chan)-mrph(i_chan));
+    end
+    
+    % remove phase (delay) outliers
+    %
+    hdr.qual.ph_out_thr=3; % always set to "3" per Kathy & Gabriele Oct 12 2012
+    SDph=std(ph(n_bad_points:end,:));
+    ph_out_thr=hdr.qual.ph_out_thr; 
+    for i_chan=1:boxy_hdr.n_chans
+        outliers=find(abs(ph(:,i_chan))>ph_out_thr*SDph(i_chan));
+        if length(outliers)>0
+            if outliers(1)==1;outliers=outliers(2:end);end % can't interp 1st pt
+            if outliers(end)==boxy_hdr.n_points % can't interp last pt
+                outliers=outliers(1:end-1);
+            end
+            n_ph_out(i_chan)=length(outliers);
+            for i_pt=1:n_ph_out(i_chan)
+                j_pt=outliers(i_pt);
+                ph(j_pt,i_chan)=(ph(j_pt-1,i_chan)+ph(j_pt+1,i_chan))/2;
+            end
+        end
+    end
+    
+	% Compute means
+	%
+	mrdc=mean(dc);
+%     dcneg=find(mrdc<0); % find chs w -dc
+%     mrdc(dcneg)=0; % and set them to 0
+	mrac=mean(ac);
+	mrph=mean(ph);
+    
+    % stash raw means in hdr
+    %
+    hdr.mrdc=mrdc;
+    hdr.mrac=mrac;
+    hdr.mrph=mrph;
+    
+    % apply DC & AC quality thresholds if desired
+    % zeroing 'bad' chans
+    %
+    good_norm_mrac_chs=[];
+    good_norm_mrdc_chs=[];
+    if qual.InitQualChks_flg==1
+        good_norm_mrac_chs=find(mrac>=hdr.qual.norm_mrac_thr);
+        bad_mrac_chs=find(mrac<hdr.qual.norm_mrac_thr);
+        good_norm_mrdc_chs=find(mrdc>=hdr.qual.norm_mrdc_thr);
+        bad_mrdc_chs=find(mrdc<hdr.qual.norm_mrdc_thr);
+
+        % set bad mrdc & mrac chans = 0
+        %
+        dc(:,bad_mrdc_chs)=0;
+        ac(:,bad_mrdc_chs)=0;
+        ph(:,bad_mrdc_chs)=0;
+        dc(:,bad_mrac_chs)=0;
+        ac(:,bad_mrac_chs)=0;
+        ph(:,bad_mrac_chs)=0;
+    else % use all chans
+        good_norm_mrac_chs=1:boxy_hdr.n_chans;
+        good_norm_mrdc_chs=1:boxy_hdr.n_chans;
+    end % if InitQualChks_flg
+    
+    % Normalize, and remove mean ('1') for filtering
+	%
+	fprintf(logid,'Rescaling (''Normalizing'')\n');
+	for i_chan=1:boxy_hdr.n_chans
+%         if mrdc(i_chan) ~= 0
+			dc(:,i_chan)=dc(:,i_chan)/mrdc(i_chan)-1;
+%         end
+%         if mrac(i_chan) ~= 0
+            ac(:,i_chan)=ac(:,i_chan)/mrac(i_chan)-1;
+%         end
+        ph(:,i_chan)=(ph(:,i_chan)-mrph(i_chan));
+    end
+    
+	ph=1e12*ph/(360*mod_freq); % convert phase to ps
+    
+    % save the cleaned up and normalized data
+    % could save outliers, would need to make cell array of channels...
+    %
+    hdr_norm=hdr; % save hdr WITH UNIQUE NAME!
+	msg=sprintf('Saving %s...\n',[out_path file_name]);
+    fprintf(logid,'Saving %s\n',[out_path file_name]);
+   	save([out_path file_name],'dc','ac','ph','mrdc','mrac','mrph','SDph',...
+        'boxy_hdr','aux',...
+        'poly_coeffs','id_tags','mod_freq','hdr_norm',...
+        'good_norm_mrac_chs','good_norm_mrdc_chs','-mat');
+% 'mark','boxy_flags','aux1','aux2','digaux','aux_data',
+
+else % if various files exist or don't...
+    if (in_file_flg<0)
+        msg=sprintf('%s does not exist.\n',[in_path file_name]);
+        fprintf(logid,'%s does not exist.\n',[in_path file_name]);
+    end
+    if (out_file_flg>0)
+        msg=sprintf('%s already exists!\n',[out_path file_name]);
+        fprintf(logid,'%s already exists!\n',[out_path file_name]);
+    end
+end
+
+if logid ~= 1; fclose(logid); end
diff --git a/dev/matlab_files/compare_python_matlab.m b/dev/matlab_files/compare_python_matlab.m
new file mode 100644
index 00000000000..63e8b7da1aa
--- /dev/null
+++ b/dev/matlab_files/compare_python_matlab.m
@@ -0,0 +1,104 @@
+% fprintf(logid,'Fixing phase wrap\n');
+
+filename = fullfile(pwd, 'dev' , 'matlab_files', 'data_matlab.mat')
+data = load(filename)
+data = data.data;
+
+% remove first few bad points
+
+n_bad_points = 12
+for index = 1:n_bad_points
+    data(:,index)=data(:,n_bad_points+1);
+end
+
+save(fullfile(pwd, 'dev' , 'matlab_files', 'data_first_bad_python.mat'),'data')
+
+% fix phase wrap
+
+for i_chan=1:size(data,1)
+    if mean(data(i_chan,1:50)) < 180
+        wrapped_pts=find(data(i_chan,:)>270);
+        data(i_chan,wrapped_pts)=data(i_chan,wrapped_pts)-360;
+    else
+        wrapped_pts=find(data(i_chan,:)<90);
+        data(i_chan,wrapped_pts)=data(i_chan,wrapped_pts)+360;
+    end
+end
+
+save(fullfile(pwd, 'dev' , 'matlab_files', 'data_unwrap_python.mat'),'data')
+
+
+% data = load("C:\Users\spork\Desktop\data_unwrap_matlab.mat")
+% data = data.data;
+
+% detrend data
+
+y=1:size(data,2);
+x=y;
+
+for i_chan=1:size(data,1)
+    poly_coeffs = polyfit(x,data(i_chan,:),3);   % 3 => 3rd order
+    tmp_ph = data(i_chan,:)-polyval(poly_coeffs,x);
+    data(i_chan,:)=tmp_ph;
+end
+
+save(fullfile(pwd, 'dev' , 'matlab_files', 'data_detrend_python.mat'),'data')
+
+% remove mean here so STD threshold makes sense
+%
+
+% data = load("C:\Users\spork\Desktop\data_detrend_matlab.mat")
+% data = data.data;
+
+mrph=mean(data,2);
+for i_chan=1:size(data,1)
+    data(i_chan,:)=(data(i_chan,:)-mrph(i_chan));
+end
+
+save(fullfile(pwd, 'dev' , 'matlab_files', 'data_mean_python.mat'),'data')
+
+
+% remove phase (delay) outliers
+
+% data = load("C:\Users\spork\Desktop\data_mean_matlab.mat")
+% data = data.data;
+
+ph_out_thr=3; % always set to "3" per Kathy & Gabriele Oct 12 2012
+SDph=std(data,[],2);
+
+for i_chan=1:size(data,1)
+    outliers=find(abs(data(i_chan,:))>ph_out_thr*SDph(i_chan));
+    if length(outliers)>0
+        if outliers(1)==1;
+            outliers=outliers(2:end);
+        end % can't interp 1st pt
+        if length(outliers)>0
+            if outliers(end)==size(data,2) % can't interp last pt
+                outliers=outliers(1:end-1);
+            end
+            n_ph_out(i_chan)=length(outliers);
+            for i_pt=1:n_ph_out(i_chan)
+                j_pt=outliers(i_pt);
+                data(i_chan, j_pt)=(data(i_chan,j_pt-1)+data(i_chan,j_pt+1))/2;
+            end
+        end
+    end
+end
+
+save(fullfile(pwd, 'dev' , 'matlab_files', 'data_outliers_python.mat'),'data')
+
+% normalise data
+mrph=mean(data,2);
+
+for i_chan=1:size(data,1)
+    data(i_chan,:)=(data(i_chan,:)-mrph(i_chan));
+end
+
+save(fullfile(pwd, 'dev' , 'matlab_files', 'data_norm_python.mat'),'data')
+
+% convert to picoseconds
+
+mod_freq = 1.1e8;
+data=1e12*data/(360*mod_freq); % convert phase to ps
+
+save(fullfile(pwd, 'dev' , 'matlab_files', 'data_picosec_python.mat'),'data')
diff --git a/dev/matlab_files/data_detrend_python.mat b/dev/matlab_files/data_detrend_python.mat
new file mode 100644
index 00000000000..b10ef2fc00a
Binary files /dev/null and b/dev/matlab_files/data_detrend_python.mat differ
diff --git a/dev/matlab_files/data_first_bad_python.mat b/dev/matlab_files/data_first_bad_python.mat
new file mode 100644
index 00000000000..8fae1859ad4
Binary files /dev/null and b/dev/matlab_files/data_first_bad_python.mat differ
diff --git a/dev/matlab_files/data_matlab.mat b/dev/matlab_files/data_matlab.mat
new file mode 100644
index 00000000000..40e8182fc80
Binary files /dev/null and b/dev/matlab_files/data_matlab.mat differ
diff --git a/dev/matlab_files/data_mean_python.mat b/dev/matlab_files/data_mean_python.mat
new file mode 100644
index 00000000000..34c6c75ac0c
Binary files /dev/null and b/dev/matlab_files/data_mean_python.mat differ
diff --git a/dev/matlab_files/data_norm_python.mat b/dev/matlab_files/data_norm_python.mat
new file mode 100644
index 00000000000..7283400ef95
Binary files /dev/null and b/dev/matlab_files/data_norm_python.mat differ
diff --git a/dev/matlab_files/data_outliers_python.mat b/dev/matlab_files/data_outliers_python.mat
new file mode 100644
index 00000000000..6776ee664ad
Binary files /dev/null and b/dev/matlab_files/data_outliers_python.mat differ
diff --git a/dev/matlab_files/data_picosec_python.mat b/dev/matlab_files/data_picosec_python.mat
new file mode 100644
index 00000000000..4ac30bc33d3
Binary files /dev/null and b/dev/matlab_files/data_picosec_python.mat differ
diff --git a/dev/matlab_files/data_unwrap_python.mat b/dev/matlab_files/data_unwrap_python.mat
new file mode 100644
index 00000000000..9544dce13ef
Binary files /dev/null and b/dev/matlab_files/data_unwrap_python.mat differ
diff --git a/dev/matlab_files/phase_wrap_data.pickle b/dev/matlab_files/phase_wrap_data.pickle
new file mode 100644
index 00000000000..c8ea5548dd2
Binary files /dev/null and b/dev/matlab_files/phase_wrap_data.pickle differ
diff --git a/dev/test_phase_wrap.py b/dev/test_phase_wrap.py
new file mode 100644
index 00000000000..75a59672192
--- /dev/null
+++ b/dev/test_phase_wrap.py
@@ -0,0 +1,161 @@
+import os
+import pickle
+import scipy
+from scipy import io
+import numpy as np
+
+crnt_dir = os.getcwd()
+
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 
+                        'phase_wrap_data.pickle')
+
+# filename_mat = os.path.join(crnt_dir, 'dev' , 'matlab_files', 
+#                             'all_datamatlab.mat')
+
+# with open(filename, 'wb') as f:
+#     pickle.dump(all_data, f)
+    
+# scipy.io.savemat(filename_mat, mdict=dict(data=all_data))
+	
+with open(filename, 'rb') as f:
+    all_data = pickle.load(f)
+    
+###phase unwrapping###
+thresh = 1e-10 #determine to which decimal place we will compare
+
+# print('Fixing first few bad points')
+
+n_bad_points = 12
+for i_point in range(n_bad_points):
+    all_data[:,i_point] = all_data[:,n_bad_points] 
+    
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 'data_first_bad_python.mat')
+first_bad_data = scipy.io.loadmat(filename)
+
+test6 = abs(first_bad_data['data'] - all_data) <= thresh
+if test6.all():
+    print("SUCCESS! First bad matches Matlab to " + str(thresh) + " decimal points!")
+else:
+    print("FAILURE! First bad DOES NOT match Matlab to " + str(thresh) + " decimal points!")
+
+# print('Fixing phase wrap')
+
+for i_chan in range(np.size(all_data, axis=0)):
+    if np.mean(all_data[i_chan,:50]) < 180:
+        wrapped_points = all_data[i_chan, :] > 270
+        all_data[i_chan, wrapped_points] -= 360
+    else:
+        wrapped_points = all_data[i_chan,:] < 90
+        all_data[i_chan, wrapped_points] += 360
+      
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 'data_unwrap_python.mat')
+unwrapped_data = scipy.io.loadmat(filename)
+
+test1 = abs(unwrapped_data['data'] - all_data) <= thresh
+if test1.all():
+    print("SUCCESS! Phase wrap matches Matlab to " + str(thresh) + " decimal points!")
+else:
+    print("FAILURE! Phase wrap DOES NOT match Matlab to " + str(thresh) + " decimal points!")
+        
+# print('Detrending phase data')
+# scipy.io.savemat(file_name = r"C:\Users\spork\Desktop\all_dataunwrap_matlab.mat",
+#               mdict=dict(data=all_data))
+
+y = np.linspace(1, np.size(all_data, axis=1),
+                np.size(all_data, axis=1))
+x = np.transpose(y)
+for i_chan in range(np.size(all_data, axis=0)):
+    poly_coeffs = np.polyfit(x,all_data[i_chan, :] ,3)
+    tmp_ph = all_data[i_chan, :] - np.polyval(poly_coeffs,x)
+    all_data[i_chan, :] = tmp_ph
+    
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 'data_detrend_python.mat')
+detrend_data = scipy.io.loadmat(filename)
+    
+test2 = abs(detrend_data['data'] - all_data) <= thresh
+if test2.all():
+    print("SUCCESS! Detrend matches Matlab to " + str(thresh) + " decimal points!")
+else:
+    print("FAILURE! Detrend DOES NOT match Matlab to " + str(thresh) + " decimal points!")
+
+# print('Removing phase mean')
+# scipy.io.savemat(file_name = r"C:\Users\spork\Desktop\all_datadetrend_matlab.mat",
+#               mdict=dict(data=all_data))
+
+mrph = np.mean(all_data,axis=1);
+for i_chan in range(np.size(all_data, axis=0)):
+    all_data[i_chan,:]=(all_data[i_chan,:]-mrph[i_chan])
+    
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 'data_mean_python.mat')
+mean_data = scipy.io.loadmat(filename)
+
+test3 = abs(mean_data['data'] - all_data) <= thresh
+if test3.all():
+    print("SUCCESS! Subtracting mean matches Matlab to " + str(thresh) + " decimal points!")
+else:
+    print("FAILURE! Subtracting mean DOES NOT match Matlab to " + str(thresh) + " decimal points!")
+
+# print('Removing phase outliers')
+# scipy.io.savemat(file_name = r"C:\Users\spork\Desktop\all_datamean_matlab.mat",
+#               mdict=dict(data=all_data))
+
+ph_out_thr=3;
+sdph=np.std(all_data,1, ddof = 1); #set ddof to 1 to mimic matlab
+n_ph_out = np.zeros(np.size(all_data, axis=0), dtype= np.int8)
+
+for i_chan in range(np.size(all_data, axis=0)):
+    outliers = np.where(np.abs(all_data[i_chan,:]) > 
+                (ph_out_thr*sdph[i_chan]))
+    outliers = outliers[0]
+    if len(outliers) > 0:
+        if outliers[0] == 0:
+            outliers = outliers[1:]
+        if len(outliers) > 0:
+            if outliers[-1] == np.size(all_data, axis=1) - 1:
+                outliers = outliers[:-1]
+            n_ph_out[i_chan] = int(len(outliers))
+            for i_pt in range(n_ph_out[i_chan]):
+                j_pt = outliers[i_pt]
+                all_data[i_chan,j_pt] = (
+                    (all_data[i_chan,j_pt-1] + 
+                      all_data[i_chan,j_pt+1])/2)
+                
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 'data_outliers_python.mat')
+outlier_data = scipy.io.loadmat(filename)
+
+test4 = abs(outlier_data['data'] - all_data) <= thresh
+if test4.all():
+    print("SUCCESS! Removed outliers match Matlab to " + str(thresh) + " decimal points!")
+else:
+    print("FAILURE! Removed outliers DO NOT match Matlab to " + str(thresh) + " decimal points!")
+    
+    
+# now let's normalise our data #
+for i_chan in range(len(all_data)):
+    all_data[i_chan,:] = (all_data[i_chan,:] - np.mean(all_data[i_chan,:]))
+    
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 'data_norm_python.mat')
+outlier_data = scipy.io.loadmat(filename)
+
+test7 = abs(outlier_data['data'] - all_data) <= thresh
+if test4.all():
+    print("SUCCESS! Norm data matches Matlab to " + str(thresh) + " decimal points!")
+else:
+    print("FAILURE! Norm data DOES NOT match Matlab to " + str(thresh) + " decimal points!")
+
+
+#convert phase to pico seconds
+# scipy.io.savemat(file_name = r"C:\Users\spork\Desktop\data_picosec_matlab.mat",
+#               mdict=dict(data=all_data))
+mtg_mdf = 1.1e8
+for i_chan in range(np.size(all_data, axis=0)):
+    all_data[i_chan,:] = ((1e12*all_data[i_chan,:])/(360*mtg_mdf))
+    
+filename = os.path.join(crnt_dir, 'dev' , 'matlab_files', 'data_picosec_python.mat')
+picosec_data = scipy.io.loadmat(filename)
+
+test5 = abs(picosec_data['data'] - all_data) <= thresh
+if test5.all():
+    print("SUCCESS! Pico seconds match Matlab to " + str(thresh) + " decimal points!")
+else:
+    print("FAILURE! Pico seconds DO NOT match Matlab to " + str(thresh) + " decimal points!")
\ No newline at end of file