Python numpy.unwrap() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

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 

Example 2

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) 

Example 3

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') 

Example 4

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) 

Example 5

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))) ] 

Example 6

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 = {} 

Example 7

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 

Example 8

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 

Example 9

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 

Example 10

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 

Example 11

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) 

Example 12

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) 

Example 13

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 

Example 14

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 

Example 15

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 

Example 16

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) 

Example 17

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) 

Example 18

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 

Example 19

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 

Example 20

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) 

Example 21

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() 

Example 22

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 

Example 23

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 

Example 24

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 

Example 25

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 

Example 26

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 

Example 27

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 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注