%  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]');