This was the answer I expect.
power = (abs(wave)).^2 ;        % compute wavelet power spectrum
incoi=(period(:)*(1./coi)<1);
Powd1 = zeros(numel(power(:, 1)), numel(power(1, :)));
for k = 1:numel(power(1, :))
    for j = 1:numel(power(:, 1))
    pow1 = power(j, k);
    incoi1 = incoi(j, k);
        if incoi1 == 1
            Powd1(j,k) = pow1;
        else            
            Powd1(j,k) = nan;
        end   
    end
end


