function [Sig, Perm] = CBPtest_FW(data, dataComp, c_alpha, p_alpha, n_random)

Data(:,1,:,:) = data;

if ~exist('c_alpha')
	cluster_a = 0.05;
else
	cluster_a = c_alpha;
end

if ~exist('p_alpha')
	pvlaue_a = 0.025;
else
	pvlaue_a = p_alpha;
end

if ~exist('n_random')
	bootNum = 5000;
else
	bootNum = n_random;
end

%%

nTime1 = size(data, 2);
nTime2 = size(data, 3);
subj  = size(data, 1); % Number of NoSect

if isempty(dataComp)
	baseData = zeros(subj, 1, nTime1, nTime2);
else
	baseData = repmat(dataComp, 1, nTime1, nTime2);
end
%% Parameters
allsubj_S.label     = cell(1,1);
allsubj_S.label{1}  = 'MLC11';
allsubj_S.powspctrm = Data;
allsubj_S.time      = 1:nTime2;
allsubj_S.freq      = 1:nTime1;
allsubj_S.dimord    = 'subj_chan_freq_time';

allsubj_Z.label     = cell(1,1);
allsubj_Z.label{1}  = 'MLC11';
allsubj_Z.powspctrm = baseData;
allsubj_Z.time      = 1:nTime2;
allsubj_Z.freq      = 1:nTime1;
allsubj_Z.dimord    = 'subj_chan_freq_time';

cfg                  = [];
cfg.channel          = 1;
cfg.latency          = 'all';
cfg.frequency        = 'all';
cfg.method           = 'montecarlo';
cfg.statistic        = 'depsamplesT';
cfg.correctm         = 'cluster';
cfg.clusteralpha     = cluster_a;
cfg.clusterstatistic = 'maxsum';
cfg.clusterthreshold = 'nonparametric_individual';
cfg.neighbours       = [];
cfg.clustertail      = 0;
cfg.alpha            = pvlaue_a;
cfg.correcttail      = 'alpha';
cfg.tail             = 0;
cfg.numrandomization = bootNum;

% Design Datarix
design = zeros(2,2*subj);
for i = 1:subj
  design(1,i) = i;
end
for i = 1:subj
  design(1,subj+i) = i;
end
design(2,1:subj)        = 1;
design(2,subj+1:2*subj) = 2;

cfg.design = design;
cfg.uvar  = 1;
cfg.ivar  = 2;

[stat] = ft_freqstatistics(cfg, allsubj_S, allsubj_Z);

Sig = stat;

%% Visualization
if size(stat.posclusters,2) == 0 && size(stat.negclusters,2) ==0
    return
end
% 

if size(stat.posclusters,2) ~= 0 && size(stat.negclusters,2) ==0 
    
    Prob_pos = [];
    
    for i = 1:size(stat.posclusters,2)
        Prob_pos = [Prob_pos; stat.posclusters(i).prob];
    end
    
    Idx = find(Prob_pos < 0.01);
    Idx = max(Idx);
    mat = squeeze(stat.posclusterslabelmat);
    
    if isempty(Idx)
        mat = repmat(0,nTime1,nTime2);
    else
        
        for j = 1:nTime2
            nIdx = find(mat(:,j) > Idx);
            mat(nIdx,j) = 0;
            
            pIdx = find(mat(:,j) <= Idx & mat(:,j) > 0);
            mat(pIdx,j) = 1;
        end
    end
    sigWind = mat;

end
% 
if size(stat.posclusters,2) == 0 && size(stat.negclusters,2) ~=0
    Prob_neg = [];
    
    for i = 1:size(stat.negclusters,2)
        Prob_neg = [Prob_neg; stat.negclusters(i).prob];
    end
    
    Idx = find(Prob_neg < 0.01);
    Idx = max(Idx);
    mat = squeeze(stat.negclusterslabelmat);
    
    if isempty(Idx)
        mat = repmat(0,nTime1,nTime2);
    else
        for j = 1:nTime2
            nIdx = find(mat(:,j) > Idx);
            mat(nIdx,j) = 0;
            
            pIdx = find(mat(:,j) <= Idx & mat(:,j) > 0);
            mat(pIdx,j) = 1;
        end
    end
    sigWind = mat;
end

if size(stat.posclusters,2) ~= 0 && size(stat.negclusters,2) ~=0
    Prob_pos = []; Prob_neg = [];
    
    for i = 1:size(stat.posclusters,2)
        Prob_pos = [Prob_pos; stat.posclusters(i).prob];
    end
    
    Idx = find(Prob_pos < 0.01);
    Idx = max(Idx);
    mat1 = squeeze(stat.posclusterslabelmat);
    
    if isempty(Idx)
        mat1 = repmat(0,nTime1,nTime2);
    else
        for j = 1:nTime2
            nIdx = find(mat1(:,j) > Idx);
            mat1(nIdx,j) = 0;
            
            pIdx = find(mat1(:,j) <= Idx & mat1(:,j) > 0);
            mat1(pIdx,j) = 1;
        end
    end
    
    for i = 1:size(stat.negclusters,2)
        Prob_neg = [Prob_neg; stat.negclusters(i).prob];
    end
    
    Idx = find(Prob_neg < 0.01);
    Idx = max(Idx);
    mat = squeeze(stat.negclusterslabelmat);
    
    if isempty(Idx)
        mat = repmat(0,nTime1,nTime2);
    else
        
        for j = 1:nTime2
            nIdx = find(mat(:,j) > Idx);
            mat(nIdx,j) = 0;
            
            pIdx = find(mat(:,j) <= Idx & mat(:,j) > 0);
            mat(pIdx,j) = 1;
        end
    end
    
    for i = 1:nTime1
        for j = 1:nTime2
            if mat(i,j) ~= mat1(i,j)
                mat(i,j) = 1;
            end
        end
    end
    sigWind = mat;
end
% 
Perm = sigWind;
