% %% plot CCA 
% clear all;
% load sbj1_sessitonA1;
% EEG = pop_select( EEG,'channel',{'C10' 'C12' 'C30' 'C31'});%right ear
% %EEG = pop_select( EEG,'channel',{'EXG2' 'E4' 'E8' 'E12'});%left ear
% %EEG = pop_select( EEG,'channel',{'A11' 'A12' 'A14' 'A15'});
% %EEG = pop_select( EEG,'channel',{'G15' 'G17' 'G20' 'G23'});%neck
% %EEG = pop_select( EEG,'channel',{'F26' 'F27' 'F28' 'F30'});%osccipital
% %EEG = pop_select( EEG,'channel',{'F29' 'F30' 'F31' 'F32'});%osccipital
% data = EEG.data;
% figure();
% for i=1:6 %for each epoche
%     result(:,i) = CCA(data(1,:,i),data(2,:,i),data(3,:,i),data(4,:,i),256);
% end
% 
% plot((32:55)/4,mean(result,2),'r','LineWidth',2);%plot mean

%% plot FFT 
clear all;
    EEG = pop_loadset('sbj4_sessitonA3.set');
    %%code below is used to analysis the data which cross sessions
    EEG2 = pop_loadset('sbj4_sessitonB4.set');
    EEG3 = pop_loadset('sbj4_sessitonC1.set');
    EEG4 = pop_loadset('sbj4_sessitonD10.set');
    
    %data_downSample = EEG.data(1:263*1024*6);
    %EEG.data = reshape(data_downSample,[263 1024 6]);
    
    
    EEG.data = EEG.data(:,:,1:6);
    EEG2.data = EEG2.data(:,:,1:6);
    EEG3.data = EEG3.data(:,:,1:6);
    EEG4.data = EEG4.data(:,:,1:6);
    EEG.trials = 6;
    for i=1:255
        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'];
%subject1
%badChannels = ['C24 ','E6 ','F1 ','G16 ','H25 ','H26 ','H27 ','H28 ','H29 ','H30 ','H31 ','H32 '];
%subject2
badChannels = ['A15 ']
specialChannelToPrint = ['C7 ' 'F24 ' 'E26 ' 'G13 '];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 = 1;
subplot_tag = plot_each_channelArea_figure;
channelName = 'A1';%in case didn't be assigned in channelArea
resolution = 4;% 4 = 0.25hz
freqFrom = 8;   
freqTo = 15;

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 = ([channelArea(channelAreaIndex),num2str(channelIndex)]);%'EXG1';%
        if(strfind(badChannels, [channelName,' '])~=0)
            ;%reach badchannel
        else
            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:',channelArea(channelAreaIndex)]);
                
                %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:',channelArea(channelAreaIndex)]);
                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;
                end
                xlabel('Hz');
                title(['FFT, Channel Area:',channelArea(channelAreaIndex)]);
                
                %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;
                end
                xlabel('Hz');
                ylabel('%');
                title(['SNR, Channel Area:',channelArea(channelAreaIndex)]);
                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
         if(strfind(specialChannelToPrint, [channelName,' '])~=0)
            specialChannel(specialChannelToPrintIndex,:) = snr_result(1,:);
            specialChannelToPrintIndex = specialChannelToPrintIndex+1;
        end
    end
    [B, index] = sort(maxSNRArray,'descend');
    fprintf('*******%c*******\r',channelArea(channelAreaIndex));
    fprintf('%i,',index);
    fprintf('\r****************\r');
    maxSNRArray(:,:) = 0;
    maxSNRChannelOverArea = strcat(maxSNRChannelOverArea,',', channelArea(channelAreaIndex), int2str(index(1)));
end
maxSNRChannelOverArea
figure();
%for(i=1:4)
    plot((x_axis-1)/resolution,specialChannel(3,:),'r',(x_axis-1)/resolution,specialChannel(1,:),'k',(x_axis-1)/resolution,specialChannel(2,:),'g',(x_axis-1)/resolution,specialChannel(4,:),'b');
    %hold on;
%end
xlabel('Hz');
ylabel('%');
title(['SNR']);
legend('occipital area', 'behind-ear area', 'face area', 'neck area');

%code below is for topoplot()
% figure();topoplot(snr_result_all,'location_test.loc','electrodes','on');
% if(FFT_first_then_getMean)%FFT
%     title([EEG.setname ', FFT first then getMean']);
% else
%     title([EEG.setname ', getMean first then FFT']);
% end
% colorbar();



