--- a
+++ b/code/preprocessing/lib/EOG_cut.m
@@ -0,0 +1,82 @@
+function [ EOG_segments, target_fs ] = EOG_cut( eog_signal, raw_fs, target_fs, raw_unit)
+%EOG_CUT Summary of this function goes here
+%   Detailed explanation goes here
+
+    %% 1. configurations
+   
+    
+    preprocess_filter_range = [0.3, 80]; %Hz
+    output_filter_range = [0.3, 12];
+    notch_frequency = 50;
+    filter_order = 22528;
+   
+    
+    cut_range = [-1,1]; %s output epochs range
+	
+    r_square_thresh = 0.6;
+    
+    
+    %% 2. preprocess
+    % 2.1 filter
+    eog_for_detection = filter_notch(eog_signal, raw_fs, notch_frequency); 
+    eog_for_detection=filter_data(eog_for_detection,raw_fs,preprocess_filter_range(1),preprocess_filter_range(2),filter_order);
+    % 2.2 resample
+    eog_for_detection = double(eog_for_detection);
+    eog_for_detection = resample(eog_for_detection', target_fs, raw_fs)';
+    eog_for_detection = double(detrend(eog_for_detection));
+    % 2.3 filter for output data if condigured
+    output_eog = filter_data(eog_for_detection, target_fs, [], output_filter_range,filter_order);
+	output_eog = double(output_eog);
+    
+
+    %% 3. cut segments
+    segment_len = (cut_range(2) - cut_range(1))*target_fs;
+    sample_num = length(eog_for_detection);
+    
+    segment_num = floor(sample_num/segment_len);
+    
+    
+    cut_locs = randi([target_fs+1, sample_num-target_fs], 1, segment_num);
+    
+    
+    segments_for_detection = nan(segment_num, segment_len);
+    EOG_segments = nan(segment_num, segment_len);
+    
+    g = fittype({'1/x','1'});
+    frq_interval=[2 79];%frq_interval=[2 48];
+    nfft = 1024;
+    df = [1 100];
+    for iter_segment = 1:segment_num
+        start_range = cut_locs(iter_segment) + cut_range(1)*target_fs;
+        stop_range = cut_locs(iter_segment) + cut_range(2)*target_fs;
+        range = start_range+1:stop_range;
+        
+        segments_for_detection(iter_segment, :) = eog_for_detection(range);
+        EOG_segments(iter_segment, :) = output_eog(range);
+        [freq, sp] = my_PSD(EOG_segments(iter_segment, :),nfft,target_fs,@hanning,50,df,'off');
+        vec = find(freq > frq_interval(1) & freq < frq_interval(2));
+        [x,f] = fit(freq(vec),sp(vec),g);
+        r_square(iter_segment) = f.rsquare;
+    end
+    
+    %% 4. preliminary exclusion based on 1/f fit
+    segment_to_remove = find(r_square<=r_square_thresh);
+    EOG_segments(segment_to_remove, :) = [];
+    
+    
+    %% 5. convert to uV
+    c = nan;
+    switch raw_unit
+        case 'uV'
+            c = 1;
+        case 'mV'
+            c = 1000;
+        case 'V'
+            c = 1000000;
+    end
+    if(isnan(c))
+        error('undefined input data unit!');
+    else
+        EOG_segments = EOG_segments.*c;
+    end
+end