function powspect = autopowspect(maxfreq, numfreqs, blocklen, autohisto, speedbins, balcounts)

%% maxfreq -- compute power spectrum from 0 - maxfreq
%% numfreqs -- number of frequency bins in the spectrum
%% blocklen -- length of each data block (in seconds)
%% ratehisto -- rate histogram data in matrix format (each row is the histogram for a movement episode of length blocklen)
%% speedbins -- an ordered list (length is same as number of rows in ratehisto) of which bin of the speed distribution each movement episide falls into
%% balcounts -- the counts for each speed bin in the balanced speed distribution (these define the "speed quotas" for the spectral analysis)

%% returns the power spectrum in 'powspect'

global iterat
f = (numfreqs/blocklen)*(0:(2^18))/2^19;  %frequency bins of the power spectrum
peakwidth = 1.5; %bandwidth (in Hz) on either side of the theta peak across which to intgrate for computing expected frequency value

nonzerobins=find(balcounts>0); %find the speed bin numbers with non-zero quotas

powspect=zeros(1,2^19); %initialize the power spectrum to zero for cumulative averaging
    
for i=1:100 %we will run the FFT a total of 100 times, randomly refilling the speed quotas each time
    
    ssignal=[]; %variable for accumulating the speed-balanced rate histogram on each pass through the data
    for j=nonzerobins %loop through the bins of the speed distribution which have nonzero quotas
        %%%first, try to fill the speed quota by sampling rows from the rate histogram WITHOUT replacement
        bsignal=[]; %variable for accumulating the rows of the rate histogram that fall into speed bin j
        thisbinepochs=autohisto(find(speedbins==j),:); %extract the rows of the autocorrelogram with running speeds that fall into speed bin j
        bsignal=[bsignal; thisbinepochs]; %accumulate the rows
        %%%second, if necessary, fill what remains of the speed quota by resampling rows from the rate histogram as little replacement as possible
        while (size(bsignal,1)<balcounts(j)) %while the quota for speed bin j has not yet been filled
            randthisbin=[thisbinepochs rand(size(thisbinepochs,1),1)]; %append a column of random numbers to the accumulated rate histogram
            randthisbin=sortrows(randthisbin,size(randthisbin,2)); %use the random number column to randomly order the rows
            randthisbin=randthisbin(:,1:size(thisbinepochs,2)); %remove he random number column now that its job is done
            if (size(bsignal,1)+size(randthisbin,1))<=balcounts(j) %if we have not yet exceeded our desired speed quota
                bsignal=[bsignal; randthisbin]; %then add all of the randomly ordered autocorrelogram rows to our accumulating total
            else %otherwise
                bsignal=[bsignal; randthisbin(1:(balcounts(j)-size(bsignal,1)),:)]; %add only as many rows as we need to fill the quota and be done
            end
        end
        ssignal=[ssignal; bsignal]; %accumulate the rows from speed bin j into the balanced rate histogram
    end
    
    Y=fft(sum(ssignal),2^19);
    Pyy = Y.* conj(Y) / 2^19;      
    powspect = powspect+Pyy/100;
    
 end