Python pylab.linspace() Examples

def plot_wt_layout(wt_layout, borders=None, depth=None):
    fig = plt.figure(figsize=(6,6), dpi=2000)
    fs = 14
    ax = plt.subplot(111)

    if depth is not None:
        N = 100
        X, Y = plt.meshgrid(plt.linspace(depth[:,0].min(), depth[:,0].max(), N), 
                            plt.linspace(depth[:,1].min(), depth[:,1].max(), N))
        Z = plt.griddata(depth[:,0],depth[:,1],depth[:,2],X,Y, interp='linear')
        plt.contourf(X,Y,Z, label='depth [m]')
        plt.colorbar().set_label('water depth [m]')
    #ax.plot(wt_layout.wt_positions[:,0], wt_layout.wt_positions[:,1], 'or', label='baseline position')
    
    ax.scatter(wt_layout.wt_positions[:,0], wt_layout.wt_positions[:,1], wt_layout._wt_list('rotor_diameter'), label='baseline position')

    if borders is not None:
        ax.plot(borders[:,0], borders[:,1], 'r--', label='border')
        
    ax.set_xlabel('x [m]'); 
    ax.set_ylabel('y [m]')
    ax.axis('equal');
    ax.legend(loc='lower left') 

def window_lanczos(N):
    r"""Lanczos window also known as sinc window.

    :param N: window length

    .. math:: w(n) = sinc \left(  \frac{2n}{N-1} - 1 \right)

    .. plot::
        :width: 80%
        :include-source:

        from spectrum import window_visu
        window_visu(64, 'lanczos')

    .. seealso:: :func:`create_window`, :class:`Window`
    """
    if N ==1:
        return ones(1)

    n = linspace(-N/2., N/2., N)
    win = sinc(2*n/(N-1.))
    return win 

def meyeraux(x):
    r"""Compute the Meyer auxiliary function

    The Meyer function is

    .. math:: y = 35 x^4-84 x^5+70 x^6-20 x^7

    :param array x:
    :return: the waveform

    .. plot::
        :include-source:
        :width: 80%

        from spectrum import meyeraux
        from pylab import linspace, plot
        t = linspace(0, 1, 1000)
        plot(t, meyeraux(t))

    """

    return  35*x**4-84.*x**5+70.*x**6-20.*x**7 

def window_poisson(N, alpha=2):
    r"""Poisson tapering window

    :param int N: window length

    .. math:: w(n) = \exp^{-\alpha \frac{|n|}{N/2} }

    with :math:`-N/2 \leq n \leq N/2`.

    .. plot::
        :width: 80%
        :include-source:

        from spectrum import window_visu
        window_visu(64, 'poisson')
        window_visu(64, 'poisson', alpha=3)
        window_visu(64, 'poisson', alpha=4)

    .. seealso:: :func:`create_window`, :class:`Window`
    """
    n = linspace(-N/2., (N)/2., N)
    w = exp(-alpha * abs(n)/(N/2.))
    return w 

def window_riemann(N):
    r"""Riemann tapering window

    :param int N: window length

    .. math:: w(n) = 1 - \left| \frac{n}{N/2}  \right|^2

    with :math:`-N/2 \leq n \leq N/2`.

    .. plot::
        :width: 80%
        :include-source:

        from spectrum import window_visu
        window_visu(64, 'riesz')

    .. seealso:: :func:`create_window`, :class:`Window`
    """
    n = linspace(-N/2., (N)/2., N)
    w = sin(n/float(N)*2.*pi) / (n / float(N)*2.*pi)
    return w 

def window_riesz(N):
    r"""Riesz tapering window

    :param N: window length

    .. math:: w(n) = 1 - \left| \frac{n}{N/2}  \right|^2

    with :math:`-N/2 \leq n \leq N/2`.

    .. plot::
        :width: 80%
        :include-source:

        from spectrum import window_visu
        window_visu(64, 'riesz')

    .. seealso:: :func:`create_window`, :class:`Window`
    """
    n = linspace(-N/2., (N)/2., N)

    w = 1 - abs(n/(N/2.))**2.
    return w 

def plot_window(self):
        """Plot the window in the time domain

        .. plot::
            :width: 80%
            :include-source:

            from spectrum.window import Window
            w = Window(64, name='hamming')
            w.plot_window()

        """
        from pylab import plot, xlim, grid, title, ylabel, axis
        x = linspace(0, 1, self.N)
        xlim(0, 1)
        plot(x, self.data)
        grid(True)
        title('%s Window (%s points)' % (self.name.capitalize(), self.N))
        ylabel('Amplitude')
        axis([0, 1, 0, 1.1]) 

def plot_CDR_correlation(self, doplot=True):
        """
        Displays correlation between sampling time points and CDR. It returns the two
        parameters of the linear fit, Pearson's r, p-value and standard error. If optional argument 'doplot' is
        False, the plot is not displayed.
        """
        pel2, tol = self.get_gene(self.rootlane, ignore_log=True)
        pel = numpy.array([pel2[m] for m in self.pl])*tol
        dr2 = self.get_gene('_CDR')[0]
        dr = numpy.array([dr2[m] for m in self.pl])
        po = scipy.stats.linregress(pel, dr)
        if doplot:
            pylab.scatter(pel, dr, s=9.0, alpha=0.7, c='r')
            pylab.xlim(min(pel), max(pel))
            pylab.ylim(0, max(dr)*1.1)
            pylab.xlabel(self.rootlane)
            pylab.ylabel('CDR')
            xk = pylab.linspace(min(pel), max(pel), 50)
            pylab.plot(xk, po[1]+po[0]*xk, 'k--', linewidth=2.0)
            pylab.show()
        return po 

def plot_wind_rose(wind_rose):
    fig = plt.figure(figsize=(12,5), dpi=1000)

    # Plotting the wind statistics
    ax1 = plt.subplot(121, polar=True)
    w = 2.*np.pi/len(wind_rose.frequency)
    b = ax1.bar(np.pi/2.0-np.array(wind_rose.wind_directions)/180.*np.pi - w/2.0, 
                np.array(wind_rose.frequency)*100, width=w)

    # Trick to set the right axes (by default it's not oriented as we are used to in the WE community)
    mirror = lambda d: 90.0 - d if d < 90.0 else 360.0 + (90.0 - d)
    ax1.set_xticklabels([u'%d\xb0'%(mirror(d)) for d in linspace(0.0, 360.0,9)[:-1]]);
    ax1.set_title('Wind direction frequency');

    # Plotting the Weibull A parameter
    ax2 = plt.subplot(122, polar=True)
    b = ax2.bar(np.pi/2.0-np.array(wind_rose.wind_directions)/180.*np.pi - w/2.0, 
                np.array(wind_rose.A), width=w)
    ax2.set_xticklabels([u'%d\xb0'%(mirror(d)) for d in linspace(0.0, 360.0,9)[:-1]]);
    ax2.set_title('Weibull A parameter per wind direction sectors'); 

def create_figure_ma():
    psd = test_ma()
    psd = cshift(psd, len(psd)/2) # switch positive and negative freq
    plot(linspace(-0.5, 0.5, 4096), 10 * log10(psd/max(psd)))
    ylim([-50,0])
    savefig('psd_ma.png') 

def create_figure_arma():
    psd = test_arma()
    psd = cshift(psd, len(psd)/2) # switch positive and negative freq
    plot(linspace(-0.5, 0.5, 4096), 10 * log10(psd/max(psd)))
    ylim([-50,0])
    savefig('psd_arma.png') 

def make_presentation():
    import matplotlib
    matplotlib.use('agg')
    import pylab as p
    doc = document(cache=False)

    with slide("Matplotlib figure"):
        fig = p.figure()
        x = p.linspace(0,2*p.pi)

        p.plot(x, p.sin(x), '--')

        figure(fig)


    with slide("Mpl animation"):

        anim_figs = []
        for i in range(20):
            fig = p.figure()   
            x =  p.linspace(0,2*p.pi)
            p.plot(x, p.sin(x+i))
            p.plot(x, p.sin(x+i+p.pi))
            p.close(fig) 
            anim_figs += [ fig ]


        animatesvg( anim_figs )

    with slide("Test gif"):
        figure('./test.gif', width=300)

    return doc 

def create_figure():
    res = test_minvar()
    psd = res[0]
    f = pylab.linspace(-0.5, 0.5, len(psd))
    pylab.plot(f, 10 * pylab.log10(psd/max(psd)), label='minvar 15')
    pylab.savefig('psd_minvar.png') 

def create_figure():
    psd = test_yule()
    psd = cshift(psd, len(psd)/2) # switch positive and negative freq

    plot(linspace(-0.5, 0.5, 4096),
               10 * log10(psd/max(psd)))
    axis([-0.5,0.5,-60,0])
    savefig('psd_yulewalker.png') 

def plot_rootlane_correlation(self):
        """
        Displays correlation between sampling time points and graph distance to root node. It returns the two
        parameters of the linear fit, Pearson's r, p-value and standard error.
        """
        pylab.scatter(self.pel, self.dr, s=9.0, alpha=0.7, c='r')
        pylab.xlim(min(self.pel), max(self.pel))
        pylab.ylim(0, max(self.dr)+1)
        pylab.xlabel(self.rootlane)
        pylab.ylabel('Distance to root node')
        xk = pylab.linspace(min(self.pel), max(self.pel), 50)
        pylab.plot(xk, self.po[1]+self.po[0]*xk, 'k--', linewidth=2.0)
        pylab.show()
        return self.po 

def voigt(x,c1,w1,c2,w2):
	""" Voigt function: convolution of Lorentzian and Gaussian.
		Convolution implemented with the FFT convolve function in scipy.
		NOT NORMALISED """
	
	### Create larger array so convolution doesn't screw up at the edges of the arrays
	# this assumes nicely behaved x-array...
	# i.e. x[0] == x.min() and x[-1] == x.max(), monotonically increasing
	dx = (x[-1]-x[0])/len(x)
	xp_min = x[0] - len(x)/3 * dx
	xp_max = x[-1] + len(x)/3 * dx
	xp = linspace(xp_min,xp_max,3*len(x))
	
	L = lorentzian(xp,c1,w1)
	G = gaussian(xp,c2,w2)
	
	#convolve
	V = conv(L,G,mode='same')
	
	#normalise to unity height !!! delete me later !!!
	V /= V.max()
	
	#create interpolation function to convert back to original array size
	fn_out = interp(xp,V)
	
	return fn_out(x) 

def create_figure():
    psd = test_correlog()
    f = linspace(-0.5, 0.5, len(psd))

    psd = cshift(psd, len(psd)/2)
    plot(f, 10*log10(psd/max(psd)))
    savefig('psd_corr.png') 

def create_figure():
    psd = test_eigenfre_music()
    f = linspace(-0.5, 0.5, len(psd))
    plot(f, 10 * log10(psd/max(psd)), '--',label='MUSIC 15')
    savefig('psd_eigenfre_music.png')

    psd = test_eigenfre_ev()
    f = linspace(-0.5, 0.5, len(psd))
    plot(f, 10 * log10(psd/max(psd)), '--',label='EV 15')
    savefig('psd_eigenfre_ev.png') 

def plot(self, **kargs):
        """Plot the data set, using the sampling information to set the x-axis
        correctly."""
        from pylab import plot, linspace, xlabel, ylabel, grid
        time = linspace(1*self.dt, self.N*self.dt, self.N)
        plot(time, self.data, **kargs)
        xlabel('Time')
        ylabel('Amplitude')
        grid(True) 

def contour_plot(func):
    rose = func()
    XS, YS = plt.meshgrid(np.linspace(-2, 2, 20), np.linspace(-2,2, 20));
    ZS = np.array([rose(x1=x, x2=y).f_xy for x,y in zip(XS.flatten(),YS.flatten())]).reshape(XS.shape);
    plt.contourf(XS, YS, ZS, 50);
    plt.colorbar() 

def mexican(lb, ub, n):
    r"""Generate the mexican hat wavelet

    The Mexican wavelet is:

    .. math:: w[x] = \cos{5x}  \exp^{-x^2/2}

    :param lb: lower bound
    :param ub: upper bound
    :param int n: waveform data samples
    :return: the waveform

    .. plot::
        :include-source:
        :width: 80%

        from spectrum import mexican
        from pylab import plot
        plot(mexican(0, 10, 100))

    """
    if n <= 0:
        raise ValueError("n must be strictly positive")

    x = numpy.linspace(lb, ub, n)
    psi = (1.-x**2.) * (2./(numpy.sqrt(3.)*pi**0.25)) * numpy.exp(-x**2/2.)
    return psi 

def morlet(lb, ub, n):
    r"""Generate the Morlet waveform


    The Morlet waveform is defined as follows:

    .. math:: w[x] = \cos{5x}  \exp^{-x^2/2}

    :param lb: lower bound
    :param ub: upper bound
    :param int n: waveform data samples


    .. plot::
        :include-source:
        :width: 80%

        from spectrum import morlet
        from pylab import plot
        plot(morlet(0,10,100))

    """
    if n <= 0:
        raise ValueError("n must be strictly positive")

    x = numpy.linspace(lb, ub, n)
    psi = numpy.cos(5*x) * numpy.exp(-x**2/2.)
    return psi 

def window_parzen(N):
    r"""Parsen tapering window (also known as de la Valle-Poussin)

    :param N: window length

    Parzen windows are piecewise cubic approximations
    of Gaussian windows. Parzen window sidelobes fall off as :math:`1/\omega^4`.

    if :math:`0\leq|x|\leq (N-1)/4`:

    .. math:: w(n) = 1-6 \left( \frac{|n|}{N/2} \right)^2 +6 \left( \frac{|n|}{N/2}\right)^3

    if :math:`(N-1)/4\leq|x|\leq (N-1)/2`

    .. math:: w(n) = 2 \left(1- \frac{|n|}{N/2}\right)^3


    .. plot::
        :width: 80%
        :include-source:

        from spectrum import window_visu
        window_visu(64, 'parzen')


    .. seealso:: :func:`create_window`, :class:`Window`
    """
    from numpy import  where, concatenate

    n = linspace(-(N-1)/2., (N-1)/2., N)
    n1 = n[where(abs(n)<=(N-1)/4.)[0]]
    n2 = n[where(n>(N-1)/4.)[0]]
    n3 = n[where(n<-(N-1)/4.)[0]]


    w1 = 1. -6.*(abs(n1)/(N/2.))**2 + 6*(abs(n1)/(N/2.))**3
    w2 = 2.*(1-abs(n2)/(N/2.))**3
    w3 = 2.*(1-abs(n3)/(N/2.))**3

    w = concatenate((w3, w1, w2))
    return w 

def window_gaussian(N, alpha=2.5):
    r"""Gaussian window

    :param N: window length

    .. math:: \exp^{-0.5 \left( \sigma\frac{n}{N/2} \right)^2}

    with :math:`\frac{N-1}{2}\leq n \leq \frac{N-1}{2}`.

    .. note:: N-1 is used to be in agreement with octave convention. The ENBW of
         1.4 is also in agreement with [Harris]_

    .. plot::
        :width: 80%
        :include-source:

        from spectrum import window_visu
        window_visu(64, 'gaussian', alpha=2.5)



    .. seealso:: scipy.signal.gaussian, :func:`create_window`
    """
    t = linspace(-(N-1)/2., (N-1)/2., N)
    #t = linspace(-(N)/2., (N)/2., N)
    w = exp(-0.5*(alpha * t/(N/2.))**2.)
    return w 

def _getF(self):
        if self.__response is None:
            self.compute_response()
        self.__frequencies = linspace(-0.5, 0.5, len(self.__response))
        return self.__frequencies 

def plot_range_fixed_y(self, ymin="auto", ymax="auto", ysteps=21, xmin="auto", xmax="auto", xsteps=200, clear=True, x_derivative=0):
        if ymin=="auto": ymin=self.ymin
        if ymax=="auto": ymax=self.ymax
        
        self.plot_fixed_y(_pylab.linspace(ymin, ymax, ysteps), x_derivative, xsteps, xmin, xmax, False, clear)
        _s.format_figure() 

def plot_range_fixed_x(self, xmin="auto", xmax="auto", xsteps=21, ymin="auto", ymax="auto", ysteps=200, clear=True, x_derivative=0):
        if xmin=="auto": xmin=self.xmin
        if xmax=="auto": xmax=self.xmax
        
        self.plot_fixed_x(_pylab.linspace(xmin, xmax, xsteps), x_derivative, ysteps, ymin, ymax, False, clear)
        _s.format_figure() 

def chirp(t, f0=0., t1=1., f1=100., form='linear', phase=0):
    r"""Evaluate a chirp signal at time t.

    A chirp signal is a frequency swept cosine wave.

    .. math:: a = \pi  (f_1 - f_0) / t_1
    .. math:: b = 2  \pi  f_0
    .. math:: y = \cos\left( \pi\frac{f_1-f_0}{t_1}  t^2 + 2\pi f_0 t + \rm{phase} \right)

    :param array t: times at which to evaluate the chirp signal
    :param float f0: frequency at time t=0 (Hz)
    :param float t1: time t1
    :param float f1: frequency at time t=t1 (Hz)
    :param str form: shape of frequency sweep in ['linear', 'quadratic', 'logarithmic']
    :param float phase: phase shift at t=0

    The parameter **form** can be:

        * 'linear'      :math:`f(t) = (f_1-f_0)(t/t_1) + f_0`
        * 'quadratic'   :math:`f(t) = (f_1-f_0)(t/t_1)^2 + f_0`
        * 'logarithmic' :math:`f(t) = (f_1-f_0)^{(t/t_1)} + f_0`

    Example:

    .. plot::
        :include-source:
        :width: 80%

        from spectrum import chirp
        from pylab import linspace, plot
        t = linspace(0, 1, 1000)
        y = chirp(t, form='linear')
        plot(y)
        y = chirp(t, form='quadratic')
        plot(y, 'r')

    """
    valid_forms = ['linear', 'quadratic', 'logarithmic']
    if form not in valid_forms:
        raise ValueError("Invalid form. Valid form are %s"
            % valid_forms)
    t = numpy.array(t)
    phase = 2. * pi * phase / 360.
    if form == "linear":
        a = pi * (f1 - f0)/t1
        b = 2. * pi * f0
        y = numpy.cos(a * t**2 + b*t + phase)
    elif form == "quadratic":
        a = (2/3. * pi * (f1-f0)/t1/t1)
        b = 2. * pi * f0
        y = numpy.cos(a*t**3 + b * t + phase)
    elif form == "logarithmic":
        a = 2. * pi * t1/numpy.log(f1-f0)
        b = 2. * pi * f0
        x = (f1-f0)**(1./t1)
        y = numpy.cos(a * x**t + b * t + phase)

    return y 

def window_tukey(N, r=0.5):
    """Tukey tapering window (or cosine-tapered window)

    :param N: window length
    :param r: defines the ratio between the constant section and the cosine
      section. It has to be between 0 and 1.

    The function returns a Hanning window for `r=0` and a full box for `r=1`.

    .. plot::
        :width: 80%
        :include-source:

        from spectrum import window_visu
        window_visu(64, 'tukey')
        window_visu(64, 'tukey', r=1)

    .. math:: 0.5 (1+cos(2pi/r (x-r/2))) for 0<=x<r/2

    .. math:: 0.5 (1+cos(2pi/r (x-1+r/2))) for x>=1-r/2


    .. seealso:: :func:`create_window`, :class:`Window`
    """
    assert r>=0 and r<=1 , "r must be in [0,1]"
    if N==1:
        return ones(1)

    if r == 0:
        return ones(N)
    elif r == 1:
        return window_hann(N)
    else:
        from numpy import flipud, concatenate, where

        ## cosine-tapered window
        x = linspace(0, 1, N)
        x1 = where(x<r/2.)
        w = 0.5*(1+cos(2*pi/r*(x[x1[0]]-r/2)))
        w = concatenate((w, ones(N-len(w)*2), flipud(w)))


        return w 

def create_all_psd():

    f = pylab.linspace(0, 1, 4096)

    pylab.figure(figsize=(12,8))

    # MA model
    p = spectrum.pma(xx, 64,128); p(); p.plot()
    """
    #ARMA 15 order
    a, b, rho = spectrum.arma_estimate(data, 15,15, 30)
    psd = spectrum.arma2psd(A=a,B=b, rho=rho)
    newpsd = tools.cshift(psd, len(psd)//2) # switch positive and negative freq
    pylab.plot(f, 10 * pylab.log10(newpsd/max(newpsd)), label='ARMA 15,15')
    """
    # YULE WALKER
    p = spectrum.pyule(xx, 7 , NFFT=4096, scale_by_freq=False); p.plot()
    # equivalent to
    # plot([x for x in p.frequencies()] , 10*log10(p.psd)); grid(True)

    #burg method
    p = spectrum.pburg(xx, 7, scale_by_freq=False);  p.plot()

    #pcovar
    p = spectrum.pcovar(xx, 7, scale_by_freq=False);  p.plot()

    #pmodcovar
    p = spectrum.pmodcovar(xx, 7, scale_by_freq=False); p.plot()

    # correlogram
    p = spectrum.pcorrelogram(xx, lag=60, NFFT=512, scale_by_freq=False); p.plot()

    # minvar
    p = spectrum.pminvar(xx, 7, NFFT=256, scale_by_freq=False); p.plot()

    # pmusic
    p = spectrum.pmusic(xx, 10,4, scale_by_freq=False); p.plot()

    # pmusic
    p = spectrum.pev(xx, 10, 4, scale_by_freq=False); p.plot()

    # periodogram
    p = spectrum.Periodogram(xx, scale_by_freq=False); p.plot()

    #
    legend( ["MA 32", "pyule 7", "pburg 7", "pcovar", "pmodcovar", "correlogram",
                "minvar", "pmusic", "pev", "periodgram"])


    pylab.ylim([-80,80])