Input: 
            Input - 3d array, with dimensions [nint, N_d0, nbl, nchan, npol]
            f_min,f_max - are the base-band begin and end frequencies.
                The frequencies can be entered in both MHz and GHz, units are factored out in all uses.
            maxDT - the maximal delay (in time bins) of the maximal dispersion.
                Appears in the paper as N_{\Delta}
                A typical input is maxDT = N_f
            dataType - To naively use FFT, one must use floating point types.
                Due to casting, use either complex64 or complex128.
            iteration num - Algorithm works in log2(Nf) iterations, each iteration changes all the sizes (like in FFT)
        Output: 
            5d array, with dimensions [nint, N_d1, nbl, nchan/2, npol]
        where N_d1 is the maximal number of bins the dispersion curve travels at one output frequency band
        
        For details, see algorithm 1 in Zackay & Ofek (2014)

def FDMT_iteration(datain, maxDT, nchan0, f_min, f_max, iteration_num, dataType):
    """
        Input: 
            Input - 3d array, with dimensions [nint, N_d0, nbl, nchan, npol]
            f_min,f_max - are the base-band begin and end frequencies.
                The frequencies can be entered in both MHz and GHz, units are factored out in all uses.
            maxDT - the maximal delay (in time bins) of the maximal dispersion.
                Appears in the paper as N_{\Delta}
                A typical input is maxDT = N_f
            dataType - To naively use FFT, one must use floating point types.
                Due to casting, use either complex64 or complex128.
            iteration num - Algorithm works in log2(Nf) iterations, each iteration changes all the sizes (like in FFT)
        Output: 
            5d array, with dimensions [nint, N_d1, nbl, nchan/2, npol]
        where N_d1 is the maximal number of bins the dispersion curve travels at one output frequency band
        
        For details, see algorithm 1 in Zackay & Ofek (2014)
    """

    nint, dT, nbl, nchan, npol = datain.shape
#    output_dims = list(input_dims)
    
    deltaF = 2**(iteration_num) * (f_max - f_min)/float(nchan0)
    dF = (f_max - f_min)/float(nchan0)
    # the maximum deltaT needed to calculate at the i'th iteration
    deltaT = int(np.ceil((maxDT-1) *(1./f_min**2 - 1./(f_min + deltaF)**2) / (1./f_min**2 - 1./f_max**2)))
    logger.debug("deltaT = {0}".format(deltaT))
    logger.debug("N_f = {0}".format(nchan0/2**(iteration_num)))
    
    dataout = np.zeros((nint, deltaT+1, nbl, nchan/2, npol), dataType)

    logger.debug('input_dims = {0}'.format(datain.shape))
    logger.debug('output_dims = {0}'.format(dataout.shape))
    
    # No negative D's are calculated => no shift is needed
    # If you want negative dispersions, this will have to change to 1+deltaT,1+deltaTOld
    # Might want to calculate negative dispersions when using coherent dedispersion, to reduce the number of trial dispersions by a factor of 2 (reducing the complexity of the coherent part of the hybrid)
    ShiftOutput = 0
    ShiftInput = 0

    F_jumps = nchan/2

    # For some situations, it is beneficial to play with this correction.
    # When applied to real data, one should carefully analyze and understand the effect of 
    # this correction on the pulse he is looking for (especially if convolving with a specific pulse profile)
    if iteration_num>0:
        correction = dF/2.
    else:
        correction = 0

    for i_F in range(F_jumps):
        f_start = (f_max - f_min)/float(F_jumps) * (i_F) + f_min
        f_end = (f_max - f_min)/float(F_jumps) *(i_F+1) + f_min
        f_middle = (f_end - f_start)/2. + f_start - correction

        # it turned out in the end, that putting the correction +dF to f_middle_larger (or -dF/2 to f_middle, and +dF/2 to f_middle larger)
        # is less sensitive than doing nothing when dedispersing a coherently dispersed pulse.
        # The confusing part is that the hitting efficiency is better with the corrections (!?!).

        f_middle_larger = (f_end - f_start)/2 + f_start + correction
        deltaTLocal = int(np.ceil((maxDT-1) *(1./f_start**2 - 1./(f_end)**2) / (1./f_min**2 - 1./f_max**2)))
        logger.debug('deltaT {0} deltaTLocal {1}'.format(deltaT, deltaTLocal))

        for i_dT in range(deltaTLocal+1):
            dT_middle = int(round(i_dT * (1./f_middle**2 - 1./f_start**2)/(1./f_end**2 - 1./f_start**2)))
            dT_middle_index = dT_middle + ShiftInput
            
            dT_middle_larger = int(round(i_dT * (1./f_middle_larger**2 - 1./f_start**2)/(1./f_end**2 - 1./f_start**2)))
            
            dT_rest = i_dT - dT_middle_larger
            dT_rest_index = dT_rest + ShiftInput
            
            logger.debug('{0}:{1}, {2}+{3}, {4} <= {5}, {6}'.format(i_T_min, i_T_max, i_dT, ShiftOutput, i_F, dT_middle_index, 2*i_F))

            # out of bounds data?
            i_T_min = 0
            i_T_max = dT_middle_larger
            dataout[i_T_min:i_T_max, i_dT + ShiftOutput, :, i_F, :] = datain[i_T_min:i_T_max, dT_middle_index, :, 2*i_F, :]

            # fully dedispersed data
            i_T_min = dT_middle_larger
            i_T_max = nint
            dataout[i_T_min:i_T_max, i_dT + ShiftOutput, :, i_F, :] = datain[i_T_min:i_T_max, dT_middle_index, :, 2*i_F, :] + datain[i_T_min - dT_middle_larger:i_T_max-dT_middle_larger, dT_rest_index, :, 2*i_F+1, :]
    
    return dataout