Python numpy 模块,range() 实例源码

我们从Python开源项目中,提取了以下23个代码示例,用于说明如何使用numpy.range()

项目:imagepy    作者:Image-Py    | 项目源码 | 文件源码
def match(img1, img2):
    if img1.ndim == 2:
        temp = np.histogram(img1, np.arange(257))[0]
        if img2.ndim == 2:
            hist = np.histogram(img2, np.arange(257))[0]
            ahist = like(temp, hist)
            img2[:] = ahist[img2]
        if img2.ndim == 3:
            for i in range(3):
                hist = np.histogram(img2[:,:,i], np.range(257))[0]
                ahist = like(temp, hist)
                img2[:,:,i] = ahist[img2[:,:,i]]
    elif img1.ndim == 3:
        if img2.ndim == 2:
            temp = np.histogram(img1, np.arange(257))[0]
            hist = np.histogram(img2, np.arange(257))[0]
            ahist = like(temp, hist)
            img2[:] = ahist[img2]
        if img2.ndim == 3:
            for i in range(3):
                temp = np.histogram(img1[:,:,i], np.arange(257))[0]
                hist = np.histogram(img2[:,:,i], np.arange(257))[0]
                ahist = like(temp, hist)
                img2[:,:,i] = ahist[img2[:,:,i]]
项目:imagepy    作者:Image-Py    | 项目源码 | 文件源码
def run(self, ips, imgs, para = None):
        ips1 = WindowsManager.get(para['img1']).ips
        ips2 = WindowsManager.get(para['img2']).ips
        ips2.snapshot()

        img = ips1.img
        imgs = ips2.imgs

        sl1, sl2 = ips1.get_nslices(), ips2.get_nslices()
        cn1, cn2 = ips1.get_nchannels(), ips2.get_nchannels()
        if not(ips1.img.dtype == np.uint8 and ips2.img.dtype == np.uint8):
            IPy.alert('Two image must be type of 8-bit or rgb!')
            return

        for i in range(sl2):
            self.progress(i, sl2)
            match(img, imgs[i])
        ips2.update = 'pix'
项目:actinf    作者:x75    | 项目源码 | 文件源码
def rh_system_generate_sweep_input(self):
        """ActiveInferenceExperiment.rh_system_generate_sweep_input

        generate system inputs on a grid to sweep the system
        """

        # create meshgrid over proprio dimensions
        sweepsteps = 21 # 11
        dim_axes = [np.linspace(self.environment.conf.m_mins[i], self.environment.conf.m_maxs[i], sweepsteps) for i in range(self.environment.conf.m_ndims)]
        full_axes = np.meshgrid(*tuple(dim_axes), indexing='ij')

        # print "dim_axes", dim_axes
        # print "full_axes", len(full_axes)
        # print "full_axes", full_axes

        for i in range(len(full_axes)):
            print i, full_axes[i].shape
            print i, full_axes[i].flatten()

        # return proxy
        self.X_system_sweep = np.vstack([full_axes[i].flatten() for i in range(len(full_axes))]).T
项目:actinf    作者:x75    | 项目源码 | 文件源码
def rh_learn_proprio(self, iter_start = 0, iter_end = 1000):
        # hook: load_run_data
        if iter_start == 0 and os.path.exists(self.mdl_pkl):
            print "found trained model at %s, skipping learning and using that" % self.mdl_pkl
            # load data from previous run
            self.load_run_data()
            return

        # current outermost experiment loop
        for i in range(iter_start, iter_end):
            for k, v in self.rh_learn_proprio_hooks.items():
                # print "k = %s, v = %s" % (k, v)
                v(i)

    ################################################################################
    # proprio learning model variant 2 using more of prediction error
项目:actinf    作者:x75    | 项目源码 | 文件源码
def rh_model_plot(self):
        """prepare and plot model outputs over input variations from sweep"""
        assert hasattr(self, "X_model_sweep")
        assert hasattr(self, "Y_model_sweep")

        print "%s.rh_plot_model sweepsteps = %d" % (self.__class__.__name__, self.X_model_sweep.shape[0])
        print "%s.rh_plot_model environment = %s" % (self.__class__.__name__, self.environment)
        print "%s.rh_plot_model environment proprio dims = %d" % (self.__class__.__name__, self.environment.conf.m_ndims)

        # scatter_data_raw   = np.hstack((self.X_model_sweep[:,1:], self.Y_model_sweep))
        # scatter_data_cols  = ["X%d" % i for i in range(1, self.X_model_sweep.shape[1])]
        # scatter_data_cols += ["Y%d" % i for i in range(self.Y_model_sweep.shape[1])]
        # print "scatter_data_raw", scatter_data_raw.shape
        # # df = pd.DataFrame(scatter_data_raw, columns=["x_%d" % i for i in range(scatter_data_raw.shape[1])])
        # df = pd.DataFrame(scatter_data_raw, columns=scatter_data_cols)

        title = "%s, input/output sweep of model %s at time %d" % (self.mode, self.model, -1)

        # plot_scattermatrix(df)
        # plot_scattermatrix_reduced(df)
        plot_colormeshmatrix_reduced(self.X_model_sweep, self.Y_model_sweep, ymin = -1.0, ymax = 1.0, title = title)

    ################################################################################
项目:actinf    作者:x75    | 项目源码 | 文件源码
def plot_scattermatrix(df, title = "plot_scattermatrix"):
    """plot a scattermatrix of dataframe df"""
    if df is None:
        print "plot_scattermatrix: no data passed"
        return

    from pandas.tools.plotting import scatter_matrix
    # df = pd.DataFrame(X, columns=['x1_t', 'x2_t', 'x1_tptau', 'x2_tptau', 'u_t'])
    # scatter_data_raw = np.hstack((np.array(Xs), np.array(Ys)))
    # scatter_data_raw = np.hstack((Xs, Ys))
    # print "scatter_data_raw", scatter_data_raw.shape

    pl.ioff()
    # df = pd.DataFrame(scatter_data_raw, columns=["x_%d" % i for i in range(scatter_data_raw.shape[1])])
    sm = scatter_matrix(df, alpha=0.2, figsize=(10, 10), diagonal='hist')
    fig = sm[0,0].get_figure()
    fig.suptitle(title)
    if SAVEPLOTS:
        fig.savefig("fig_%03d_scattermatrix.pdf" % (fig.number), dpi=300)
    fig.show()
    # pl.show()
项目:actinf    作者:x75    | 项目源码 | 文件源码
def test_models(args):
    from actinf_models import ActInfModel
    from actinf_models import ActInfKNN
    from actinf_models import ActInfSOESGP

    idim = 4
    odim = 2
    numdatapoints = 10

    for aimclass in [ActInfModel, ActInfKNN, ActInfSOESGP, ActInfSTORKGP, ActInfGMM, ActInfHebbianSOM]:
        print("Testing aimclass = %s" % (aimclass,))
        aim = aimclass(idim = idim, odim = odim)

        X = np.random.uniform(-0.1, 0.1, (numdatapoints, 1, idim))
        y = np.random.uniform(-0.1, 0.1, (numdatapoints, 1, odim))

        for i in range(numdatapoints-1):
            print("Fitting model with X = %s, Y = %s" % (X[i].shape, y[i].shape))
            aim.fit(X[i], y[i])
        y_ = aim.predict(X[i+1])
        print("Prediction error = %s" % (y_ - y[i+1]))
项目:DL_exercises    作者:lyuwenyu    | 项目源码 | 文件源码
def main():

    model = DCGAN()

    data_load = data_loader()

    with tf.Session() as sess:

        sess.run( tf.initialize_all_variables() )
        writer = tf.train.SummaryWriter('./logs', sess.graph)

        for step in range(100):

            z, image = data_load.load_next_batch()

            feed_dict = { model.z: z, model.image: image}

            sess.run([model.d_opt, model.d_opt], feed_dict = feed_dict)

            sess.run([model.g_opt], feed_dict = feed_dict)


            writer.add_summary(summary, step)
项目:decoding-brain-challenge-2016    作者:alexandrebarachant    | 项目源码 | 文件源码
def covariances(X, estimator='cov'):
    est = _check_est(estimator)
    Nt, Ne, Ns = X.shape
    covmats = numpy.zeros((Nt, Ne, Ne))
    for i in range(Nt):
        covmats[i, :, :] = est(X[i, :, :])
    return covmats
项目:decoding-brain-challenge-2016    作者:alexandrebarachant    | 项目源码 | 文件源码
def covariances_EP(X, P, estimator='cov'):
    est = _check_est(estimator)
    Nt, Ne, Ns = X.shape
    Np, Ns = P.shape
    covmats = numpy.zeros((Nt, Ne + Np, Ne + Np))
    for i in range(Nt):
        covmats[i, :, :] = est(numpy.concatenate((P, X[i, :, :]), axis=0))
    return covmats
项目:decoding-brain-challenge-2016    作者:alexandrebarachant    | 项目源码 | 文件源码
def coherence(X, nfft=256, fs=2, noverlap=0):
    """Compute coherence."""
    n_chan = X.shape[0]
    ij = []
    for i in range(n_chan):
        for j in range(i+1, n_chan):
            ij.append((i, j))
    Cxy, Phase, freqs = mlab.cohere_pairs(X, ij, NFFT=nfft, Fs=fs,
                                          noverlap=noverlap)
    coh = numpy.zeros((n_chan, n_chan, len(freqs)))
    for i in range(n_chan):
        coh[i, i] = 1
        for j in range(i+1, n_chan):
            coh[i, j] = coh[j, i] = Cxy[(i, j)]
    return coh
项目:actinf    作者:x75    | 项目源码 | 文件源码
def rh_system_plot(self):
        """prepare and plot system outputs over input variations from sweep"""
        assert hasattr(self, "X_system_sweep")
        assert hasattr(self, "Y_system_sweep")

        print "%s.rh_plot_system sweepsteps = %d" % (self.__class__.__name__, self.X_system_sweep.shape[0])
        print "%s.rh_plot_system environment = %s" % (self.__class__.__name__, self.environment)
        print "%s.rh_plot_system environment proprio dims = %d" % (self.__class__.__name__, self.environment.conf.m_ndims)

        scatter_data_raw   = np.hstack((self.X_system_sweep, self.Y_system_sweep))
        scatter_data_cols  = ["X%d" % i for i in range(self.X_system_sweep.shape[1])]
        scatter_data_cols += ["Y%d" % i for i in range(self.Y_system_sweep.shape[1])]
        print "scatter_data_raw", scatter_data_raw.shape
        # df = pd.DataFrame(scatter_data_raw, columns=["x_%d" % i for i in range(scatter_data_raw.shape[1])])
        df = pd.DataFrame(scatter_data_raw, columns=scatter_data_cols)

        title = "%s: i/o behvaiour for %s, in = X, out = Y" % (self.mode, self.environment_str,)

        # plot_scattermatrix(df)
        plot_scattermatrix_reduced(df, title = title)

    ################################################################################
    # model sweep hooks


    # map a model
项目:actinf    作者:x75    | 项目源码 | 文件源码
def plot_scattermatrix_reduced(df, title = "plot_scattermatrix_reduced"):
    input_cols  = [i for i in df.columns if i.startswith("X")]
    output_cols = [i for i in df.columns if i.startswith("Y")]
    Xs = df[input_cols]
    Ys = df[output_cols]

    numsamples = df.shape[0]
    print "plot_scattermatrix_reduced: numsamples = %d" % numsamples

    # numplots = Xs.shape[1] * Ys.shape[1]
    # print "numplots = %d" % numplots

    gs = gridspec.GridSpec(Ys.shape[1], Xs.shape[1])
    pl.ioff()
    fig = pl.figure()
    fig.suptitle(title)
    # alpha = 1.0 / np.power(numsamples, 1.0/(Xs.shape[1] - 0))
    alpha = 0.2
    print "alpha", alpha
    cols = ["k", "b", "r", "g", "c", "m", "y"]
    for i in range(Xs.shape[1]):
        for j in range(Ys.shape[1]):
            # print "i, j", i, j, Xs, Ys
            ax = fig.add_subplot(gs[j, i])
            ax.plot(Xs.as_matrix()[:,i], Ys.as_matrix()[:,j], "ko", alpha = alpha)
            ax.set_xlabel(input_cols[i])
            ax.set_ylabel(output_cols[j])
    if SAVEPLOTS:
        fig.savefig("fig_%03d_scattermatrix_reduced.pdf" % (fig.number), dpi=300)
    fig.show()
项目:actinf    作者:x75    | 项目源码 | 文件源码
def plot_colormeshmatrix_reduced(
        X, Y, ymin = None, ymax = None,
        title = "plot_colormeshmatrix_reduced"):

    print "plot_colormeshmatrix_reduced X.shape", X.shape, "Y.shape", Y.shape
    # input_cols  = [i for i in df.columns if i.startswith("X")]
    # output_cols = [i for i in df.columns if i.startswith("Y")]
    # Xs = df[input_cols]
    # Ys = df[output_cols]

    # numsamples = df.shape[0]
    # print "plot_scattermatrix_reduced: numsamples = %d" % numsamples

    # # numplots = Xs.shape[1] * Ys.shape[1]
    # # print "numplots = %d" % numplots

    cbar_orientation = "vertical" # "horizontal"
    gs = gridspec.GridSpec(Y.shape[2], X.shape[2]/2)
    pl.ioff()
    fig = pl.figure()
    fig.suptitle(title)
    # # alpha = 1.0 / np.power(numsamples, 1.0/(Xs.shape[1] - 0))
    # alpha = 0.2
    # print "alpha", alpha
    # cols = ["k", "b", "r", "g", "c", "m", "y"]
    for i in range(X.shape[2]/2):
        for j in range(Y.shape[2]):
            # print "i, j", i, j, Xs, Ys
            ax = fig.add_subplot(gs[j, i])
            pcm = ax.pcolormesh(X[:,:,i], X[:,:,X.shape[2]/2+i], Y[:,:,j], vmin = ymin, vmax = ymax)
            # ax.plot(Xs.as_matrix()[:,i], Ys.as_matrix()[:,j], "ko", alpha = alpha)
            ax.set_xlabel("goal")
            ax.set_ylabel("error")
            cbar = fig.colorbar(mappable = pcm, ax=ax, orientation=cbar_orientation)
            ax.set_aspect(1)
    if SAVEPLOTS:
        fig.savefig("fig_%03d_colormeshmatrix_reduced.pdf" % (fig.number), dpi=300)
    fig.show()
项目:decoding_challenge_cortana_2016_3rd    作者:kingjr    | 项目源码 | 文件源码
def covariances(X, estimator='cov'):
    est = _check_est(estimator)
    Nt, Ne, Ns = X.shape
    covmats = numpy.zeros((Nt, Ne, Ne))
    for i in range(Nt):
        covmats[i, :, :] = est(X[i, :, :])
    return covmats
项目:decoding_challenge_cortana_2016_3rd    作者:kingjr    | 项目源码 | 文件源码
def covariances_EP(X, P, estimator='cov'):
    est = _check_est(estimator)
    Nt, Ne, Ns = X.shape
    Np, Ns = P.shape
    covmats = numpy.zeros((Nt, Ne + Np, Ne + Np))
    for i in range(Nt):
        covmats[i, :, :] = est(numpy.concatenate((P, X[i, :, :]), axis=0))
    return covmats
项目:decoding_challenge_cortana_2016_3rd    作者:kingjr    | 项目源码 | 文件源码
def coherence(X, nfft=256, fs=2, noverlap=0):
    """Compute coherence."""
    n_chan = X.shape[0]
    ij = []
    for i in range(n_chan):
        for j in range(i+1, n_chan):
            ij.append((i, j))
    Cxy, Phase, freqs = mlab.cohere_pairs(X, ij, NFFT=nfft, Fs=fs,
                                          noverlap=noverlap)
    coh = numpy.zeros((n_chan, n_chan, len(freqs)))
    for i in range(n_chan):
        coh[i, i] = 1
        for j in range(i+1, n_chan):
            coh[i, j] = coh[j, i] = Cxy[(i, j)]
    return coh
项目:DL_exercises    作者:lyuwenyu    | 项目源码 | 文件源码
def __init__(self, batch_size, rangex=1.0, z_dim=10,img_h=64, img_w=64, imgs_dir='./'):

        self.batch_size = batch_size
        self.img_h = img_h
        self.img_w = img_w

        self.rangex = rangex
        self.z_dim = z_dim

        #self.imgs_dir = imgs_dir
        self.images_ = glob.glob('{}/*.jpg'.format(imgs_dir))
        self.perm_ = np.random.permutation(np.range(len(self.images_)))
        self.cur_ = 0
项目:DL_exercises    作者:lyuwenyu    | 项目源码 | 文件源码
def image_loader(self):

        if self.cur_ + self.batch_size >= len(self.images_):
            self.cur_ = 0
            self.perm_ = np.random.permutation(np.range(len(self.images_)))

        img_batch = np.zeros([self.batch_size, self.img_h, self.img_w, 3])

        for i, idx in enumerate(self.perm_[self.cur_: self.cur_+self.batch_size]):
            img = transform.resize(io.imread(self.images_[idx]), (self.img_h, self.img_w))
            img_batch[i, ...] = img

        self.cur_ += self.batch_size

        return img_batch
项目:decoding-brain-challenge-2016    作者:alexandrebarachant    | 项目源码 | 文件源码
def cospectrum(X, window=128, overlap=0.75, fmin=None, fmax=None, fs=None,
               phase_correction=False):
    Ne, Ns = X.shape
    number_freqs = int(window / 2)

    step = int((1.0 - overlap) * window)
    step = max(1, step)

    number_windows = (Ns - window) / step + 1
    # pre-allocation of memory
    fdata = numpy.zeros((number_windows, Ne, number_freqs), dtype=complex)
    win = numpy.hanning(window)

    # Loop on all frequencies
    for window_ix in range(int(number_windows)):

        # time markers to select the data
        # marker of the beginning of the time window
        t1 = int(window_ix * step)
        # marker of the end of the time window
        t2 = int(t1 + window)
        # select current window and apodize it
        cdata = X[:, t1:t2] * win

        # FFT calculation
        fdata[window_ix, :, :] = numpy.fft.fft(
            cdata, n=window, axis=1)[:, 0:number_freqs]

        # if(phase_correction):
        # fdata = fdata.*(exp(-sqrt(-1)*t1*( numpy.range(window)
        # ).T/window*2*pi)*numpy.ones((1,Ne))

    # Adjust Frequency range to specified range (in case it is a parameter)
    if fmin is not None:
        f = numpy.arange(0, 1, 1.0 / number_freqs) * (fs / 2.0)
        Fix = (f >= fmin) & (f <= fmax)
        fdata = fdata[:, :, Fix]

    # fdata = fdata.real
    Nf = fdata.shape[2]
    S = numpy.zeros((Ne, Ne, Nf), dtype=complex)

    for i in range(Nf):
        S[:, :, i] = numpy.dot(
            fdata[:, :, i].conj().T, fdata[:, :, i]) / number_windows

    return S
项目:actinf    作者:x75    | 项目源码 | 文件源码
def sample_discrete_from_extero(self, i):
        self.mm.fit(self.S_ext, self.M_prop_pred)
        ext_err = np.sum(np.abs(self.goal_ext - self.S_ext))
        if i % self.goal_sample_interval == 0 or \
            ((i - self.goal_sample_time) > 5 and ext_err > 0.1):
            # update e2p
            EP = np.hstack((np.asarray(self.e2p.X_), np.asarray(self.e2p.y_)))
            # print "EP[%d] = %s" % (i, EP)
            EP = EP[10:] # knn bootstrapping creates additional datapoints
            # if i % 100 == 0:
            # re-fit gmm e2p
            # self.mm.fit(np.asarray(self.e2p.X_)[10:], np.asarray(self.e2p.y_)[10:])
            # self.mm.fit(np.asarray(self.e2p.X_)[10:], np.asarray(self.e2p.y_)[10:])

            # print "EP, cen_lst, cov_lst, p_k, logL", EP, self.cen_lst, self.cov_lst, self.p_k, self.logL
            ref_interval = 1
            self.cond = EP[(i+ref_interval) % EP.shape[0]] # X_[i,:3]
            self.cond[2:] = np.nan
            self.cond_ = np.random.uniform(-1, 1, (5, ))
            # randomly fetch an exteroceptive state that we have seen already (= reachable)
            self.goal_ext = EP[np.random.choice(range(self.numsteps/2)),:2].reshape((1, self.dim_ext))
            # self.cond_[:2] = self.goal_ext
            # self.cond_[2:] = np.nan
            # print "self.cond", self.cond
            # print "self.cond_", self.cond_

            # predict proprioceptive goal from exteroceptive one
            # if hasattr(self.mm, "cen_lst"):
            #     self.goal_prop = self.mm.sample(self.cond_)
            # else:
            #     self.goal_prop = self.mm.sample(self.goal_ext)
            self.goal_prop = self.mm.sample(self.goal_ext)

            self.goal_sample_time     = i

            # (cen_con, cov_con, new_p_k) = gmm.cond_dist(self.cond_, self.cen_lst, self.cov_lst, self.p_k)
            # self.goal_prop = gmm.sample_gaussian_mixture(cen_con, cov_con, new_p_k, samples = 1)

            # # discrete goal
            # self.goal_prop = np.random.uniform(self.environment.conf.m_mins, self.environment.conf.m_maxs, (1, self.odim))
            print "new goal_prop[%d] = %s" % (i, self.goal_prop)
            print "    goal_ext[%d] = %s" % (i, self.goal_ext)
            print "e_pred = %f" % (np.linalg.norm(self.E_prop_pred, 2))
            print "ext_er = %f" % (ext_err)

    # def lh_sample_error_gradient(self, i):
    #     """sample the local error gradient"""
    #     # hook: goal sampling
    #     if i % self.goal_sample_interval == 0:
    #         self.goal_prop = np.random.uniform(self.environment.conf.m_mins * 0.95, self.environment.conf.m_maxs * 0.95, (1, self.odim))
    #         print "new goal[%d] = %s" % (i, self.goal_prop)
    #         print "e_pred = %f" % (np.linalg.norm(self.E_prop_pred, 2))
项目:actinf    作者:x75    | 项目源码 | 文件源码
def rh_model_sweep_generate_input_grid_a(self):
        sweepsteps = 11 # 21
        # extero config
        dim_axes = [np.linspace(self.environment.conf.s_mins[i], self.environment.conf.s_maxs[i], sweepsteps) for i in range(self.environment.conf.s_ndims)]
        # dim_axes = [np.linspace(self.environment.conf.s_mins[i], self.environment.conf.s_maxs[i], sweepsteps) for i in range(self.mdl.idim)]
        print "rh_model_sweep_generate_input_grid: s_ndims = %d, dim_axes = %s" % (self.environment.conf.s_ndims, dim_axes,)
        full_axes = np.meshgrid(*tuple(dim_axes), indexing='ij')
        print "rh_model_sweep_generate_input_grid: full_axes = %s, %s" % (len(full_axes), full_axes,)

        for i in range(len(full_axes)):
            print i, full_axes[i].shape
            print i, full_axes[i].flatten()

        # return proxy
        error_grid = np.vstack([full_axes[i].flatten() for i in range(len(full_axes))])
        print "error_grid", error_grid.shape

        # draw state / goal configurations
        X_accum = []
        states = np.linspace(-1, 1, sweepsteps)
        # for state in range(1): # sweepsteps):
        for state in states:
            # randomize initial position
            # self.M_prop_pred = self.environment.compute_motor_command(np.random.uniform(-1.0, 1.0, (1, self.odim)))
            # draw random goal and keep it fixed
            # self.goal_prop = self.environment.compute_motor_command(np.random.uniform(-1.0, 1.0, (1, self.odim)))
            self.goal_prop = self.environment.compute_motor_command(np.ones((1, self.odim)) * state)
            # self.goal_prop = np.random.uniform(self.environment.conf.m_mins, self.environment.conf.m_maxs, (1, self.odim))
            GOALS = np.repeat(self.goal_prop, error_grid.shape[1], axis = 0) # as many goals as error components
            # FIXME: hacks for M1/M2
            if self.mdl.idim == 3:
                X = GOALS
            elif self.mdl.idim == 6:
                X = np.hstack((GOALS, error_grid.T))
            else:
                X = np.hstack((GOALS, error_grid.T))
            X_accum.append(X)

        X_accum = np.array(X_accum)

        # don't need this?
        # X_accum = X_accum.reshape((X_accum.shape[0] * X_accum.shape[1], X_accum.shape[2]))

        print "X_accum.shape = %s, mdl.idim = %d, mdl.odim = %d" % (X_accum.shape, self.mdl.idim, self.mdl.odim)
        # print X_accum
        X = X_accum
        # X's and pred's indices now mean: slowest: goal, e1, e2, fastest: e3
        self.X_model_sweep = X.copy()
        return X
        # print "self.X_model_sweep.shape", self.X_model_sweep.shape
项目:decoding_challenge_cortana_2016_3rd    作者:kingjr    | 项目源码 | 文件源码
def cospectrum(X, window=128, overlap=0.75, fmin=None, fmax=None, fs=None,
               phase_correction=False):
    Ne, Ns = X.shape
    number_freqs = int(window / 2)

    step = int((1.0 - overlap) * window)
    step = max(1, step)

    number_windows = (Ns - window) / step + 1
    # pre-allocation of memory
    fdata = numpy.zeros((number_windows, Ne, number_freqs), dtype=complex)
    win = numpy.hanning(window)

    # Loop on all frequencies
    for window_ix in range(int(number_windows)):

        # time markers to select the data
        # marker of the beginning of the time window
        t1 = int(window_ix * step)
        # marker of the end of the time window
        t2 = int(t1 + window)
        # select current window and apodize it
        cdata = X[:, t1:t2] * win

        # FFT calculation
        fdata[window_ix, :, :] = numpy.fft.fft(
            cdata, n=window, axis=1)[:, 0:number_freqs]

        # if(phase_correction):
        # fdata = fdata.*(exp(-sqrt(-1)*t1*( numpy.range(window)
        # ).T/window*2*pi)*numpy.ones((1,Ne))

    # Adjust Frequency range to specified range (in case it is a parameter)
    if fmin is not None:
        f = numpy.arange(0, 1, 1.0 / number_freqs) * (fs / 2.0)
        Fix = (f >= fmin) & (f <= fmax)
        fdata = fdata[:, :, Fix]

    # fdata = fdata.real
    Nf = fdata.shape[2]
    S = numpy.zeros((Ne, Ne, Nf), dtype=complex)

    for i in range(Nf):
        S[:, :, i] = numpy.dot(
            fdata[:, :, i].conj().T, fdata[:, :, i]) / number_windows

    return S