%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is the implementation of Wavelet (WT) based Occular Artificat (OA)
% removal from single channel Electroencephalography (EEG) as described
% in the paper:
% 
%    S. Khatun, R. Mahajan, and B. I. Morshed, "Comparative Study of 
%    Wavelet Based Unsupervised Ocular Artifact Removal Techniques for 
%    Single Channel EEG Data�, IEEE Journal of Translational Engineering 
%    in Health and Medicine (JTEHM), vol. 4, no. 1, pp. 1-8, Dec. 2016.
%
% We request you to cite the above paper if you utilize this algorithm.
%
% Important Parameters:
%           sig    (input)    Input single channel EEG signals, 
%                             provided here as an input file "t.csv"
%           wname  (input)    Wavelet basis function utilized for cleaning
%           sig1   (output)   Reconstructed EEG signal (after OA removal)
%
% This version employs functions which are defined in some 
% other source code file.
%
% This work is licensed under a Creative Commons. 
% Attribution-NonCommercial 3.0 Unported License  
% ( http://creativecommons.org/licenses/by-nc/3.0/ )
%
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY EXPRESSED OR IMPLIED WARRANTIES
% OF ANY KIND, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THIS SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THIS SOFTWARE.
%
% BM 5/25/2016
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc; close all; clear all;

%% Examples of Necessary Parameters

sig = importdata('t.csv'); % Input raw EEG signal with ocular artifacts
wname = 'coif3'; % Wavelet basis function selection
Level = 9; % Levels of wavelet decomposition (this value corresponds to 
% sampling rate is 256 sps)

%% Wavelet Decomposition

[C L] = wavedec(sig,Level,wname);

%% Thresholding for OA removal
      
        j=L(1)+1;
        for i=1:1:7   
            x = C(j:1:(j+L(i+1)-1));  
            T = 1.5*std(x);             
            bb = 1;
            for k=j:1:(j+L(i+1)-1)
                  if(C(k)>T)
                        C(k)=0;
                    elseif(C(k)<-T)
                        C(k)=0;
                  end
            end
        j=j+L(i+1);
        end
        
%% Reconstruction and generating the clean EEG signal

sig1 = waverec(C,L,wname);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%