python类sinc()的实例源码

def rootRaisedCosine(t):
    bit_period = 1/BIT_FREQUENCY

    if (t== bit_period/(4*BETA)):
        return (BETA/(np.pi*np.sqrt(2*bit_period)) * \
                ((np.pi + 2)*np.sin(np.pi/(4*BETA)) + (np.pi - 2)*np.cos(np.pi/(4*BETA))))
    else:
        return (4 * BETA / np.pi / np.sqrt(bit_period) * \
            (np.cos((1 + BETA) * np.pi * t / bit_period) + \
             (1 - BETA) * np.pi / (4 * BETA) * np.sinc((1-BETA)*t/bit_period)) / \
                (1 - (4*BETA*t/bit_period)**2))

def get_CML(self, q, t):
        """
        Calculate C, M, L forming the elements of T matrix
        :param q: a shifted coordinate grid
        :param t: time
        :return: tuple C, M, L
        """
        assert q is self.x_plus or q is self.x_minus, \
            "the shifted coordinate (q) must be either x_plus or x_minus"

        # get the difference of adiabatic potential curves
        Vg_minus_Ve = (self._Vg_plus_Ve_x_plus if q is self.x_plus else self._Vg_minus_Ve_x_minus)

        Veg = self.Veg(q, t)

        D = Veg**2 + 0.25*Vg_minus_Ve**2
        np.sqrt(D, out=D)

        S = np.sinc(D * self.dt / np.pi)
        S *= self.dt

        C = D * self.dt
        np.cos(C, out=C)

        M = S * Vg_minus_Ve
        M *= 0.5

        L = S * Veg

        return C, M, L

def ef_cascade(screenpos, i, nletters):
    v = np.array([0, -1])
    d = lambda t : 1 if t < 0 else abs(np.sinc(t) / (1 + t**4))
    return lambda t: screenpos + v * 400 * d(t - 0.15 * i)

def collect_pixel(self, pixel_i, k,j,i):
        #print pixel_i, k,j,i
        t0 = time.time()
        #px_data = np.random.rand()
        #px_data = t0 - self.prev_px
        x0,y0 = self.pos
        x_set = self.stage.settings['x_position']
        y_set = self.stage.settings['y_position']
        x_hw = self.stage.settings.x_position.read_from_hardware(send_signal=False)
        y_hw = self.stage.settings.y_position.read_from_hardware(send_signal=False)
        if np.abs(x_hw - x0) > 1:
            self.log.debug('='*60)
            self.log.debug('pos      {} {}'.format(x0, y0))
            self.log.debug('settings {} {}'.format(x_set, y_set))
            self.log.debug('hw       {} {}'.format(x_hw, y_hw))            
            self.log.debug('settings value delta {} {}'.format(x_set-x0, y_set-y0))
            self.log.debug('read_hw  value delta {} {}'.format(x_hw-x0, y_hw-y0))
            self.log.debug('='*60)

        x = x_hw
        y = y_hw

        px_data = np.sinc((x-50)*0.05)**2 * np.sinc(0.05*(y-50))**2 #+ 0.05*np.random.random()
        #px_data = (x-xhw)**2 + ( y-yhw)**2
        #if px_data > 1:
        #    print('hw', x, xhw, y, yhw)
        self.display_image_map[k,j,i] = px_data
        if self.settings['save_h5']:
            self.test_data[k,j,i] = px_data 
        time.sleep(self.settings['pixel_time'])
        #self.prev_px = t0

def Sinc(freq=440, amp=1.0, offset=0):
    """Makes a Sinc function.

    freq: float frequency in Hz
    amp: float amplitude, 1.0 is nominal max
    offset: float phase offset in radians

    returns: Sinusoid object
    """
    return Sinusoid(freq, amp, offset, func=np.sinc)

def intensitiesFFFaster(nq, dq, dIntegrand, keys, ffDict, dr):
    """ uses atomistic form factors """
    dIntensity = np.zeros((nq, 2), dtype=float)
    partInt = np.zeros((nq, len(dIntegrand[0])))
    nameList = [k.split(",") for k in keys]
    # print "nameList=",nameList
    qList = np.zeros(nq)
    qList[:] = [float(i * dq) for i in range(nq)]
    dIntensity[:, 0] = qList[:]
    partInt[:, 0] = qList[:]
    rList = dIntegrand[:, 0]
    # for j in range(0,len(dIntegrand)):
    #    r=dIntegrand[j,0]
    #    sinc=j0(q*r)
    formFacProd = np.zeros((nq, len(dIntegrand[0])))
    for i in range(nq):
        sincList = np.sinc(rList * qList[i] / math.pi) * dr
        for k in range(1, len(dIntegrand[0])):
            # print k
            formFacProd[i, k] = ff.fiveGaussian(ffDict[nameList[k - 1][0]], qList[i])\
                * ff.fiveGaussian(ffDict[nameList[k - 1][1]], qList[i])
            partInt[i, k] += (sincList[:] * dIntegrand[:, k]
                              ).sum() * formFacProd[i, k]
    for i in range(nq):
        dIntensity[i, 1] = partInt[i, 1:].sum()
    return partInt, dIntensity

def dirichlet(x):
    return numpy.sinc(x)

def dirichlet(x):
    return numpy.sinc(x)

def initLanczos(self, filtOrder):
        self.filtOrder = filtOrder

        if self.filtOrder % 2 != 0:
             raise Exception('Invalid filtOrder: ' + str(self.filtOrder) +
                 ' Must be an even integer.')

        radius = self.filtOrder // 2
        win = np.sinc(np.linspace(-radius, radius, self.filtOrder+1) / float(radius)) # lanczos
        #win = spsig.hamming(self.filtOrder+1) # sinc-hamming

        # this should be automated somehow XXX - idfah
        if self.filtOrder <= 6:
            cutoff = 2*0.570

        elif self.filtOrder <= 8:
            cutoff = 2*0.676

        elif self.filtOrder <= 12:
            cutoff = 2*0.781

        elif self.filtOrder <= 16:
            cutoff = 2*0.836

        elif self.filtOrder <= 32:
            cutoff = 2*0.918

        elif self.filtOrder <= 64:
            cutoff = 2*0.959

        # need to fix for multiple pool sizes XXX - idfah
        cutoff /= float(self.poolSize)

        taps = cutoff * np.linspace(-radius, radius, self.filtOrder+1, dtype=self.dtype)
        impulseResponse = cutoff * np.sinc(taps) * win

        self.filters = []
        nReadoutLayers = 1 if self.nHidden is None else 2
        for ni, no in self.layerDims[:-nReadoutLayers]:
            noEmb = no*(self.filtOrder+1) # no outs after filter embedding

            filtMat = np.zeros(noEmb*2, dtype=self.dtype)
            filtMat[noEmb-1::-no] = impulseResponse

            # filters strided for embedding
            sz = filtMat.itemsize
            filtMat = npst.as_strided(filtMat, (no,noEmb), strides=(sz,sz))[::-1].T
            self.filters.append(filtMat.copy())

def rrc(t):
    '''
    Input: T, evaluation point (seconds)
    Output: value of root-raised-cosine at time T
    '''
    # Delay between two bits
    bit_period = 1/BIT_FREQUENCY
    # Total amount of bits to transmit
    nb_bits = len(LIST_OF_BITS)
    # To be returned (sum of contributions)
    s = 0.0
    # Max value of rrc
    m = 4*BETA/np.pi/np.sqrt(bit_period) + (1-BETA)/np.sqrt(bit_period) + sum(abs(2*rootRaisedCosine(i*bit_period)) for i in range(1, TRUNCATION))

    if(t < - TRUNCATION * bit_period  or t >= (nb_bits + TRUNCATION) * bit_period):
        # T out of support
        r = 0.0
    else:
        # Bits that will affect function at time T
        relevant_bits = np.zeros(2*TRUNCATION+1)
        for i in range(2*TRUNCATION+1):
            j = t/bit_period + i - TRUNCATION
            j = int(j) if int(j) <= j else int(j) - 1
            if(j >= 0 and j < nb_bits):
                relevant_bits[i] = -1 if LIST_OF_BITS[j] == '0' else 1

        for i in range(2*TRUNCATION+1):
            tt = t/bit_period
            tt = t - int(tt)*bit_period if int(tt) <= tt else t - (int(tt)-1)*bit_period
            if(t == bit_period * (1 / 4 / BETA + (i - TRUNCATION))):
                # L'Hospital's rule because of potential discontinuity
                s += relevant_bits[i] * BETA / np.pi / np.sqrt(2*bit_period) * 1 / m * \
                     ((np.pi + 2) * np.sin(np.pi/4/BETA) + \
                      (np.pi - 2) * np.cos(np.pi/4/BETA))
            else:
                # General case formula
                s += relevant_bits[i] * 4*BETA/np.pi/np.sqrt(bit_period) * 1 / m * \
                     (np.cos((1 + BETA) * np.pi * ((tt / bit_period - (i-TRUNCATION)))) + \
                      (1 - BETA) * np.pi / 4 / BETA * \
                      np.sinc((1 - BETA) * (tt / bit_period - (i-TRUNCATION))))/ \
                      (1 - (4*BETA*(tt / bit_period - (i-TRUNCATION)))**2)

    return s

###   ###   ###   ###   ###   ###   ###

def noise_power(self, freq):
        """Returns a function to calculate the noise PS at the given freq.
        """

        z = freq_to_z(freq)

        beam_size = self.beam_size(freq)
        A_pix = beam_size**2
        A_survey = self.f_sky * 4 * np.pi

        tau = A_pix / A_survey * units.year * self.num_year * self.beam_num

        # Calculate the comoving size of a frequency bin (at the given freq)
        d = self.proper_distance
        dxf = (d(freq_to_z(freq - self.freq_width)) - d(z))

        # Define the window function in k-space for the parallel and perpendicular directions.
        # Use a sinc function for parallel as it is the FT of a top-hat bin. This is probably a bad choice.
        def window_par(kpar):
            y = kpar * dxf / (4 * np.pi)
            return np.sinc(y) * (np.abs(y) < 1.0)

        # Azimuthally average over the X and Y window functions to produce an
        # overall k_perp window function. Do this by averaging for a set number
        # of points k values and then generating an interpolating function to
        # appeoximate the full result.
        def _int(phi, k):
            # Integrand to average over
            x = (3e2 * k * d(z)) / (freq * 2 * np.pi)

            xx = x * np.cos(phi)
            xy = x * np.sin(phi)
            return (self.window_x(xx) * self.window_y(xy))**2

        def _w_xy_average(k):
            # Full averaged window function
            return scipy.integrate.fixed_quad(_int, 0, 2 * np.pi, args=(k,), n=1024)[0]

        # Generate a log interpolated approximation
        k_val = np.linspace(0, self.kmax, 256)
        int_val = np.array([_w_xy_average(k)**0.5 for k in k_val])
        _w_perp_interp = scipy.interpolate.interp1d(k_val, np.log(int_val))

        def window_perp(kperp):
            return np.exp(_w_perp_interp(kperp))

        # Calculate the comoving volume of a single pixel (beam)
        V_pix = A_pix * d(z)**2 * dxf

        # Receiver temperature contribution to instrumental Stokes I
        T_recv_I = self.T_recv / 2**0.5

        return inv_noise_ps_21cm(T_recv_I + self.T_sky(freq), tau, V_pix,
                                 self.freq_width, window_par, window_perp)

def focus_multiprocessing(self, row):
        """
        Focus SAR image with TDBP algorithm. NOTE: Image must be range compressed.
        It uses local squint angle (antenna coordinate system) and distances to target to focus.

        Parameters
        ----------
        row: int.
          image row to be focus.

        Returns
        -------
        list of numpy complex.
          List containing entire focused row (numpy complex data) calculated in parallel mode.
        """

        # Light speed.
        c = 300000000.0
        # SAR bandwidth, central frequency and lambda.
        sar_B = self.param.get_float_parameter("Radar/B")
        sar_f0 = self.param.get_float_parameter("Radar/f0")
        sar_lambda = c/sar_f0

        nt_fast_time = self.simulated_image.traj.nt
        nt_slow_time = self.simulated_image.Nt

        # Partial row calculated in parallel mode focusing.
        partial_row = np.empty(self.ny, dtype=np.complex128)

        x_foc_ind = row
        for y_foc_ind in range(self.ny):
            foc_lin_ind = x_foc_ind*self.ny + y_foc_ind

            # Synthetic range compressed data (matched 2D filter).
            # Antenna Enclosure (lobe).
            doppler_amplitude_lin = (np.sinc(self.local_squint_ref_traj[foc_lin_ind, :]/self.radar_beamwidth*0.886 ))**2
            doppler_amplitude = np.tile(doppler_amplitude_lin, [self.simulated_image.Nt, 1])

            # Range amplitude: range positions in raw data of backscattered signal. These are the sincs with range 
            # migration (range compressed image).
            range_amplitude = np.sinc( sar_B*( (np.tile(self.simulated_image.t_axis_fast_time, [nt_fast_time, 1])).transpose()
                                               - np.tile(2*self.distances_ref_traj[foc_lin_ind, :]/c, [nt_slow_time, 1]) ) )

            # Limit bandwidth to threshold given by a window. Use only 3dB of antenna lobe for azimuth, limited by squint threshold.
            doppler_threshold_win = np.absolute( np.tile(self.local_squint_ref_traj[foc_lin_ind, :], [nt_slow_time, 1]) ) < self.squint_threshold
            raw_amplitude = doppler_amplitude*range_amplitude*doppler_threshold_win

            # Phase of backscattered signal (2*pi*2*r/lambda).
            raw_phase = np.exp(-1j*4*np.pi/sar_lambda*np.tile(self.distances_ref_traj[foc_lin_ind, :], [nt_slow_time, 1]))

            # Get module of raw_amplitude (for every xn, yn).
            mod_raw_amplitude = np.sum(abs(raw_amplitude)**2)
            # Repeat over x,y (slow time and fast time) to normalize.
            mod_raw_amplitude = np.tile(mod_raw_amplitude, [nt_slow_time, nt_fast_time])

            # Get raw odographer with raw_amplitude and raw_phase, i.e. with amplitude and phase information, and normalize.
            raw_to_foc = (np.conjugate(raw_phase))*raw_amplitude/mod_raw_amplitude

            partial_row[y_foc_ind] = np.sum(self.new_raw*raw_to_foc)

        return list(partial_row)

def generate_img(self, param):
        """
        Generate range compressed image.

        Parameters
        ----------
        param: object (ConfigurationManager instance).
          ConfigurationManager instance to read parameters from file.

        Returns
        -------
        -.
        """

        # Light speed.
        c = 300000000.0

        # Load beamwidth, bandwidth and central frequency to use locally.
        sar_bmw = (param.get_float_parameter("Radar/beamwidth")*np.pi)/180.
        sar_B = param.get_float_parameter("Radar/B")
        sar_f0 = param.get_float_parameter("Radar/f0")

        # Get angles squint and look of view with respect to the antenna coordinate system.
        #self.get_angles_antenna()
        self.local_look, self.local_squint = Utils.get_angles_antenna(self.traj, self.nom_target)

        # Set fast time axis.
        start = 2*(min(self.distances))/c - self.fast_time_pixel_margin_mono*self.radar_dt
        end = 2*(max(self.distances))/c + self.fast_time_pixel_margin_mono*self.radar_dt
        step = self.radar_dt
        self.t_axis_fast_time = np.arange(start, end, step)

        # Number of elements in fast time axis.
        self.Nt = np.size(self.t_axis_fast_time)

        self.freq_axis_fftshift = Utils.freq_axis(self.radar_dt, self.Nt, False, True)

        sar_lambda = c/sar_f0

        # Doppler amplitude (envolvente de la antena).
        doppler_amplitude = (np.sinc( (np.tile(self.local_squint, [self.Nt, 1]))/sar_bmw*(2*0.443) ))**2

        # Range amplitude: range positions in raw data of backscattered signal.
        Nd = np.size(self.distances)
        range_amplitude = np.sinc( sar_B*( (np.tile(self.t_axis_fast_time, [Nd, 1])).transpose() - np.tile(2*self.distances/c, [self.Nt, 1]) ) )

        # Signal phase received: 2*pi*2*r/lambda.
        signal_phase = np.exp(-1j*4*np.pi/sar_lambda*np.tile(self.distances, [self.Nt, 1]))

        # Generate range compressed simulated image.
        self.image = doppler_amplitude*range_amplitude*signal_phase