function [sig_cell,sig_cell_tot] = FH_sig_tests(srce,trigdirname,sig_test,p_indx,latency,stimlist);
% [sig_cell,sig_cell_tot] = FH_sig_tests(srce,trigdirname,sig_test,p_indx,latency,stimlist)
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Significance test for F-test and Hotelling test
% For a specified sqd file and its corresponding trigger files, List of frequencies
% for all stimuli, and false positive percentage, we compute
% significance test based on F-test and/or the Hotelling test.  
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% srce : sqd file (if not specified, a sqd file is requested)
% trigdirname : path for trigger files that contains all 20 trigger files
% named in ascending order trig_160 to trig_179
% sig_test : 3 conditions 1: F-test, 2: Hotelling test, 3: both F and Hotelling test
% p_indx : Number of false positives on average per 157 channels
% latency : Post trigger delay for Steady State Response
% stimlist : List of stimulus frequencies 
%
% sig_cell : cell of significant channels for all trigger files
% sig_cell_tot : Confusion matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Author: Nayef Ahmar aenayef@glue.umd.edu
%
% for the Computational Sensorimotor Systems Lab (CSSL) UMCP
% http://www.isr.umd.edu/Labs/CSSL/
% Version 0.91 Dec. 1 2005
%
% Latest version at
% http://www.isr.umd.edu/Labs/CSSL/simonlab/resources/resources.html
%

if nargin == 0,
    [fn, pn] = uigetfile('*.sqd','Select your SQD file source');
    if sum(fn==0)&&sum(pn==0),return,end
    srce = fullfile(pn,fn);
end
if nargin < 6, stimlist = [1.5 1.5 1.5 1.5 3.5 3.5 3.5 3.5 7.5 7.5 7.5 7.5 15.5 15.5 15.5 15.5 31.5 31.5 31.5 31.5];end ;
if nargin < 5, latency = 0.3;end ;
if nargin < 4, p_indx = 1;end ;
if nargin < 3, sig_test = 1;end ;
if nargin < 2, trigdirname = uigetdir( 'Pick a Directory for trigger files'); end;

info = sqdread(srce,'info');%extract info
sf = info.SampleRate;%sampling frequency in HZ
t0 = clock;
num_stim = length(stimlist);
freqlist = unique(stimlist);
numfreq = length(freqlist);
num_false_pve = p_indx*16;
test_trig_sam_cube = zeros(2,num_stim,numfreq,157); %@initialize cube of all data
inc = 1;%matlab start at zero hence add one
pretrigger = 0;
posttrigger = 0;
stim_duration = 2;%stimulus duration in seconds
len = length(textread([trigdirname,'/trig_160.txt']));
epoc_samples = stim_duration*sf;
channel = 1:157;
fftw('planner','patient');%Choose fastest FFT for
h = waitbar(0,'Please wait...');

for trigch = 160:179,%[160:175,178:181],
    triggerfile = [trigdirname,'/trig_',int2str(trigch),'.txt'];
    trigger = textread(triggerfile);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %Extract and concatenate sample paths of filtered data for all channels
    %maxtrig =  trigger(len)*sf + latency*sf + epoc_samples + posttrigger*sf ;
    %data_stim = zeros(1:maxtrig,157);
    data_stim = [];
    for j = 1:len,
        lowtrig = round(trigger(j)*sf - pretrigger*sf + latency*sf + 1) ;%trigger -(pretrigger) + 1(matlab starts at 1 vs meg160 at 0)
        hightrig = round(lowtrig + pretrigger*sf + epoc_samples + posttrigger*sf -1);%  (pretrigger) + 1s recorded data +  (posttrigger)
        datain = getdata(info,'Channels',0:156,'Samples',[lowtrig hightrig]);
        data_stim = [data_stim; datain];
    end
    for sami = 1:numfreq,
        sam = freqlist(sami);%[1.5,3.5,7.5,15.5,31.5],
        amf = sam*stim_duration;
        am_freq = amf*len + inc;

        for chan = channel,
            data_ch = data_stim(:,chan);

            if sig_test ~= 2,
                fft_data_ch = fft(data_ch);
                fft_data_ch(1) = 0;
                fft_data_amf(chan) = fft_data_ch(am_freq,1);

                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %F test for hidden periodicity(Simplified)
                denom_f = sum(abs(fft_data_ch(am_freq-60:am_freq-1)).^2)+sum(abs(fft_data_ch(am_freq+1:am_freq+60)).^2);
                num_f = 120*abs(fft_data_ch(am_freq))^2;
                R_f_test = num_f/denom_f;
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                R_f_test_chlist(chan,:) = R_f_test;
            end

            if sig_test ~= 1,
                data_ch = data_ch - mean(data_ch);%
                %Hotelling
                for k = 1:len,
                    lb = (k-1)*epoc_samples + 1;
                    ub = k*epoc_samples;
                    datafft_1pr = fft(data_ch(lb:ub));
                    sam_complex_50(k) = datafft_1pr(amf+inc);
                end

                sam_50epoch = sam_complex_50;
                sam_50epoch = sam_50epoch(:);
                sam_realandimg(:,1) = real(sam_50epoch);
                sam_realandimg(:,2) = imag(sam_50epoch);
                [nr,nc] = size(sam_realandimg);
                mu0 = [0;0];
                sam_realandimg_mean=mean(sam_realandimg);
                dev = sam_realandimg-kron(ones(nr,1),sam_realandimg_mean);
                s=dev'*dev/(nr-1);
                sinv=inv(s);
                wrk = sam_realandimg_mean'-mu0;
                tst = nr*wrk'*sinv*wrk;
                R_hot = tst*(nr-nc)/(nc*(nr-1));
                % tst1 = tst*(nr-nc)/(nc*(nr-1));
                % p_value= 1- fcdf(tst1,nc,nr-nc);
                % tst_matrix(sam_i,trig_i-159,chan_i) = tst;
                % p_value_matrix(sam_i,trig_i-159,chan_i) = p_value;
                R_hotelling_chlist(chan,:) = R_hot;
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end%channel
        %To extract all possible values of all scenarios(#trig, #sam) to compute joint alpha
        trigi = trigch-159;%@
        if sig_test ~= 2, test_trig_sam_cube(1,trigi,sami,:)= R_f_test_chlist; end      %@%F-test
        if sig_test ~= 1, test_trig_sam_cube(2,trigi,sami,:)= R_hotelling_chlist; end   %@%Hotelling
        sam
    end%for sam
    trigch
    waitbar((trigch-159)/num_stim,h);
end%for trigch

sig_cell = cell(num_stim,2);
sig_cell_tot = cell(numfreq,num_stim,2);

for trig_in = 1:num_stim,%5:8,%1:20,
    for stimf_in = 1:numfreq,%2,%1:5,
        [freq_bnd,bound] = find(freqlist(stimf_in)~=stimlist);

        %F-Test
        if sig_test ~= 2,
            flist_fp = test_trig_sam_cube(1,bound,stimf_in,:);
            flist_fp = sort(flist_fp(:));
            f_thresh = flist_fp(length(flist_fp)-num_false_pve+1);
            flist_tp = test_trig_sam_cube(1,trig_in,stimf_in,:);
            [flist_sig,ind_F] = find(flist_tp>f_thresh);
            roc_f(num_false_pve,trig_in,stimf_in) = length(flist_sig);
            sig_cell_tot(stimf_in,trig_in,1) = {[ind_F]}; 
        end
        %Hotelling
        if sig_test ~= 1,
            hlist_fp = test_trig_sam_cube(2,bound,stimf_in,:);
            hlist_fp = sort(hlist_fp(:));
            h_thresh = hlist_fp(length(hlist_fp)-num_false_pve+1);
            hlist_tp = test_trig_sam_cube(2,trig_in,stimf_in,:);
            [hlist_sig,ind_H] = find(hlist_tp>h_thresh);
            roc_h(num_false_pve,trig_in,stimf_in) = length(hlist_sig);
            sig_cell_tot(stimf_in,trig_in,2) = {[ind_H]}; 
        end
    end
    [aa,stim_ind] = find(freqlist == stimlist(trig_in));
    sig_cell(trig_in,1) = sig_cell_tot(stim_ind,trig_in,1);
    sig_cell(trig_in,2) = sig_cell_tot(stim_ind,trig_in,2);
end
sig_cell
sig_cell_tot
deltat = etime(clock,t0)
close(h);
%save test_trig_sam_cube_0380_1_24 test_trig_sam_cube% save F and H values
