epoch = eeg(:, epoch_start:epoch_end);

            % Calculate power spectra spectra:
            [spectra, freqs] = spectopo(epoch, 0, 128);

            % Save spectra:
            all_spectra(i,j,:,:) = spectra;
            
            % Get average power for each frequency band for each channel:

            % Theta:
            theta_start_sample = find(freqs==4);
            theta_end_sample = find(freqs==8);

            % Alpha:
            alpha_start_sample = find(freqs==8);
            alpha_end_sample = find(freqs==12);

            % Beta:
            beta_start_sample = find(freqs==12);
            beta_end_sample = find(freqs==30);

            % Delta:
            delta_start_sample = find(freqs==0.5);
            delta_end_sample = find(freqs==4);

            % Gamma:
            gamma_start_sample = find(freqs==30);
            gamma_end_sample = find(freqs==40);

            % Get power in each band for each channel:
            for k = 1:size(spectra, 1)

                disp(['Channel ', num2str(k)]);

                theta_spectrum = spectra(k, theta_start_sample:theta_end_sample);
                [avg_theta, std_theta] = normfit(theta_spectrum);

                alpha_spectrum = spectra(k, alpha_start_sample:alpha_end_sample);
                [avg_alpha, std_alpha] = normfit(alpha_spectrum);

                beta_spectrum = spectra(k, beta_start_sample:beta_end_sample);
                [avg_beta, std_beta] = normfit(beta_spectrum);

                delta_spectrum = spectra(k, delta_start_sample:delta_end_sample);
                [avg_delta, std_delta] = normfit(delta_spectrum);

                gamma_spectrum = spectra(k, gamma_start_sample:gamma_end_sample);
                [avg_gamma, std_gamma] = normfit(gamma_spectrum);
                
            end