% Ultrafast Fibre Laser Experiment, University of Toronto.
% Written by: Gajanan Kuganesan Date: Spring 2004.
% Programmed for MATLAB 6.1, (available on Faraday via Nortel Labs).
%
% The following code imports the interferometric autocorrelation data and
% processed it to obtain the Full-Width-Half-Maximum (FHM) of the Intensity
% Autocorrelation. This is accomplished through a combination of both Mean
% and Median Filtering.
%
% Imported Data: Data can be imported from LABVIEW software (oscilloscope)
% available on the UFL computer. Data from LABVIEW should be saved as an
% Microsoft Excel file (*.xls), the setup data (first few lines) should be
% recorded and then elminated to be compatiable with the below code.
% That is the data imported should contain only 2 column in an *.xls file;
% the first column contains the oscilloscope voltage data and the second
% column contains the oscilloscope time data.

%*****************************************************************************

clear; %clear previous data
rawdata = importdata('c:\ultra1b\mes1.xls'); %Imports data as formatted above

ch1div = 20; %channel 1 volts/division, in milliVolts
ch1pos = -101.6; %channel 1 position w.r.t. rawdata zero
opwr = 1; %optical power of cable in milliWatts

% Arrange raw data into Voltage and Time vectors
for i = 1:length(rawdata)
voltage(i) = rawdata(i,1);
time(i) = rawdata(i,2);
end;

%*****************************************************************************

% Set up Vectors pased on Settings
voltage = (voltage*ch1div - ch1pos)/opwr; % Voltage Normalized by Optical Pwr
time = time*1E+3; % Time Set to Milliseconds
bkgd = mean(voltage(1:200)); % Calculated Average Background Signal
nvoltage = 2*voltage/bkgd; % Normalized Signal, Based on Background

amplitude = nvoltage; % Amplitude is analyzed signal in code
maxpeak = max(amplitude); % Maximum Peak Value
peakval = find(amplitude == max(amplitude));% Location of Max Peak (index value)

%*****************************************************************************

%FILTERING PROCESSS
fsize = 3; % Size of MEAN filter
fsize2 = 2; % Size of MEDIAN filter
N = 10; % Points Analyzed (about max value) for MEAN filter

data = sqrt(nvoltage); % Square-root to obtain Intensity
% (NOT Intensity Squared as measured)
avg = data; % Data to be Filtered.

% MEAN FILTER, about Peak Signal
for i = peakval-N:peakval+N
ctr = 0;
sum = 0;
for j = -fsize:fsize
if (((i+j)>1) & ((i+j)<length(data)))
ctr = ctr + 1;
sum = sum + data(i+j);
end;
if (ctr ~= 0)
avg(i) = sum/ctr;
end;
end;
end;

% MEDIAN FILTER
avg = medfilt1(avg, fsize2);

voltage1 = avg.^2; % Squared Signal to obtain Intensity Squared
fit = nvoltage - voltage1; % Data Fit, Original Data Minus Filtered Signal

%*****************************************************************************

% The following plot displays both the orignal signal and the filtered signal.
% It should be used as a check to confirm that the processed signal has a
% close (near identical) signal to that of the Intensity Autocorrelation.
figure; plot(time,nvoltage, '-r');
hold on; plot(time,voltage1, '-b');
title('Orignal Signal vs Processed Signal');
ylabel('Normalized Autocorrelation (Voltage)');
xlabel('Time Scale (Lag) [ms]');

%*****************************************************************************

% voltage1 is the signal with the background subtracted (for FWHM calculation)
% Background value calculated from mean value of initial signal.
voltage1 = voltage1 - mean(voltage1(1:200));

%*****************************************************************************

% The following plot displays the original signal, as well as the normalized
% data (normalized from the maximum observed signal).

%subplot(2,1,1); plot(time,voltage);
%title('Original Autocorrelation (Raw Data)');
%ylabel('Voltage (V)'); xlabel('Speaker Lag [ms]');

%subplot(2,1,2); plot(mtime,nvoltage);
%title('Calibrated Autocorrelation');
%ylabel('Normalized Voltage (to background)')
%xlabel('Time Scale (Lag) [ms]');

%*****************************************************************************

% The following plot displays the signal that has been filtered from the
% original data via Mean & Median Filtering
%figure; plot(time,fit);
%title ('Removed Signal');
%ylabel('Normalized Autocorrelation')
%xlabel('Time Scale (Lag) [ms]');

%*****************************************************************************

% The following plot displays the processed data due to Filtering
figure; plot(time, voltage1);
title('Final Processed Data');
ylabel('Normalized Autocorrelation')
xlabel('Time Scale (Lag) [ms]');

%*****************************************************************************

% Finds FWHM (not optimized for extreme spikes)
amplitude = voltage1;
maxpeak = max(amplitude);
peakindex = find(amplitude == max(amplitude));
peakrange = find((amplitude >= 0.5*max(amplitude)));
lowval = time(peakrange(1));
highval = time(peakrange(length(peakrange)));
FWHM = highval - lowval;
sig = voltage1;

%*****************************************************************************

% Displays the FWHM, Maximum Voltage, Median Filter and Mean Filter
% sizes on the previous plot.
text(time(100),0.5*max(sig),'FWHM:');
text(time(500),0.5*max(sig),num2str(FWHM));
text(time(100),0.75*max(sig),'Max Voltage:');
text(time(750),0.75*max(sig),num2str(max(voltage1)));
text(time(100),0.85*max(sig),'Median Filter Size:');
text(time(850),0.85*max(sig),num2str(fsize2));
text(time(100),0.95*max(sig),'Mean Filter Size:');
text(time(750),0.95*max(sig),num2str(fsize));

%*****************************************************************************

% The following plot displays the Power Spectrum (Frequency Domain)
% of the obtained data in dB. The large peaks at high frequencies
% (Niquist frequency) is evidence that the data is aliased.
%N=length(time);
%freq = -1/(2*0.1):1/(N*0.1):(N-1)/(2*N*0.1);
%ft = fftshift(fft(voltage));
%ps = ft.*conj(ft);
%figure; plot(freq,20*log10(ps/max(ps)));
%title('Power Spectrum'); ylabel('Power [dB]'); xlabel('Frequency [kHz]');