def detect_face_12net(cls_prob,roi,out_side,scale,width,height,threshold):
    in_side = 2*out_side+11
    stride = 0
    if out_side != 1:
        stride = float(in_side-12)/(out_side-1)
    (x,y) = np.where(cls_prob>=threshold)
    boundingbox = np.array([x,y]).T
    bb1 = np.fix((stride * (boundingbox) + 0 ) * scale)
    bb2 = np.fix((stride * (boundingbox) + 11) * scale)
    boundingbox = np.concatenate((bb1,bb2),axis = 1)
    dx1 = roi[0][x,y]
    dx2 = roi[1][x,y]
    dx3 = roi[2][x,y]
    dx4 = roi[3][x,y]
    score = np.array([cls_prob[x,y]]).T
    offset = np.array([dx1,dx2,dx3,dx4]).T
    boundingbox = boundingbox + offset*12.0*scale
    rectangles = np.concatenate((boundingbox,score),axis=1)
    rectangles = rect2square(rectangles)
    pick = []
    for i in range(len(rectangles)):
	x1 = int(max(0     ,rectangles[i][0]))
	y1 = int(max(0     ,rectangles[i][1]))
	x2 = int(min(width ,rectangles[i][2]))
	y2 = int(min(height,rectangles[i][3]))
	sc = rectangles[i][4]
	if x2>x1 and y2>y1:
	    pick.append([x1,y1,x2,y2,sc])
    return NMS(pick,0.5,'iou')

def symmetrized_dot_cloud(waveform, angle=60, top=50, lag=1):
    """
    same as symmetrized dots but returns a 2d point cloud instead of a boolean image. note the
    return value is a tuple (x, y) of int ndarrays, which is equivalent to (cols, rows).
    """

    # normalize the waveform to 0-top 
    waveform = normalize(waveform, top)

    # roll the waveform by lag to get the angles we'll use
    roll_waveform = np.roll(waveform, lag)

    # now compute the angles
    first_angle = np.concatenate(
        [base_angle + roll_waveform for base_angle in np.arange(0, 360, angle)])
    second_angle = np.concatenate(
        [base_angle - roll_waveform for base_angle in np.arange(0, 360, angle)])
    all_angles = np.deg2rad(np.concatenate([first_angle, second_angle]))

    # tile the original waveform until it matches the length of the angles
    num_repeats = int(all_angles.size / waveform.size)
    waveform = np.tile(waveform, num_repeats)

    # now make the point cloud
    x_cols = waveform * np.cos(all_angles) + top
    y_rows = waveform * np.sin(all_angles) + top

    return (np.fix(x_cols).astype(np.int), np.fix(y_rows).astype(np.int))

def dec2Dec2(Dec):
	'''
	change degree to dec2
	dec2 is the degree: np.array([ XXdegree, XXminute, XXsecond])
	if dec2 is negative, it should be like np.array([-XXdegree, -XXminute, -XXsecond])

	examples:
		for single degree
			input: np.array([degree])

			dec2Dec2(np.array([45.51]))
			output = 
			array([[ 45.,  30.,  36.]])
			it is the dec2 representing 45 degree 30 minute 36 second

		for array degrees
			input: np.array([degree1, degree2, degree3....])

			dec2Dec2(np.array([-256.32, -135.5, 10.3234, 256.3333333333]))
			output = 
			array([[-256.        ,  -19.        ,  -12.        ],
			       [-135.        ,  -30.        ,    0.        ],
			       [  10.        ,   19.        ,   24.24      ],
			       [ 256.        ,   19.        ,   59.99999988]])
			they are 4 dec2s
	'''
	Dec = np.float64(Dec)
	degree = np.fix(Dec)
	#print 'degree=',degree
	_point = (Dec-degree)*60.0
	minute = np.fix(_point)
	#print 'minute=',minute
	arc = (_point-minute)*60.0
	#print 'arc=',arc
	return np.array([degree,minute,arc]).T   

def generate_bounding_box(imap, reg, scale, threshold):
    """Use heatmap to generate bounding boxes"""
    # pylint: disable=too-many-locals
    stride = 2
    cellsize = 12

    imap = np.transpose(imap)
    d_x1 = np.transpose(reg[:, :, 0])
    d_y1 = np.transpose(reg[:, :, 1])
    d_x2 = np.transpose(reg[:, :, 2])
    d_y2 = np.transpose(reg[:, :, 3])
    dim_y, dim_x = np.where(imap >= threshold)
    if dim_y.shape[0] == 1:
        d_x1 = np.flipud(d_x1)
        d_y1 = np.flipud(d_y1)
        d_x2 = np.flipud(d_x2)
        d_y2 = np.flipud(d_y2)
    score = imap[(dim_y, dim_x)]
    reg = np.transpose(np.vstack([d_x1[(dim_y, dim_x)], d_y1[(dim_y, dim_x)],
                                  d_x2[(dim_y, dim_x)], d_y2[(dim_y, dim_x)]]))
    if reg.size == 0:
        reg = np.empty((0, 3))
    bbox = np.transpose(np.vstack([dim_y, dim_x]))
    q_1 = np.fix((stride * bbox + 1) / scale)
    q_2 = np.fix((stride * bbox + cellsize - 1 + 1) / scale)
    boundingbox = np.hstack([q_1, q_2, np.expand_dims(score, 1), reg])
    return boundingbox, reg

def epoch_plot(t, y, period, name='none'):

    phase = (t/period)*1.0 - np.fix(t/period)
    sortidxs = np.argsort(phase)
    sorted_phase = phase[sortidxs]
    sorted_y = y[sortidxs]
    sorted_dates = t[sortidxs]
    pp = np.concatenate([sorted_phase, sorted_phase+1.0], 1)
    yy = np.concatenate([sorted_y, sorted_y], 1)
    epochs = np.fix(sorted_dates/period)
    epoch_list = set(epochs)
    mi = min(y)-0.15
    ma2 = heapq.nlargest(2, y)[1]
    plt.clf()
    x = 0
    for i in epoch_list:
        period_date_idxs = np.asarray(np.where(epochs == i))
        period_mags = sorted_y[period_date_idxs]
        period_dates = sorted_dates[period_date_idxs]
        period_phases = sorted_phase[period_date_idxs]
        pp = np.concatenate([period_phases, period_phases + 1.0], 1)
        yy = np.concatenate([period_mags, period_mags], 1)
        plt.scatter(pp, yy, color=cm.jet(1.*x/len(epoch_list)), alpha=0.7, edgecolors='none')
        x +=1
    print ma2, mi
    plt.ylim(ma2, mi)
    #plt.gca().invert_yaxis()
    plt.text(-0.4, mi+0.1, "period = {} days".format(period), fontsize=14)
    plt.xlabel("phase")
    plt.ylabel("V-band magnitude")
    if name !='none':
        #plt.title("{}".format(name))
        plt.savefig("../output/{}_epoch.png".format(name))
    plt.show()

def LST(y, m, d, ut1, EL): 
	'''
	calculate the Local sidereal time
	input: year,month,day(do not have to be integer),ut1,Longitude of the site (all input can be array)
	output: Local sidereal time in degree
	examples:
		LST(2012,10,29,19+53./60 -8 ,116.4) # the Local sidereal time in Beijing(+8) at 2012-10-29 19:53:00
		output = 
		332.86374268340001 # (degree)
	'''
	def J0(year, month, day):
		'''temp function of LST'''
		j0 = 367.0*year - np.fix(7.0*(year + np.fix((month + 9)/12.0))/4.0)+ np.fix(275.0*month/9) + day + 1721013.5
		return j0

	y = np.float64(y)
	m = np.float64(m)
	d = np.float64(d)
	ut1 = np.float64(ut1)
	EL = np.float64(EL)

	dd = np.fix(d)
	time = (d-dd)*24.0
	ut = (ut1 + time)%24
	j0 = J0(y, m, dd)
	j = (j0 - 2451545.0)/36525.0
	g0 = 100.4606184 + (36000.77004*j) + (0.000387933*j**2)- 2.583e-8*j**3
	g0 = g0%360
	gst = g0 + (360.98564724*ut/24.0)
	lst = gst + EL
	lst = lst - 360.0*np.fix(lst/360.0)
	return lst

def idftreal(A, N, M):
    '''
    Calculates multiple 1D inverse DFT from complex to real

    Input:
        A - input complex vectors in column form (samples from zero to Nyquist)
        N - number of output samples (= number of implied complex input samples)
        M - number of input vectors

    Based on idftreal.m by R.G. Pratt
    '''

    a = np.zeros((N, M))
    n = np.arange(N).reshape((N,1))

    # Set maximum non-Nyquist frequency index (works for even or odd N)
    imax = np.int(np.fix((N+1)//2)-1)
    k1 = np.arange(np.fix(N//2)+1)               # Freq indices from zero to Nyquist
    k2 = np.arange(1, imax+1)                   # Freq indices except zero and Nyquist
    nk1 = n * k1.T
    nk2 = n * k2.T
    w = np.exp(-2j*np.pi / N)
    W = w**nk1
    W2 = w**nk2
    W[:,1:imax+1] += W2                         # Add two matrices properly shifted
    a = np.dot(W, A[:np.fix(N//2)+1,:M]).real    # (leads to doubling for non-Nyquist)

    return a

def cos_pattern(AO_resolution, freq):
# Experiment spec: sinewave on angles of stepper motor:
# 60 degree amplitude (max), 15 Hz (max)
# from Karin and Holger 12-apr-2011 15:00:00

    step = 360.0/5000      #...step(/microstep) angle (degree)
    amp = 60       /2  #...sine wave amplitude (peak to peak degree)
#    freq = 1      #...sine wave frequency (Hz)

    Y = numpy.arange(0,amp,step)   
    X = (1/(2*numpy.pi*freq))*numpy.arcsin(Y/amp)

#    AO_resolution = 2.5e-6          # (us) resolution: 400k Hz maximum, the period of which is 2.5 us
#    AO_resolution = 5e-6          # (us) resolution: 200k Hz maximum, the period of which is 5 us
#    AO_resolution = 1e-4          # (us) resolution: 10k Hz maximum, the period of which is 100 us

    Xs = numpy.fix(X / AO_resolution )

    SM_sig = numpy.zeros(numpy.fix(Xs[numpy.size(Xs)-1])+1)
    for i in range(numpy.size(Xs)): SM_sig[Xs[i]]=1

    rev_SM_sig = numpy.zeros(numpy.fix(Xs[numpy.size(Xs)-1])+1)
    for i in range(numpy.size(SM_sig)): rev_SM_sig[i]=SM_sig[numpy.size(SM_sig)-1-i]

#    duration = 2 * X[numpy.size(X)-1]    # (second) duration of the signal train

    wavelet = 5 * numpy.append( rev_SM_sig , SM_sig)
    wavelet = numpy.append(wavelet, wavelet)
    wavelet = numpy.append(wavelet, numpy.array([0]))

    return wavelet

def nlfer(signal, pitch, parameters):

    #---------------------------------------------------------------
    # Set parameters.
    #---------------------------------------------------------------
    N_f0_min = np.around((parameters['f0_min']*2/float(signal.new_fs))*pitch.nfft)
    N_f0_max = np.around((parameters['f0_max']/float(signal.new_fs))*pitch.nfft)

    window = hanning(pitch.frame_size+2)[1:-1]
    data = np.zeros((signal.size))  #Needs other array, otherwise stride and
    data[:] = signal.filtered     #windowing will modify signal.filtered

    #---------------------------------------------------------------
    # Main routine.
    #---------------------------------------------------------------
    samples = np.arange(int(np.fix(float(pitch.frame_size)/2)),
                        signal.size-int(np.fix(float(pitch.frame_size)/2)),
                        pitch.frame_jump)

    data_matrix = np.empty((len(samples), pitch.frame_size))
    data_matrix[:, :] = stride_matrix(data, len(samples),
                                    pitch.frame_size, pitch.frame_jump)
    data_matrix *= window

    specData = np.fft.rfft(data_matrix, pitch.nfft)

    frame_energy = np.abs(specData[:, N_f0_min-1:N_f0_max]).sum(axis=1)
    pitch.set_energy(frame_energy, parameters['nlfer_thresh1'])
    pitch.set_frames_pos(samples)

def get_cofe_target(ut,lat,lon,target):
#function to use ephem to find az and el of specified target for COFE
#parameters: UT, LAt Lon Target
    cofe=ephem.Observer()
    cofe.elevation=0.0
    year=2013
    month=10
    day=04
    az=[]
    el=[]
    for u,la,lo in zip(ut,lat,lon):
        if (u >24):
            u=u-24
            day=05
        hour=int(np.fix(u))
        minute=(u-hour)*60
        iminute=int(np.fix(minute))
        second=int(np.fix((minute-iminute)*60))
        datestring=str(year)+'/'+str(month)+'/'+str(day)+' '+str(hour)+':'+str(iminute)+':'+str(second)
        datestring
        cofe.date=datestring
        cofe.lon=str(rtd*lo)
        cofe.lat=str(rtd*la)
        pos=ephem.__getattribute__(target)(cofe)
        az.append(pos.az)
        el.append(pos.alt)
        
    return np.array(az),np.array(el)

 def plot_imfs(self, time, tstart=None, tend=None, tunit='s', savefig_name=''):
     time = time - time[0]
     dt = time[1] - time[0]
     if tstart != None:
         tstart = int(np.fix(tstart / dt))
     if tend != None:
         tend = int(np.fix(tend / dt))
     num_subplot = self.num_imf + 1
     fig = plt.figure()
     ax = fig.add_subplot(num_subplot, 1, 1)
     ax.plot(time[tstart:tend], self.input_signal[tstart:tend])
     ax.set_ylabel('Data')
     ax.get_yaxis().set_label_coords(-0.1,0.5)
     for i in range(self.num_imf - 1):
         ax = fig.add_subplot(num_subplot, 1, i + 2)
         ax.plot(time[tstart:tend], self.imfs[:,i][tstart:tend])
         ax.set_ylabel(r'$C_{' + str(i + 1) + '}$')
         ax.get_yaxis().set_label_coords(-0.1,0.5)
     ax = fig.add_subplot(num_subplot, 1, num_subplot)
     ax.plot(time[tstart:tend], self.imfs[:,-1][tstart:tend])
     ax.get_yaxis().set_label_coords(-0.1,0.5)
     ax.set_ylabel('Trend')
     ax.set_xlabel('Time (' + tunit + ')')
     if savefig_name == '':
         plt.show()
     else:
         plt.savefig(savefig_name, format='eps', dpi=1000)

def generateBoundingBox(imap, reg, scale, t):
    # use heatmap to generate bounding boxes
    stride = 2
    cellsize = 12

    imap = np.transpose(imap)
    dx1 = np.transpose(reg[:, :, 0])
    dy1 = np.transpose(reg[:, :, 1])
    dx2 = np.transpose(reg[:, :, 2])
    dy2 = np.transpose(reg[:, :, 3])
    y, x = np.where(imap >= t)
    if y.shape[0] == 1:
        dx1 = np.flipud(dx1)
        dy1 = np.flipud(dy1)
        dx2 = np.flipud(dx2)
        dy2 = np.flipud(dy2)
    score = imap[(y, x)]
    reg = np.transpose(np.vstack([dx1[(y, x)], dy1[(y, x)], dx2[(y, x)], dy2[(y, x)]]))
    if reg.size == 0:
        reg = np.empty((0, 3))
    bb = np.transpose(np.vstack([y, x]))
    q1 = np.fix((stride * bb + 1) / scale)
    q2 = np.fix((stride * bb + cellsize - 1 + 1) / scale)
    boundingbox = np.hstack([q1, q2, np.expand_dims(score, 1), reg])
    return boundingbox, reg

def intermediate_profile(x, xhinge, delta):
    """Generate an intermediate profile of some quantity.
    Ferron et. al. 1998

    Returns the 'top down' and 'bottom up' profiles as well as the average
    of the two.

    """

    xf = np.flipud(x)

    xtd = np.zeros_like(x)
    xbu = np.zeros_like(x)

    ntd = np.fix(x[0]/delta - xhinge/delta)
    nbu = np.fix(xf[0]/delta - xhinge/delta)

    xtd[0] = xhinge + ntd*delta
    xbu[0] = xhinge + nbu*delta

    for i in xrange(len(x) - 1):
        ntd = np.fix(x[i+1]/delta - xtd[i]/delta)
        nbu = np.fix(xf[i+1]/delta - xbu[i]/delta)

        xtd[i+1] = xtd[i] + ntd*delta
        xbu[i+1] = xbu[i] + nbu*delta

    xbu = np.flipud(xbu)

    xav = (xtd + xbu)/2.

    return xtd, xbu, xav

    def _simulate_image(self):
        """
        Generates the fake output.
        """
        with self._acquisition_init_lock:
            pos = self.align.position.value
            logging.debug("Simulating image shift by %s", pos)
            ac, bc = pos.get("a"), pos.get("b")
            ang = math.radians(135)
            # AB->XY
            xc = -(ac * math.sin(ang) + bc * math.cos(ang))
            yc = -(ac * math.cos(ang) - bc * math.sin(ang))
            pixelSize = self.fake_img.metadata[model.MD_PIXEL_SIZE]
            self.fake_img.metadata[model.MD_ACQ_DATE] = time.time()
            x_pxs = xc / pixelSize[0]
            y_pxs = yc / pixelSize[1]

            # Image shifted based on LensAligner position
            z = 1j  # imaginary unit
            self.deltar = x_pxs
            self.deltac = y_pxs
            nr, nc = self.fake_img.shape
            array_nr = numpy.arange(-numpy.fix(nr / 2), numpy.ceil(nr / 2))
            array_nc = numpy.arange(-numpy.fix(nc / 2), numpy.ceil(nc / 2))
            Nr = fft.ifftshift(array_nr)
            Nc = fft.ifftshift(array_nc)
            [Nc, Nr] = numpy.meshgrid(Nc, Nr)
            sim_img = fft.ifft2(fft.fft2(self.fake_img) * numpy.power(math.e,
                            z * 2 * math.pi * (self.deltar * Nr / nr + self.deltac * Nc / nc)))
            output = model.DataArray(abs(sim_img), self.fake_img.metadata)
            return output

def interval2semitone(interval, accidentals):

  # Need it to be int
  import numpy as np
  interval = int(interval)
  
  # semitone equivalents for intervals
  #           interval   1,2,3,4,5,6,7         
  semitonetranslation = [0,2,4,5,7,9,11]

  semitone = 0

  success = True
  
  if (interval > 0) and (interval < 50):

    # semitone value is the number of semitones equivalent to the interval
    # added to the number of accidentals (sharps are positive, flats are
    # negative) and the number of octaves above the reference note to
    # account for extensions 

    des_index = int(np.mod(interval,8) + np.fix(interval/8))-1
    semitone = int(semitonetranslation[des_index] + accidentals + 12*np.fix(interval/8))

  else:
    success = False    
    print 'Error in interval2semitone: out of range interval'

  return semitone, success

def jd2gps(jd):
    """
    % JD2GPS  Converts Julian date to GPS week number (since
    %   1980.01.06) and seconds of week. 
    % Usage:   [gpsweek,sow,rollover]=jd2gps(jd)
    % Input:   jd       - Julian date
    % Output:  gpsweek  - GPS week number
    %          sow      - seconds of week since 0 hr, Sun.
    %          rollover - number of GPS week rollovers (modulus 1024)

    % Copyright (c) 2011, Michael R. Craymer
    % All rights reserved.
    % Email: [email protected]
    """

    jdgps = cal2jd(1980,1,6);    # beginning of GPS week numbering
    nweek = int(np.fix((jd-jdgps)/7.))
    sow = (jd - (jdgps+nweek*7)) * 3600*24
    rollover = np.fix(nweek/1024)  # rollover every 1024 weeks
    gpsweek = int(nweek)

# rollover is being returned as an array?
# should just be an int

    return gpsweek,sow,rollover 

def coherr(C,J1,J2,p=0.05,Nsp1=None,Nsp2=None):
    """
    Function to compute lower and upper confidence intervals on
    coherency (absolute value of coherence).

    C:            coherence (real or complex)
    J1,J2:        tapered fourier transforms
    p:            the target P value (default 0.05)
    Nsp1:         number of spikes in J1, used for finite size correction.
    Nsp2:         number of spikes in J2, used for finite size correction.
                  Default is None, for no correction

    Outputs:
    CI:           confidence interval for C, N x 2 array, (lower, upper)
    phi_std:      stanard deviation of phi, N array
    """
    from numpy import iscomplexobj, absolute, fix, zeros, setdiff1d, real, sqrt,\
         arctanh, tanh
    from scipy.stats import t

    J1 = _combine_trials(J1)
    J2 = _combine_trials(J2)
    N,K = J1.shape
    assert J1.shape==J2.shape, "J1 and J2 must have the same dimensions."
    assert N == C.size, "S and J lengths don't match"
    if iscomplexobj(C): C = absolute(C)

    pp = 1 - p/2
    dof = 2*K
    dof1 = dof if Nsp1 is None else fix(2.*Nsp1*dof/(2.*Nsp1+dof))
    dof2 = dof if Nsp2 is None else fix(2.*Nsp2*dof/(2.*Nsp2+dof))
    dof = min(dof1,dof2)

    Cerr = zeros((N,2))
    tcrit = t(dof-1).ppf(pp).tolist()
    atanhCxyk = zeros((N,K))
    phasefactorxyk = zeros((N,K),dtype='complex128')

    for k in xrange(K):
        indxk = setdiff1d(range(K),[k])
        J1k = J1[:,indxk]
        J2k = J2[:,indxk]
        eJ1k = real(J1k * J1k.conj()).sum(1)
        eJ2k = real(J2k * J2k.conj()).sum(1)
        eJ12k = (J1k.conj() * J2k).sum(1)
        Cxyk = eJ12k/sqrt(eJ1k*eJ2k)
        absCxyk = absolute(Cxyk)
        atanhCxyk[:,k] = sqrt(2*K-2)*arctanh(absCxyk)
        phasefactorxyk[:,k] = Cxyk / absCxyk

    atanhC = sqrt(2*K-2)*arctanh(C);
    sigma12 = sqrt(K-1)* atanhCxyk.std(1)

    Cu = atanhC + tcrit * sigma12
    Cl = atanhC - tcrit * sigma12
    Cerr[:,0] = tanh(Cl / sqrt(2*K-2))
    Cerr[:,1] = tanh(Cu / sqrt(2*K-2))
    phistd = (2*K-2) * (1 - absolute(phasefactorxyk.mean(1)))
    return Cerr, phistd

def stim(w=360):
    global SM_step
    SM_step = 360.0/5000.0      #...step(/microstep) angle (degree)
    amp = 360.0 /2              #...sine wave amplitude (peak to peak degree)
    freq = 1                  #...sine wave frequency (Hz)

    Y = numpy.arange(0,amp,SM_step)   
    X = (1/(2*numpy.pi*freq))*numpy.arcsin(Y/amp)

    #w = 180.0                   #... deg/sec
    #tLmt = (360.0/w)/4/(SM_step)
    #tLmt = 1.0/5000
    global tLmt
    tLmt = SM_step / w          #... sec/step

    dX = numpy.zeros(numpy.size(X))
    dX[0] = X[0]
    for i in range(1, numpy.size(X)-1): dX[i] = X[i] - X[i-1]

    for i in range(1,numpy.size(X)):
        if dX[i] < tLmt: dX[i] = tLmt

    XX = numpy.zeros(numpy.size(X))
    for i in range(1,numpy.size(XX)): XX[i] = XX[i-1] + dX[i]

    #    pylab.plot(X,Y);pylab.show();

    AO_resolution = 5e-6          # (us) resolution: 200k Hz maximum, the period of which is 5 us    

    Xs = numpy.fix(XX / AO_resolution )

    SM_sig = numpy.zeros(numpy.fix(Xs[numpy.size(Xs)-1])+1)    #... numpy.size(Xs)-1
    for i in range(numpy.size(Xs)): SM_sig[Xs[i]]=1

    rev_SM_sig = numpy.zeros(numpy.fix(Xs[numpy.size(Xs)-1])+1)
    for i in range(numpy.size(SM_sig)): rev_SM_sig[i]=SM_sig[numpy.size(SM_sig)-1-i]

    #    duration = 2 * X[numpy.size(X)-1]    # (second) duration of the signal train

    wavelet = 5 * numpy.append( rev_SM_sig,SM_sig )
    wavelet = numpy.append(wavelet, wavelet)
    wavelet = numpy.append(wavelet, numpy.array([0]))

    SM_dir = numpy.append(5*numpy.ones(numpy.size(rev_SM_sig)+numpy.size(SM_sig)),numpy.zeros(numpy.size(rev_SM_sig)+numpy.size(SM_sig)))
    SM_dir = numpy.append(SM_dir, numpy.array([0]))

    # for easier data processing
    sig_head = 5*numpy.ones(20)
    sig_tail = numpy.append(numpy.zeros(20),5*numpy.ones(20))
    sig_tail = numpy.append(sig_tail,numpy.zeros(1))
    SM_pulse = numpy.append(sig_head,wavelet)
    SM_pulse = numpy.append(SM_pulse,sig_tail)
    SM_dir = numpy.append(sig_head,SM_dir)
    SM_dir = numpy.append(SM_dir,sig_tail)

    stim = numpy.append(SM_pulse, SM_dir)
                                                                    
    #wavelet = 5 * numpy.ones(40000)
    return stim

def jd2date(jd):
    """This function finds the Year, month, day, hour, minute and second
    given the Julian date.

    Algorithm     : Set up starting values
                    Find the elapsed days through the year in a loop
                    Call routine to find each individual value

    Author        : Capt Dave Vallado  USAFA/DFAS  719-472-4109  26 Feb 1990
    In Ada        : Dr Ron Lisowski    USAFA/DFAS  719-472-4110  17 May 1995
    In Matlab     : LtCol Thomas Yoder USAFA/DFAS  719-333-4110  Spring 2000
    In Python     : Shankar Kulumani   GWU         630-336-6257  2017 06 15

    Inputs        :
        JD          - Julian Date                          days from 4713 B.C.

    OutPuts       :
        Yr          - Year                                 1900 .. 2100
        Mon         - Month                                   1 .. 12
        D           - Day                                     1 .. 28,29,30,31
        H           - Hour                                    0 .. 23
        M           - Minute                                  0 .. 59
        S           - Second                                0.0 .. 59.999

    Locals        :
        days        - Day of year plus fraction of a day   days
        Tu          - Julian Centuries from 1 Jan 1900
        Temp        - Temporary Long_Float value
        LeapYrs     - Number of Leap years from 1900

    Constants     :
        None.

    Coupling      :
        DayofYr2MDHMS  - Finds Month, day, hour, minute and second given Days
        and Yr

    References    :
        1988 Almanac for Computers  pg. B2
        Escobal       pg. 17-19
        Kaplan        pg. 329-330
    """

    temp = jd - 2415019.5
    tu = temp / 365.25
    yr = 1900 + np.fix(tu)
    leapyrs = np.fix((yr - 1900 - 1) * 0.25)
    days = temp - ((yr - 1900) * 365.0 + leapyrs)

    # check for beginning of year
    if days < 1.0:
        yr = yr - 1
        leapyrs = np.fix((yr - 1900 - 1) * 0.25)
        days = temp - ((yr - 1900) * 365.0 + leapyrs)

    mon, day, h, m, s = dayofyr2mdhms(yr, days)

    return (yr, mon, day, h, m, s)

def sccs_bit_sync(y,Ns):
	"""
	rx_symb_d,clk,track = sccs_bit_sync(y,Ns)

	//////////////////////////////////////////////////////
	 Symbol synchronization algorithm using SCCS
	//////////////////////////////////////////////////////
	     y = baseband NRZ data waveform
	    Ns = nominal number of samples per symbol
	
	Reworked from ECE 5675 Project
	Translated from m-code version
	Mark Wickert April 2014
	"""
	# decimated symbol sequence for SEP
	rx_symb_d = np.zeros(int(np.fix(len(y)/Ns)))
	track = np.zeros(int(np.fix(len(y)/Ns)))
	bit_count = -1
	y_abs = np.zeros(len(y))
	clk = np.zeros(len(y))
	k = Ns+1 #initial 1-of-Ns symbol synch clock phase
	# Sample-by-sample processing required
	for i in xrange(len(y)):
	    #y_abs(i) = abs(round(real(y(i))))
	    if i >= Ns: # do not process first Ns samples
	        # Collect timing decision unit (TDU) samples
	        y_abs[i] = np.abs(np.sum(y[i-Ns+1:i+1]))
	        # Update sampling instant and take a sample
	        # For causality reason the early sample is 'i',
	        # the on-time or prompt sample is 'i-1', and  
	        # the late sample is 'i-2'.
	        if (k == 0):
	            # Load the samples into the 3x1 TDU register w_hat.
	            # w_hat[1] = late, w_hat[2] = on-time; w_hat[3] = early.
	            w_hat = y_abs[i-2:i+1]
	            bit_count += 1
	            if w_hat[1] != 0:
	                if w_hat[0] < w_hat[2]:
	                    k = Ns-1
	                    clk[i-2] = 1
	                    rx_symb_d[bit_count] = y[i-2-int(np.round(Ns/2))-1]
	                elif w_hat[0] > w_hat[2]:
	                    k = Ns+1
	                    clk[i] = 1
	                    rx_symb_d[bit_count] = y[i-int(np.round(Ns/2))-1]
	                else:
	                    k = Ns
	                    clk[i-1] = 1
	                    rx_symb_d[bit_count] = y[i-1-int(round(Ns/2))-1]
	            else:
	                k = Ns
	                clk[i-1] = 1
	                rx_symb_d[bit_count] = y[i-1-int(round(Ns/2))]
	            track[bit_count] = np.mod(i,Ns)
	    k -= 1
	# Trim the final output to bit_count
	rx_symb_d = rx_symb_d[:bit_count]
	return rx_symb_d, clk, track