
%% =======NOTE========
%   analysis from A1->EXG8, skip EXG4(reference)
%%
clear all;
    %75-9Hz
%     EEG = pop_loadset('filename','sbj4_sessitonA10.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG2 = pop_loadset('filename','sbj4_sessitonB2.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG3 = pop_loadset('filename','sbj4_sessitonC3.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG4 = pop_loadset('filename','sbj4_sessitonD5.set','filepath','/data/projects/yute/SSVEP/subject4/');
    %75-10Hz
%     EEG = pop_loadset('filename','sbj4_sessitonA3.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG2 = pop_loadset('filename','sbj4_sessitonB4.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG3 = pop_loadset('filename','sbj4_sessitonC1.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG4 = pop_loadset('filename','sbj4_sessitonD10.set','filepath','/data/projects/yute/SSVEP/subject4/');
     %75-11Hz
%     EEG = pop_loadset('filename','sbj4_sessitonA7.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG2 = pop_loadset('filename','sbj4_sessitonB5.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG3 = pop_loadset('filename','sbj4_sessitonC9.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG4 = pop_loadset('filename','sbj4_sessitonD6.set','filepath','/data/projects/yute/SSVEP/subject4/');
%          %75-12Hz
%     EEG = pop_loadset('filename','sbj4_sessitonA4.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG2 = pop_loadset('filename','sbj4_sessitonB10.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG3 = pop_loadset('filename','sbj4_sessitonC5.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG4 = pop_loadset('filename','sbj4_sessitonD1.set','filepath','/data/projects/yute/SSVEP/subject4/');
         %75-13Hz
%     EEG = pop_loadset('filename','sbj4_sessitonA2.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG2 = pop_loadset('filename','sbj4_sessitonB3.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG3 = pop_loadset('filename','sbj4_sessitonC4.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG4 = pop_loadset('filename','sbj4_sessitonD2.set','filepath','/data/projects/yute/SSVEP/subject4/');
    %120-9hz
%     EEG = pop_loadset('filename','sbj4_sessitonA8.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG2 = pop_loadset('filename','sbj4_sessitonB6.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG3 = pop_loadset('filename','sbj4_sessitonC10.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG4 = pop_loadset('filename','sbj4_sessitonD7.set','filepath','/data/projects/yute/SSVEP/subject4/');
    %120-10hz
    EEG = pop_loadset('filename','sbj4_sessitonA1.set','filepath','/data/projects/yute/SSVEP/subject4/');
    EEG2 = pop_loadset('filename','sbj4_sessitonB8.set','filepath','/data/projects/yute/SSVEP/subject4/');
    EEG3 = pop_loadset('filename','sbj4_sessitonC6.set','filepath','/data/projects/yute/SSVEP/subject4/');
    EEG4 = pop_loadset('filename','sbj4_sessitonD4.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     %120-13hz
%     EEG = pop_loadset('filename','sbj4_sessitonA5.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG2 = pop_loadset('filename','sbj4_sessitonB9.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG3 = pop_loadset('filename','sbj4_sessitonC7.set','filepath','/data/projects/yute/SSVEP/subject4/');
%     EEG4 = pop_loadset('filename','sbj4_sessitonD9.set','filepath','/data/projects/yute/SSVEP/subject4/');
    
    EEG.trials =6;
    EEG.data = EEG.data(:,:,1:EEG.trials);
    EEG2.data = EEG2.data(:,:,1:EEG.trials);
    EEG3.data = EEG3.data(:,:,1:EEG.trials);
    EEG4.data = EEG4.data(:,:,1:EEG.trials);
    for i=1:EEG.nbchan
        meanFour (i,:,:) = (EEG.data(i,:,:)+EEG2.data(i,:,:)+EEG3.data(i,:,:)+EEG4.data(i,:,:))/4;
    end
    EEG.data = meanFour;
    
    
%%
channelArea = {'A' 'B' 'C' 'D' 'E' 'F' 'G' 'H' 'EXT'};
badChannels = {'H25','H26','H27','H28','H29','H30','H31','H32'};%['C24 ','E6 ','F1 ','G16 ','H25 ','H26 ','H27 ','H28 ','H29 ','H30 ','H31 ','H32 '];
specialChannelToPrint = {'C30' 'EXT1' 'F25' 'G29'};
%rank: face->ear->occipital->neck

specialChannelToPrintIndex =1;
maxSNRChannelOverArea = ''; %be used in the end
SNR_hz_we_want = 10; 
FFT_first_then_getMean = 0;%0:getMean first /  1:getFFT first
plot_each_channelArea_figure = 0;
subplot_tag = plot_each_channelArea_figure;
channelName = 'A1';%in case didn't be assigned in channelArea
resolution = 4;% 4 = 0.25hz
freqFrom = 8;   
freqTo = 25;

maxSNRArray_index=1;%initial
EEGOri = EEG;
x_axis = (freqFrom*resolution+1:freqTo*resolution+1);
%check error
if(SNR_hz_we_want<freqFrom || SNR_hz_we_want>freqTo)
    fprintf('*****Error, SNR_hz_we_want<freqFrom || SNR_hz_we_want>freqTo****\r\r');
end;
SNR_hz_we_want = SNR_hz_we_want*resolution+1;
for channelAreaIndex=1:length(channelArea)% A->H
    if(plot_each_channelArea_figure)
        figure();
    end
    for channelIndex=1:32
        channelName = strcat(channelArea(channelAreaIndex),num2str(channelIndex));%([channelArea(channelAreaIndex),num2str(channelIndex)]);%'EXG1';%
        if(strmatch(badChannels,channelName)~=0)%if(strfind(badChannels, [channelName,' '])~=0)
            ;%reach badchannel
        else
            if(strmatch('EXT9',channelName)~=0)
                break;
            end
            if(strmatch('EXT4',channelName)~=0)
                continue;
            end
            EEG = pop_select( EEGOri,'channel',channelName);
            data = EEG.data;
            %%  This block is get the FFT of each epoches, then get mean, then do SNR operation.          
            %plot FFT, first do fft in each epoche then get the mean.
            if(FFT_first_then_getMean)
                for i=1:EEG.trials%for each epoche
                    fft_result(i,:) = abs(fft(data(1,:,i)));
                    if(subplot_tag)
                        subplot(1,2,1);plot((x_axis-1)/resolution,fft_result(i,x_axis));%plot FFT
                        hold on;
                        %pause;
                    end
                end
                subplot(1,2,1);plot((x_axis-1)/resolution,mean(fft_result(:,x_axis)),'r','LineWidth',2);
                xlabel('Hz');
                title(['FFT, Channel Area:',channelName]);
                
                %plot SNR
                for i=1:EEG.trials%for each epoche
                    snr_sum = sum(fft_result(i,x_axis));
                    snr_result(i,:) = fft_result(i,x_axis)/snr_sum*100;
                    if(subplot_tag)
                        subplot(1,2,2);plot((x_axis-1)/resolution,snr_result(i,:));
                        hold on;
                        %pause;
                    end;
                end
                subplot(1,2,2);plot((x_axis-1)/resolution,mean(snr_result(:,:)),'r','LineWidth',2);
                xlabel('Hz');
                ylabel('%');
                title(['SNR, Channel Area:',channelName]);
                maxSNRArray(channelIndex)=max(mean(snr_result(:,:)));
                %pull out the SNR on 'SNR_hz_we_want'
                temp = mean(snr_result(:,:));
                snr_result_all(maxSNRArray_index,:) = temp(SNR_hz_we_want-(freqFrom*resolution));
                
                maxSNRArray_index=maxSNRArray_index+1;
            else
                %% The block below is to get the mean over epoches, then FFT or SNR operation.
                %plot FFT, first get the mean over the epoches, then do FFT
                data_temp =  zeros([1 length(data)]);%initial the data_temp
                for i=1:EEG.trials%for each epoche
                    data_temp(1,:) = data_temp(1,:)+data(1,:,i);%sum each epoche
                end
                data_temp(1,:) = data_temp(1,:)/EEG.trials;%get average
                fft_result(1,:) = abs(fft(data_temp(1,:)));
                if(subplot_tag)
                    subplot(1,2,1);plot((x_axis-1)/resolution,fft_result(1,x_axis));%plot FFT
                    hold on;
                    %pause;
                    xlabel('Hz');
                 title(['FFT, Channel Area:',channelName]);
                end
                
                
                %plot SNR= SNR_hz_we_want/mean(freqFrom->freqTo, excep SNR_hz_we_want)
                snr_sum = sum(fft_result(1,x_axis))-fft_result(1,x_axis(SNR_hz_we_want-(freqFrom*resolution)));
                snr_result(1,:) = fft_result(1,x_axis)/(snr_sum/(length(x_axis)-1));           
%                 %plot SNR=SNR_hz_we_want/(freqFrom->freqTo)  
%                 snr_sum = sum(fft_result(1,x_axis));
%                 snr_result(1,:) = fft_result(1,x_axis)/snr_sum*100;

                if(subplot_tag)
                    subplot(1,2,2);plot((x_axis-1)/resolution,snr_result(1,:));
                    hold on;
                    %pause;
                    xlabel('Hz');
                    ylabel('%');
                    title(['SNR, Channel Area:',channelName]);
                end
                
                maxSNRArray(channelIndex)=max(snr_result(1,:));
                %pull out the SNR on 'SNR_hz_we_want'
                snr_result_all(maxSNRArray_index) = snr_result(1,SNR_hz_we_want-(freqFrom*resolution));
                
                maxSNRArray_index=maxSNRArray_index+1;
            end
        end
         specialChannelToPrintIndex=0;
         specialChannelToPrintIndex = strmatch(channelName,specialChannelToPrint,'exact');
         if(specialChannelToPrintIndex<=4)
             specialChannel(specialChannelToPrintIndex,:) = snr_result(1,:);
             %specialChannelToPrintIndex = specialChannelToPrintIndex+1;
         end
    end
    [B, index] = sort(maxSNRArray,'descend');
    fprintf('*******%c*******\r',cell2mat(channelName));
    fprintf('%i,',index);
    fprintf('\r****************\r');
    maxSNRArray(:,:) = 0;
    %maxSNRChannelOverArea = strcat(maxSNRChannelOverArea,',', channelArea(channelAreaIndex), int2str(index(1)));
end
%maxSNRChannelOverArea
figure();
%for(i=1:4)
figure_SNR=plot((x_axis-1)/resolution,specialChannel(1,:),'g',(x_axis-1)/resolution,specialChannel(2,:),'k',(x_axis-1)/resolution,specialChannel(3,:),'r',(x_axis-1)/resolution,specialChannel(4,:),'b');
%hold on;
%end
title([EEG.setname, ' ',num2str((SNR_hz_we_want-1)/resolution),'Hz',' Ch= ' specialChannelToPrint{1},' ', specialChannelToPrint{2},' ', specialChannelToPrint{3},' ', specialChannelToPrint{4}]);%title(['SNR, Ch=',specialChannelToPrint]);
legend_SNR=legend('face area','behind-ear area','occipital area','neck area');
xlabel('Hz','fontsize',12);
ylabel('%','fontsize',12);
set(figure_SNR,'LineWidth',2);
set(legend_SNR,'FontSize',12);
set(gca,'FontSize',12)
box off
% code below is for topoplot()
figure();
topoplot(snr_result_all,'location_topoplot.sfp');%'/data/projects/yute/SSVEP/location_test.loc','electrodes','on');
title([EEG.setname,' ',num2str((SNR_hz_we_want-1)/resolution),'Hz', ' Ch= ' specialChannelToPrint{1},' ', specialChannelToPrint{2},' ', specialChannelToPrint{3},' ', specialChannelToPrint{4}]);%title(strcat(EEG.setname,',Ch=',specialChannelToPrint));%title([EEG.setname ',Ch=',specialChannelToPrint]);

colorbar();
save data10_120hz specialChannel;