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

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

项目:decoding_challenge_cortana_2016_3rd    作者:kingjr    | 项目源码 | 文件源码
def _inside_contour(pos, contour):
    """Aux function"""
    npos = len(pos)
    x, y = pos[:, :2].T

    check_mask = np.ones((npos), dtype=bool)
    check_mask[((x < np.min(x)) | (y < np.min(y)) |
                (x > np.max(x)) | (y > np.max(y)))] = False

    critval = 0.1
    sel = np.where(check_mask)[0]
    for this_sel in sel:
        contourx = contour[:, 0] - pos[this_sel, 0]
        contoury = contour[:, 1] - pos[this_sel, 1]
        angle = np.arctan2(contoury, contourx)
        angle = np.unwrap(angle)
        total = np.sum(np.diff(angle))
        check_mask[this_sel] = np.abs(total) > critval

    return check_mask
项目:pybot    作者:spillai    | 项目源码 | 文件源码
def plot_poses(poses, pose_type='absolute'): 

    f = plt.figure(1)
    f.clf()
    f.subplots_adjust(hspace=0.25)

    # , **{'ylabel': 'Position X(m)', 'xlim': [min(xs), max(xs)]}
    pparams = { 'linewidth':1, }
    # font = { 'family': 'normal', 'weight': 'normal', 'size': 12 }
    # mpl.rc('font', **font)

    rpyxyz = np.vstack([pose.to_rpyxyz() for pose in poses])
    ax = f.add_subplot(2, 1, 1)

    if pose_type == 'absolute': 
        for j,label in enumerate(['Roll', 'Pitch', 'Yaw']): 
            ax.plot(np.unwrap(rpyxyz[:,j], discont=np.pi), label=label)
    elif pose_type == 'relative': 
        for j,label in enumerate(['Roll', 'Pitch', 'Yaw']): 
            ax.plot(np.unwrap(rpyxyz[1:,j]-rpyxyz[:-1,j], discont=np.pi), label=label)
    else: 
        raise RuntimeError('Unknown pose_type=%s, use absolute or relative' % pose_type)
    ax.legend(loc='upper right', fancybox=False, ncol=3, 
              # bbox_to_anchor=(1.0,1.2), 
              prop={'size':13})

    ax = f.add_subplot(2, 1, 2)
    for j,label in enumerate(['X', 'Y', 'Z']): 
        ax.plot(rpyxyz[:,j+3], label=label)
    ax.legend(loc='upper right', fancybox=False, ncol=3, 
              # bbox_to_anchor=(1.0,1.2), 
              prop={'size':13})

    plt.tight_layout()
    plt.show(block=True)
项目:datuner    作者:cornell-zhang    | 项目源码 | 文件源码
def another():
# plot the figure of signals
  fig = plt.figure(1)
  ax = fig.add_subplot(311)
  ax.set_title('input signal')
  ax.plot(t[1:Npts],x[1:Npts])    # just plot part of the signal
  ax = fig.add_subplot(312)
  ax.set_title('expected signal')
  ax.plot(t[1:Npts],xo[1:Npts])
  ax = fig.add_subplot(313)
  ax.set_title('filtered signal')
  ax.plot(t[1:Npts],y[1:Npts])
  fig.savefig('pic1.png')

# plot the mag & phase response of the LPF
  w,h = signal.freqz(b,1)
  hdb = 20 * np.log10(np.abs(h))
  hphs = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
  fig2=plt.figure(2)
  ax2 = fig2.add_subplot(211)
  ax2.set_title('frequency response')
  ax2.plot(w/max(w),hdb)
  ax2 = fig2.add_subplot(212)
  ax2.set_title('phase response')
  ax2.plot(w/max(w),hphs)
  fig2.savefig('pic2.png')
项目:jiveplot    作者:haavee    | 项目源码 | 文件源码
def unwrap(plotar, ms2mappings):
    for k in plotar.keys():
        for d in plotar[k].keys():
            # get a reference to the data set
            dsref = plotar[k][d]
            # get sorted data and unwrap the phases
            (xvals, yvals) = zip(*sorted(zip(dsref.xval, dsref.yval), key=operator.itemgetter(0)))
            yvals      = numpy.unwrap(numpy.deg2rad(yvals))
            dsref.xval = numpy.array(xvals)
            dsref.yval = numpy.rad2deg(yvals)
项目:jiveplot    作者:haavee    | 项目源码 | 文件源码
def phaserate(plotar, ms2mappings):
    if plotar.plotType!='phatime':
        raise RuntimeError, "phaserate() cannot run on plot type {0}".format( plotar.plotType )
    spm = ms2mappings.spectralMap
    # iterate over all plots and all data sets within the plots
    for k in plotar.keys():
        for d in plotar[k].keys():
            # get a reference to the data set
            dsref = plotar[k][d]
            # get the full data set label - we have access to all the data set's properties (FQ, SB, POL etc)
            n     = plots.join_label(k, d)
            # fit a line through the unwrapped phase 
            unw    = numpy.unwrap(numpy.deg2rad(dsref.yval))
            coeffs = numpy.polyfit(dsref.xval, unw, 1)
            # evaluate the fitted polynomial at the x-loci
            extray = numpy.polyval(coeffs, dsref.xval)
            # here we could compute the reliability of the fit
            diff   = unw - extray
            ss_tot = numpy.sum(numpy.square(unw - unw.mean()))
            ss_res = numpy.sum(numpy.square(diff))
            r_sq   = 1.0 - ss_res/ss_tot
            # compare std deviation and variance in the residuals after fit
            std_r  = numpy.std(diff)
            var_r  = numpy.var(diff)
            f      = spm.frequencyOfFREQ_SB(n.FQ, n.SB)
            rate   = coeffs[0]
            if var_r<std_r:
                print "{0}: {1:.8f} ps/s @ {2:5.4f}MHz [R2={3:.3f}]".format(n, rate/(2.0*numpy.pi*f*1.0e-12), f/1.0e6, r_sq )
                # before plotting wrap back to -pi,pi and transform to degrees
                dsref.extra = [ drawline_fn(n, dsref.xval, numpy.rad2deg(do_wrap(extray))) ]
项目:jiveplot    作者:haavee    | 项目源码 | 文件源码
def phasedbg(plotar, ms2mappings):
    global cache
    if plotar.plotType!='phatime':
        raise RuntimeError, "phasedbg() cannot run on plot type {0}".format( plotar.plotType )
    store = len(cache)==0
    # iterate over all plots and all data sets within the plots
    for k in plotar.keys():
        for d in plotar[k].keys():
            # get a reference to the data set
            dsref = plotar[k][d]
            # get the full data set label - we have access to all the data set's properties (FQ, SB, POL etc)
            n     = plots.join_label(k, d)
            # fit a line through the unwrapped phase 
            unw    = numpy.unwrap(numpy.deg2rad(dsref.yval))
            #coeffs = numpy.polyfit(dsref.xval, unw, 1)
            coeffs = numpy.polyfit(xrange(len(dsref.yval)), unw, 1)
            # evaluate the fitted polynomial at the x-loci
            extray = numpy.polyval(coeffs, dsref.xval)
            # here we could compute the reliability of the fit
            diff   = unw - extray
            # compare std deviation and variance in the residuals after fit
            std_r  = numpy.std(diff)
            var_r  = numpy.var(diff)
            coeffs = numpy.rad2deg(coeffs)
            if var_r<std_r:
                # decide what to do
                if store:
                    cache[n] = coeffs
                else:
                    # check if current key exists in cache; if so
                    # do differencing
                    otherCoeffs = cache.get(n, None)
                    if otherCoeffs is None:
                        print "{0}: not found in cache".format( n )
                    else:
                        delta = otherCoeffs - coeffs
                        print "{0.BL} {0.SB} {0.P}: dRate={1:5.4f} dOff={2:4.1f}".format(n, delta[0], delta[1]+360.0 if delta[1]<0.0 else delta[1])
                # before plotting wrap back to -pi,pi and transform to degrees
                #dsref.extra = [ drawline_fn(n, dsref.xval, numpy.rad2deg(do_wrap(extray))) ]
    if not store:
        cache = {}
项目:arlpy    作者:org-arl    | 项目源码 | 文件源码
def freqz(b, a=1, fs=2.0, worN=None, whole=False):
    """Plot frequency response of a filter.

    This is a convenience function to plot frequency response, and internally uses
    :func:`scipy.signal.freqz` to estimate the response. For further details, see the
    documentation for :func:`scipy.signal.freqz`.

    :param b: numerator of a linear filter
    :param a: denominator of a linear filter
    :param fs: sampling rate in Hz (optional, normalized frequency if not specified)
    :param worN: see :func:`scipy.signal.freqz`
    :param whole: see :func:`scipy.signal.freqz`
    :returns: (frequency vector, frequency response vector)

    >>> import arlpy
    >>> arlpy.signal.freqz([1,1,1,1,1], fs=120000);
    >>> w, h = arlpy.signal.freqz([1,1,1,1,1], fs=120000)
    """
    import matplotlib.pyplot as plt
    w, h = _sig.freqz(b, a, worN, whole)
    f = w*fs/(2*_np.pi)
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    plt.plot(f, 20*_np.log10(abs(h)), 'b')
    plt.ylabel('Amplitude [dB]', color='b')
    plt.xlabel('Frequency [Hz]')
    plt.grid()
    ax1.twinx()
    angles = _np.unwrap(_np.angle(h))
    plt.plot(f, angles, 'g')
    plt.ylabel('Angle (radians)', color='g')
    plt.axis('tight')
    plt.show()
    return w, h
项目:cebl    作者:idfah    | 项目源码 | 文件源码
def plotPhaseResponse(self, freqs=None, scale='radians',
                         showCorners=True,
                         label='Frequency Response',
                         ax=None, **kwargs):
        """Plot the frequency response of the filter.
        """
        if ax is None:
            fig = plt.figure()
            ax = fig.add_subplot(1,1,1)

        freqs, responses = self.frequencyResponse(freqs=freqs)
        freqs = freqs * self.sampRate * 0.5 / np.pi
        responseAngles = np.unwrap(np.angle(responses))

        scale = scale.lower()
        if scale == 'radians':
            ax.set_ylabel('Phase (Radians)')
        elif scale == 'cycles':
            ax.set_ylabel('Phase (Cycles)')
            responseAngles /= 2.0*np.pi
        elif scale == 'degrees':
            ax.set_ylabel('Phase (Degrees)')
            responseAngles = responseAngles*180.0 / np.pi
        else:
            raise Exception('Invalid scale: ' + str(scale) + '.')

        lines = ax.plot(freqs, responseAngles,
                        label=label, **kwargs)

        result = {'ax': ax, 'lines': lines}

        if showCorners:
            cornerLines = ax.vlines((self.lowFreq,self.highFreq),
                np.min(responseAngles), np.max(responseAngles),
                color='violet', linestyle='--', label='Corners')
            result['cornerLines'] = cornerLines

        ax.set_xlabel('Frequency (Hz)')

        return result
项目:CRIkit2    作者:CoherentRamanNIST    | 项目源码 | 文件源码
def _calc(self, data, ret_obj, **kwargs):

        self._inst_als = _AlsCvxopt(**kwargs)

        try:
            shp = data.shape[0:-2]
            total_num = _np.array(shp).prod()

            counter = 1
            for idx in _np.ndindex(shp):
                print('Detrended iteration {} / {}'.format(counter, total_num))
                ph = _np.unwrap(_np.angle(data[idx]))
                if self.rng is None:
                    err_phase = self._inst_als.calculate(ph)
                else:
                    err_phase = self._inst_als.calculate(ph[..., self.rng])

                h = _np.zeros(err_phase.shape)
                h += _hilbert(err_phase)

                correction_factor = 1/_np.exp(h) * _np.exp(-1j*err_phase)

                if self.rng is None:
                    ret_obj[idx] *= correction_factor
                else:
                    ret_obj[idx][..., self.rng] *= correction_factor
                counter += 1
        except:
            return False
        else:
#            print(self._inst_als.__dict__)
            return True
项目:zignal    作者:ronnyandersson    | 项目源码 | 文件源码
def phase_resp(self, frequencies=None, unwrap=False):
        """Calculate the real phase response"""
        w, h    = self.complex_freq_resp(frequencies)
        phase   = np.angle(h, deg=False)
        phase   = np.unwrap(phase) if unwrap else phase
        phase   = np.rad2deg(phase)
        freqs   = rad2hz(w, self.fs)

        return freqs, phase
项目:zignal    作者:ronnyandersson    | 项目源码 | 文件源码
def plot_mag_phase(self, filename=None, plotpoints=10000, unwrap=False):
        """Produce a plot with magnitude and phase response in the same
        figure. The y-axis on the left side belongs to the magnitude
        response. The y-axis on the right side belongs to the phase
        response
        """
        unused_freq, mag = self.magnitude_resp(plotpoints)
        freq, pha = self.phase_resp(plotpoints, unwrap=unwrap)

        fig = plt.figure(1)
        ax_mag = fig.add_subplot(111)
        ax_pha = ax_mag.twinx()
        ax_mag.semilogx(freq, mag, label='magnitude', color='red',  linestyle='-')
        ax_pha.semilogx(freq, pha, label='phase',     color='blue', linestyle='--')
        ax_mag.grid(True)

        ax_mag.set_xlim(10, self.fs/2)
        #ax_mag.set_ylim(bottom=-80)    # FIXME: ad proper padding
        #ax_mag.margins(0.1)

        ax_mag.set_title('Frequency response')
        ax_mag.set_xlabel('Frequency [Hz]')
        ax_mag.set_ylabel('Magnitude [dB]')
        ax_pha.set_ylabel('Phase [deg]')

        handles1, labels1 = ax_mag.get_legend_handles_labels()
        handles2, labels2 = ax_pha.get_legend_handles_labels()

        plt.legend(handles1 + handles2, labels1 + labels2, loc='best')

        if filename is None:
            plt.show()
        else:
            try:
                plt.savefig(filename)
            finally:
                plt.close(1)
项目:stlab    作者:yausern    | 项目源码 | 文件源码
def phaseunwrap(array):
    """
    unwraps the phase from a given complex array returning a signal with 0 average phase slope
    """
    data = np.asarray(array)
    phase = np.unwrap(np.angle(data))
    avg = np.average(np.diff(phase))
    data = [x*np.exp(-1j*avg*i) for i,x in enumerate(data)]
#    print(np.average(np.diff(np.angle(data))))
    return np.asarray(data)
项目:stlab    作者:yausern    | 项目源码 | 文件源码
def find_resonance(frec,S11,margin=31,doplots=False):
    """
    Returns the resonance frequency from maximum of the derivative of the phase
    """
    #Smooth data for initial guesses
    sReS11 = np.array(smooth(S11.real,margin,3))
    sImS11 = np.array(smooth(S11.imag,margin,3))
    sS11 = np.array( [x+1j*y for x,y in zip(sReS11,sImS11) ] )
    #Make smoothed phase vector removing 2pi jumps
    sArgS11 = np.angle(sS11)
    sArgS11 = np.unwrap(sArgS11)
    sdiffang = np.diff(sArgS11)

    #Get resonance index from maximum of the derivative of the phase
    avgph =  np.average(sdiffang)
    errvec = [np.power(x-avgph,2.) for x in sdiffang]
    ires = np.argmax(errvec[margin:-margin])+margin
    f0i=frec[ires]
    print("Max index: ",ires," Max frequency: ",f0i)

    if doplots:
        plt.title('Original signal (Re,Im)')
        plt.plot(frec,S11.real)
        plt.plot(frec,S11.imag)
        plt.show()

        plt.plot(np.angle(sS11))
        plt.title('Smoothed Phase')
        plt.axis('auto')
        plt.show()
        plt.plot(sdiffang[margin:-margin])
        plt.plot(sdiffang)
        plt.title('Diff of Smoothed Phase')
        plt.show()

        plt.title('Smoothed (ph-phavg)\^2')
        plt.plot(errvec)
        plt.plot([ires],[errvec[ires]],'ro')
        plt.show()

    return f0i
项目:stlab    作者:yausern    | 项目源码 | 文件源码
def diff_find_resonance(frec,diffS11,margin=31,doplots=False):
    """
    Returns the resonance frequency from maximum of the derivative of the phase
    """
    #Smooth data for initial guesses
    sReS11 = np.array(smooth(diffS11.real,margin,3))
    sImS11 = np.array(smooth(diffS11.imag,margin,3))
    sS11 = np.array( [x+1j*y for x,y in zip(sReS11,sImS11) ] )
    #Make smoothed phase vector removing 2pi jumps
    sArgS11 = np.angle(sS11)
    sArgS11 = np.unwrap(sArgS11)
    sdiffang = np.diff(sArgS11)

    #Get resonance index from maximum of the derivative of the phase
    avgph =  np.average(sdiffang)
    errvec = [np.power(x-avgph,2.) for x in sdiffang]
    ires = np.argmax(errvec[margin:-margin])+margin
    f0i=frec[ires]
    print("Max index: ",ires," Max frequency: ",f0i)

    if doplots:
        plt.title('Original signal (Re,Im)')
        plt.plot(frec,diffS11.real)
        plt.plot(frec,diffS11.imag)
        plt.plot(frec,np.abs(diffS11))
        plt.show()

        plt.plot(np.angle(sS11))
        plt.title('Smoothed Phase')
        plt.axis('auto')
        plt.show()
        plt.plot(sdiffang[margin:-margin])
        plt.title('Diff of Smoothed Phase')
        plt.show()

        plt.title('Smoothed (ph-phavg)\^2')
        plt.plot(errvec)
        plt.plot([ires],[errvec[ires]],'ro')
        plt.show()

    return f0i
项目:orbs    作者:thomasorb    | 项目源码 | 文件源码
def unwrap_phase_map_0(self):
        """Unwrap order 0 phase map.


        Phase is defined modulo pi/2. The Unwrapping is a
        reconstruction of the phase so that the distance between two
        neighboor pixels is always less than pi/4. Then the real phase
        pattern can be recovered and fitted easily.

        The idea is the same as with np.unwrap() but in 2D, on a
        possibly very noisy map, where a naive 2d unwrapping cannot be
        done.
        """
        self.phase_map_order_0_unwraped = orb.utils.image.unwrap_phase_map0(
            np.copy(self.phase_maps[0]))

        # Save unwraped map
        phase_map_path = self._get_phase_map_path(0, phase_map_type='unwraped')

        self.write_fits(phase_map_path,
                        orb.cutils.unbin_image(
                            np.copy(self.phase_map_order_0_unwraped),
                            self.dimx_unbinned,
                            self.dimy_unbinned), 
                        fits_header=self._get_phase_map_header(
                            0, phase_map_type='unwraped'),
                        overwrite=self.overwrite)
        if self.indexer is not None:
            self.indexer['phase_map_unwraped_0'] = phase_map_path
项目:LabelFusion    作者:RobotLocomotion    | 项目源码 | 文件源码
def main():

    filename = 'posegraph.posegraph'
    plot = False
    windowLength = 11
    windowType = 'hanning'

    if not os.path.isfile(filename):
        print 'could not find posegraph file:', filename
        sys.exit(1)

    data = np.loadtxt(filename)
    poseTimes = data[:,0]
    poses = np.array(data[:,1:])
    poses = np.array([to_xyzrpy(pose) for pose in poses])
    poses[:,3:] = np.unwrap(poses[:,3:])

    for i in range(poses.shape[1]):
        poses[:,i] = smooth(poses[:,i], window_len=windowLength, window=windowType)

    if plot:
        plot(poses, poseTimes)

    poses = np.array([to_xyzquat(pose) for pose in poses])
    poses = np.hstack([np.reshape(poseTimes, (-1,1)), poses])
    np.savetxt('posegraph_smoothed.posegraph', poses)
项目:GRIPy    作者:giruenf    | 项目源码 | 文件源码
def __init__(self, real_signal, sampling):
        self._fs = 1/sampling
        self._analytic_signal =  hilbert(real_signal)           
        self._amplitude_envelope = np.abs(self._analytic_signal)
        self._instantaneous_phase = np.unwrap(np.angle(self._analytic_signal))
        self._instantaneous_frequency = (np.diff(self._instantaneous_phase) / 
                                        (2.0*np.pi) * self._fs)
        self._instantaneous_frequency = np.insert(self._instantaneous_frequency, 0, np.nan)
项目:radar    作者:amoose136    | 项目源码 | 文件源码
def unwrap(p, discont=pi, axis=-1):
    """
    Unwrap by changing deltas between values to 2*pi complement.

    Unwrap radian phase `p` by changing absolute jumps greater than
    `discont` to their 2*pi complement along the given axis.

    Parameters
    ----------
    p : array_like
        Input array.
    discont : float, optional
        Maximum discontinuity between values, default is ``pi``.
    axis : int, optional
        Axis along which unwrap will operate, default is the last axis.

    Returns
    -------
    out : ndarray
        Output array.

    See Also
    --------
    rad2deg, deg2rad

    Notes
    -----
    If the discontinuity in `p` is smaller than ``pi``, but larger than
    `discont`, no unwrapping is done because taking the 2*pi complement
    would only make the discontinuity larger.

    Examples
    --------
    >>> phase = np.linspace(0, np.pi, num=5)
    >>> phase[3:] += np.pi
    >>> phase
    array([ 0.        ,  0.78539816,  1.57079633,  5.49778714,  6.28318531])
    >>> np.unwrap(phase)
    array([ 0.        ,  0.78539816,  1.57079633, -0.78539816,  0.        ])

    """
    p = asarray(p)
    nd = len(p.shape)
    dd = diff(p, axis=axis)
    slice1 = [slice(None, None)]*nd     # full slices
    slice1[axis] = slice(1, None)
    ddmod = mod(dd + pi, 2*pi) - pi
    _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
    ph_correct = ddmod - dd
    _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
    up = array(p, copy=True, dtype='d')
    up[slice1] = p[slice1] + ph_correct.cumsum(axis)
    return up
项目:krpcScripts    作者:jwvanderbeck    | 项目源码 | 文件源码
def unwrap(p, discont=pi, axis=-1):
    """
    Unwrap by changing deltas between values to 2*pi complement.

    Unwrap radian phase `p` by changing absolute jumps greater than
    `discont` to their 2*pi complement along the given axis.

    Parameters
    ----------
    p : array_like
        Input array.
    discont : float, optional
        Maximum discontinuity between values, default is ``pi``.
    axis : int, optional
        Axis along which unwrap will operate, default is the last axis.

    Returns
    -------
    out : ndarray
        Output array.

    See Also
    --------
    rad2deg, deg2rad

    Notes
    -----
    If the discontinuity in `p` is smaller than ``pi``, but larger than
    `discont`, no unwrapping is done because taking the 2*pi complement
    would only make the discontinuity larger.

    Examples
    --------
    >>> phase = np.linspace(0, np.pi, num=5)
    >>> phase[3:] += np.pi
    >>> phase
    array([ 0.        ,  0.78539816,  1.57079633,  5.49778714,  6.28318531])
    >>> np.unwrap(phase)
    array([ 0.        ,  0.78539816,  1.57079633, -0.78539816,  0.        ])

    """
    p = asarray(p)
    nd = len(p.shape)
    dd = diff(p, axis=axis)
    slice1 = [slice(None, None)]*nd     # full slices
    slice1[axis] = slice(1, None)
    ddmod = mod(dd + pi, 2*pi) - pi
    _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
    ph_correct = ddmod - dd
    _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
    up = array(p, copy=True, dtype='d')
    up[slice1] = p[slice1] + ph_correct.cumsum(axis)
    return up
项目:dyfunconn    作者:makism    | 项目源码 | 文件源码
def analytic_signal(signal, fb, fs=128, order=3):
    """ Passband filtering and Hilbert transformation


    Parameters
    ----------
    signal: real array-like, shape(n_channels, n_samples)
        Input signal

    fb: list of length 2
        The low and high frequencies

    fs: int
        Sampling frequency

    order : int
        Filter order

    Returns
    -------
    filtered_signal: real array-like, shape(n_channels, n_samples)
        The input signal, filtered within the given frequencies

    hilberted_signal: complex array-like, shape(n_channels, n_samples)
        The Hilbert representation of the input signal

    unwrapped_phase: real array-like, shape(n_channels, n_samples)
        The unwrapped phase of the Hilbert representation


    Notes
    -----
    Internally, we use SciPy's Butterworth implementation (`scipy.signal.butter`)
    and the two-pass filter `scipy.signal.filtfilt` to achieve results identical
    to MATLAB.
    """
    fs = float(fs)

    passband = [fb[0] / (fs / 2.0), fb[1] / (fs / 2.0)]
    passband = np.ravel(passband)
    b, a = scipy.signal.butter(
        order, passband, 'bandpass', analog=False, output='ba')

    filtered_signal = scipy.signal.filtfilt(b, a, signal)
    hilberted_signal = scipy.signal.hilbert(filtered_signal)
    unwrapped_phase = np.unwrap(np.angle(hilberted_signal))

    return (filtered_signal, hilberted_signal, unwrapped_phase)
项目:pyMTL    作者:bibliolytic    | 项目源码 | 文件源码
def test_temporal_learning(flen=[5]):
    '''
    Function that computes CSP filters then uses those with the temporal filter MTL 
    idea, and confirms that the output has a spectral profile that is similar to expected.
    Generate y values from trace of filtered covariance to ensure that is not an issue
    '''
    def covmat(x,y):
        return (1/(x.shape[1]-1)*x.dot(y.T))
    D = scio.loadmat('/is/ei/vjayaram/code/python/pyMTL_MunichMI.mat')
    data = D['T'].ravel()
    labels = D['Y'].ravel()
    fig = plt.figure()
    fig.suptitle('Recovered frequency filters for various filter lengths')
    model = TemporalBRC(max_prior_iter=100)
    # compute cross-covariance matrices and generate X
    for find,freq in enumerate(flen):
        X = []
        for tind,d in enumerate(data):
            d = d.transpose((2,0,1))
            x = np.zeros((d.shape[0], freq))
            nsamples = d.shape[2]
            for ind, t in enumerate(d):
                for j in range(freq):
                    C = covmat(t[:,0:(nsamples-j)],t[:,j:])
                    x[ind,j] = np.trace(C + C.T)
            X.append(x)
            labels[tind] = labels[tind].ravel()

        model.fit_multi_task(X,labels)
        b = numerics.solve_fir_coef(model.prior.mu.flatten())[0]
        # look at filter properties
        w,h = freqz(b[1:],worN=100)
        w = w*500/2/np.pi # convert to Hz

        ax1 = fig.add_subplot(3,3,find+1)
        plt.plot(w, 20 * np.log10(abs(h)), 'b')
        plt.ylabel('Amplitude [dB]', color='b')
        plt.xlabel('Frequency [Hz]')
        ax2 = ax1.twinx()
        angles = np.unwrap(np.angle(h))
        plt.plot(w, angles, 'g')
        plt.ylabel('Angle (radians)', color='g')
        plt.grid()
        plt.axis('tight')
    plt.show()
项目:stlab    作者:yausern    | 项目源码 | 文件源码
def GetAllData(self,keep_uncal=True):
        pars,parnames = self.GetTraceNames()
        self.SetActiveTrace(pars[0])
        names = ['Frequency (Hz)']
        alltrc = [self.GetFrequency()]
        for pp in parnames:
            names.append('%sre ()' % pp)
            names.append('%sim ()' % pp)
            names.append('%sdB (dB)' % pp)
            names.append('%sPh (rad)' % pp)
        for par in pars:
            self.SetActiveTrace(par)
            yyre,yyim = self.GetTraceData()
            alltrc.append(yyre)
            alltrc.append(yyim)
            yydb = 20.*np.log10( np.abs(yyre+1j*yyim) )
            yyph = np.unwrap(np.angle(yyre+1j*yyim))
            alltrc.append(yydb)
            alltrc.append(yyph)
        Cal = self.GetCal()
        if Cal and keep_uncal:
            for pp in parnames:
                names.append('%sre unc ()' % pp)
                names.append('%sim unc ()' % pp)
                names.append('%sdB unc (dB)' % pp)
                names.append('%sPh unc (rad)' % pp)
            self.CalOff()
            for par in pars:
                self.SetActiveTrace(par)
                yyre,yyim = self.GetTraceData()
                alltrc.append(yyre)
                alltrc.append(yyim)
                yydb = 20.*np.log10( np.abs(yyre+1j*yyim) )
                yyph = np.unwrap(np.angle(yyre+1j*yyim))
                alltrc.append(yydb)
                alltrc.append(yyph)
            self.CalOn()
        final = stlabdict()
        for name,data in zip(names,alltrc):
            final[name]=data
        final.addparcolumn('Power (dBm)',self.GetPower())
        return final
项目:aws-lambda-numpy    作者:vitolimandibhrata    | 项目源码 | 文件源码
def unwrap(p, discont=pi, axis=-1):
    """
    Unwrap by changing deltas between values to 2*pi complement.

    Unwrap radian phase `p` by changing absolute jumps greater than
    `discont` to their 2*pi complement along the given axis.

    Parameters
    ----------
    p : array_like
        Input array.
    discont : float, optional
        Maximum discontinuity between values, default is ``pi``.
    axis : int, optional
        Axis along which unwrap will operate, default is the last axis.

    Returns
    -------
    out : ndarray
        Output array.

    See Also
    --------
    rad2deg, deg2rad

    Notes
    -----
    If the discontinuity in `p` is smaller than ``pi``, but larger than
    `discont`, no unwrapping is done because taking the 2*pi complement
    would only make the discontinuity larger.

    Examples
    --------
    >>> phase = np.linspace(0, np.pi, num=5)
    >>> phase[3:] += np.pi
    >>> phase
    array([ 0.        ,  0.78539816,  1.57079633,  5.49778714,  6.28318531])
    >>> np.unwrap(phase)
    array([ 0.        ,  0.78539816,  1.57079633, -0.78539816,  0.        ])

    """
    p = asarray(p)
    nd = len(p.shape)
    dd = diff(p, axis=axis)
    slice1 = [slice(None, None)]*nd     # full slices
    slice1[axis] = slice(1, None)
    ddmod = mod(dd + pi, 2*pi) - pi
    _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
    ph_correct = ddmod - dd
    _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
    up = array(p, copy=True, dtype='d')
    up[slice1] = p[slice1] + ph_correct.cumsum(axis)
    return up
项目:lambda-numba    作者:rlhotovy    | 项目源码 | 文件源码
def unwrap(p, discont=pi, axis=-1):
    """
    Unwrap by changing deltas between values to 2*pi complement.

    Unwrap radian phase `p` by changing absolute jumps greater than
    `discont` to their 2*pi complement along the given axis.

    Parameters
    ----------
    p : array_like
        Input array.
    discont : float, optional
        Maximum discontinuity between values, default is ``pi``.
    axis : int, optional
        Axis along which unwrap will operate, default is the last axis.

    Returns
    -------
    out : ndarray
        Output array.

    See Also
    --------
    rad2deg, deg2rad

    Notes
    -----
    If the discontinuity in `p` is smaller than ``pi``, but larger than
    `discont`, no unwrapping is done because taking the 2*pi complement
    would only make the discontinuity larger.

    Examples
    --------
    >>> phase = np.linspace(0, np.pi, num=5)
    >>> phase[3:] += np.pi
    >>> phase
    array([ 0.        ,  0.78539816,  1.57079633,  5.49778714,  6.28318531])
    >>> np.unwrap(phase)
    array([ 0.        ,  0.78539816,  1.57079633, -0.78539816,  0.        ])

    """
    p = asarray(p)
    nd = len(p.shape)
    dd = diff(p, axis=axis)
    slice1 = [slice(None, None)]*nd     # full slices
    slice1[axis] = slice(1, None)
    ddmod = mod(dd + pi, 2*pi) - pi
    _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
    ph_correct = ddmod - dd
    _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
    up = array(p, copy=True, dtype='d')
    up[slice1] = p[slice1] + ph_correct.cumsum(axis)
    return up
项目:deliver    作者:orchestor    | 项目源码 | 文件源码
def unwrap(p, discont=pi, axis=-1):
    """
    Unwrap by changing deltas between values to 2*pi complement.

    Unwrap radian phase `p` by changing absolute jumps greater than
    `discont` to their 2*pi complement along the given axis.

    Parameters
    ----------
    p : array_like
        Input array.
    discont : float, optional
        Maximum discontinuity between values, default is ``pi``.
    axis : int, optional
        Axis along which unwrap will operate, default is the last axis.

    Returns
    -------
    out : ndarray
        Output array.

    See Also
    --------
    rad2deg, deg2rad

    Notes
    -----
    If the discontinuity in `p` is smaller than ``pi``, but larger than
    `discont`, no unwrapping is done because taking the 2*pi complement
    would only make the discontinuity larger.

    Examples
    --------
    >>> phase = np.linspace(0, np.pi, num=5)
    >>> phase[3:] += np.pi
    >>> phase
    array([ 0.        ,  0.78539816,  1.57079633,  5.49778714,  6.28318531])
    >>> np.unwrap(phase)
    array([ 0.        ,  0.78539816,  1.57079633, -0.78539816,  0.        ])

    """
    p = asarray(p)
    nd = len(p.shape)
    dd = diff(p, axis=axis)
    slice1 = [slice(None, None)]*nd     # full slices
    slice1[axis] = slice(1, None)
    ddmod = mod(dd + pi, 2*pi) - pi
    _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
    ph_correct = ddmod - dd
    _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
    up = array(p, copy=True, dtype='d')
    up[slice1] = p[slice1] + ph_correct.cumsum(axis)
    return up
项目:GRIPy    作者:giruenf    | 项目源码 | 文件源码
def STFT(x, window_size, noverlap, time_start, Ts, mode='psd'):
    print 'Ts:', Ts
    f0 = 1.0/Ts
    print 'F Zero:', f0
    print 'Frquencia de Nyquist:', f0/2, 'Hz'
    print 'Resolução em frquencia:', f0/window_size, 'Hz'    
    print 'Resolução temporal:', Ts, 'seconds'   

    #
    window = scipy.hanning(window_size)
    stft_data = np.array([np.fft.fft(window*data) 
                for data in SlidingWindow(x, len(window), noverlap)]
    )
    #




    if len(window) % 2:
        numFreqs = stft_data.shape[1]/2+1
    else:
        numFreqs = stft_data.shape[1]/2
    freqs = np.fft.fftfreq(window_size, Ts)[:numFreqs]
    stft_data = stft_data[:, :numFreqs]
    time_values = np.arange(window_size/2, 
                         len(x)-window_size/2+1, 
                         window_size-noverlap) * Ts  + time_start
    #                        
    if mode == 'psd':
        # Also include scaling factors for one-sided densities and dividing by
        # the sampling frequency, if desired. Scale everything, except the DC
        # component and the NFFT/2 component:
        stft_data *= 2.
        # MATLAB divides by the sampling frequency so that density function
        # has units of dB/Hz and can be integrated by the plotted frequency
        # values. Perform the same scaling here.
        #if scale_by_freq:
        scale_by_freq = True    
        if scale_by_freq:    
            Fs = 1/Ts    
            stft_data /= Fs
            # Scale the spectrum by the norm of the window to compensate for
            # windowing loss; see Bendat & Piersol Sec 11.5.2.
            stft_data /= (np.abs(window)**2).sum()   
        else:
            # In this case, preserve power in the segment, not amplitude
            stft_data /= np.abs(window).sum()**2
        stft_data = stft_data * np.conjugate(stft_data) 
        stft_data = stft_data.real 

    elif mode == 'magnitude':
        stft_data = np.absolute(stft_data)   

    elif mode == 'angle' or mode == 'phase':
        stft_data = np.angle(stft_data)
        if mode == 'phase':
            stft_data = np.unwrap(stft_data, axis=0)

    return stft_data.T, freqs, time_values
项目:Alfred    作者:jkachhadia    | 项目源码 | 文件源码
def unwrap(p, discont=pi, axis=-1):
    """
    Unwrap by changing deltas between values to 2*pi complement.

    Unwrap radian phase `p` by changing absolute jumps greater than
    `discont` to their 2*pi complement along the given axis.

    Parameters
    ----------
    p : array_like
        Input array.
    discont : float, optional
        Maximum discontinuity between values, default is ``pi``.
    axis : int, optional
        Axis along which unwrap will operate, default is the last axis.

    Returns
    -------
    out : ndarray
        Output array.

    See Also
    --------
    rad2deg, deg2rad

    Notes
    -----
    If the discontinuity in `p` is smaller than ``pi``, but larger than
    `discont`, no unwrapping is done because taking the 2*pi complement
    would only make the discontinuity larger.

    Examples
    --------
    >>> phase = np.linspace(0, np.pi, num=5)
    >>> phase[3:] += np.pi
    >>> phase
    array([ 0.        ,  0.78539816,  1.57079633,  5.49778714,  6.28318531])
    >>> np.unwrap(phase)
    array([ 0.        ,  0.78539816,  1.57079633, -0.78539816,  0.        ])

    """
    p = asarray(p)
    nd = len(p.shape)
    dd = diff(p, axis=axis)
    slice1 = [slice(None, None)]*nd     # full slices
    slice1[axis] = slice(1, None)
    ddmod = mod(dd + pi, 2*pi) - pi
    _nx.copyto(ddmod, pi, where=(ddmod == -pi) & (dd > 0))
    ph_correct = ddmod - dd
    _nx.copyto(ph_correct, 0, where=abs(dd) < discont)
    up = array(p, copy=True, dtype='d')
    up[slice1] = p[slice1] + ph_correct.cumsum(axis)
    return up