python类angle()的实例源码

def getAngleInDegrees(_fft):

    _angle = np.angle(_fft, True);
    _rows, _columns = _fft.shape;

    for i in range(0, _rows):
        for j in range(0, _columns):
            if (_angle[i][j] < 0 and abs(_angle[i][j]) <= 180):
                #print(_angle[i][j]);
                #print("\n");
                _angle[i][j] = 180 + _angle[i][j];
                #print(_angle[i][j]);
                #print("\n");
            elif (_angle[i][j] < 0 and abs(_angle[i][j]) > 180):
                #print(_angle[i][j]);
                #print("\n");
                _angle[i][j] = 360 + _angle[i][j];
                #print(_angle[i][j]);
                #print("\n");
    return _angle;

def appres(F, H, sig, chg, taux, c, mu, eps, n):

    Res = np.zeros_like(F)
    Phase = np.zeros_like(F)
    App_ImpZ= np.zeros_like(F, dtype='complex_')

    for i in range(0, len(F)):

        UD, EH, Z , K = Propagate(F[i], H, sig, chg, taux, c, mu, eps, n)

        App_ImpZ[i] = EH[0, 1]/EH[1, 1]

        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
        Phase[i] = np.angle(App_ImpZ[i], deg = True)

    return Res, Phase

# Evaluate Up, Down components, E and H field, for a frequency range,
# a discretized depth range and a time range (use to calculate envelope)

def phase_plot(data, toffset, log_scale, title):
    """Plot the phase of the data in linear or log scale."""
    print("phase")

    phase = numpy.angle(data) / numpy.pi

    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(phase)

    axmx = numpy.max(phase)

    #ax.axis([-axmx, axmx, -axmx, axmx])
    ax.grid(True)
    ax.set_xlabel('time')
    ax.set_ylabel('phase')
    ax.set_title(title)

    return fig

def unknown_feature_extractor(x, sr, win_len, shift_len, barks, inner_win, inner_shift, win_type, method_version):
    x_spectrum = stft_extractor(x, win_len, shift_len, win_type)
    coef = get_fft_bark_mat(sr, win_len, barks, 20, sr//2)
    bark_spect = np.matmul(coef, x_spectrum)
    ams = np.zeros((barks, inner_win//2+1, (bark_spect.shape[1] - inner_win)//inner_shift))
    for i in range(barks):
        channel_stft = stft_extractor(bark_spect[i, :], inner_win, inner_shift, 'hanning')
        if method_version == 'v1':
            ams[i, :, :] = 20 * np.log(np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift]))
        elif method_version == 'v2':
            channel_amplitude = np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
            channel_angle = np.angle(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
            channel_angle = channel_angle - (np.floor(channel_angle / (2.*np.pi)) * (2.*np.pi))
            ams[i, :, :] = np.power(channel_amplitude, 1./3.) * channel_angle
        else:
            ams[i, :, :] = np.abs(channel_stft)
    return ams

def ams_extractor(x, sr, win_len, shift_len, barks, inner_win, inner_shift, win_type, method_version):
    x_spectrum = stft_extractor(x, win_len, shift_len, win_type)
    coef = get_fft_bark_mat(sr, win_len, barks, 20, sr//2)
    bark_spect = np.matmul(coef, x_spectrum)
    ams = np.zeros((barks, inner_win//2+1, (bark_spect.shape[1] - inner_win)//inner_shift))
    for i in range(barks):
        channel_stft = stft_extractor(bark_spect[i, :], inner_win, inner_shift, 'hanning')
        if method_version == 'v1':
            ams[i, :, :] = 20 * np.log(np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift]))
        elif method_version == 'v2':
            channel_amplitude = np.abs(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
            channel_angle = np.angle(channel_stft[:inner_win//2+1, :(bark_spect.shape[1] - inner_win)//inner_shift])
            channel_angle = channel_angle - (np.floor(channel_angle / (2.*np.pi)) * (2.*np.pi))
            ams[i, :, :] = np.power(channel_amplitude, 1./3.) * channel_angle
        else:
            ams[i, :, :] = np.abs(channel_stft)
    return ams

def griffin_lim(mag, phase_angle, n_fft, hop, num_iters):
  """Iterative algorithm for phase retrival from a magnitude spectrogram.

  Args:
    mag: Magnitude spectrogram.
    phase_angle: Initial condition for phase.
    n_fft: Size of the FFT.
    hop: Stride of FFT. Defaults to n_fft/2.
    num_iters: Griffin-Lim iterations to perform.

  Returns:
    audio: 1-D array of float32 sound samples.
  """
  fft_config = dict(n_fft=n_fft, win_length=n_fft, hop_length=hop, center=True)
  ifft_config = dict(win_length=n_fft, hop_length=hop, center=True)
  complex_specgram = inv_magphase(mag, phase_angle)
  for i in range(num_iters):
    audio = librosa.istft(complex_specgram, **ifft_config)
    if i != num_iters - 1:
      complex_specgram = librosa.stft(audio, **fft_config)
      _, phase = librosa.magphase(complex_specgram)
      phase_angle = np.angle(phase)
      complex_specgram = inv_magphase(mag, phase_angle)
  return audio

def __init__(self, qubit_names, quad="real"):
        super(PulseCalibration, self).__init__()
        self.qubit_names = qubit_names if isinstance(qubit_names, list) else [qubit_names]
        self.qubit     = [QubitFactory(qubit_name) for qubit_name in qubit_names] if isinstance(qubit_names, list) else QubitFactory(qubit_names)
        self.filename   = 'None'
        self.exp        = None
        self.axis_descriptor = None
        self.cw_mode    = False
        self.saved_settings = config.load_meas_file(config.meas_file)
        self.settings = deepcopy(self.saved_settings) #make a copy for used during calibration
        self.quad = quad
        if quad == "real":
            self.quad_fun = np.real
        elif quad == "imag":
            self.quad_fun = np.imag
        elif quad == "amp":
            self.quad_fun = np.abs
        elif quad == "phase":
            self.quad_fun = np.angle
        else:
            raise ValueError('Quadrature to calibrate must be one of ("real", "imag", "amp", "phase").')
        self.plot       = self.init_plot()

def displayResult(self):
        fig = plt.figure(101)
        plt.subplot(221)
        plt.imshow(np.abs(self.reconstruction),origin='lower')
        plt.draw()
        plt.title('Reconstruction Magnitude')
        plt.subplot(222)
        plt.imshow(np.angle(self.reconstruction),origin='lower')
        plt.draw()
        plt.title('Reconstruction Phase')
        plt.subplot(223)
        plt.imshow(np.abs(self.aperture),origin='lower')
        plt.title("Aperture Magnitude")
        plt.draw()
        plt.subplot(224)
        plt.imshow(np.angle(self.aperture),origin='lower')
        plt.title("Aperture Phase")
        plt.draw()
        fig.canvas.draw()


        #fig.tight_layout()
        # display.display(fig)
        # display.clear_output(wait=True)
        # time.sleep(.00001)

def displayResult2(self):

        fig = plt.figure() 
        ax = fig.add_subplot(2,2,1 ) 
        ax.imshow(np.abs(self.reconstruction))         
        ax.set_title('Reconstruction Magnitude')

        ax = fig.add_subplot(2,2,2 ) 
        ax.imshow(np.angle(self.reconstruction))        
        ax.set_title('Reconstruction Phase')

        ax = fig.add_subplot(2,2,3 ) 
        ax.imshow(np.abs(self.aperture))
        ax.set_title("Aperture Magnitude")

        ax = fig.add_subplot(2,2,4 ) 
        ax.imshow(np.angle(self.aperture))
        ax.set_title("Aperture Phase")
        fig.tight_layout()

def phase_to_color_wheel(complex_number):
    """Map a phase of a complexnumber to a color in (r,g,b).

    complex_number is phase is first mapped to angle in the range
    [0, 2pi] and then to a color wheel with blue at zero phase.
    """
    angles = np.angle(complex_number)
    angle_round = int(((angles + 2 * np.pi) % (2 * np.pi))/np.pi*6)
    # print(angleround)
    color_map = {0: (0, 0, 1),  # blue,
                 1: (0.5, 0, 1),  # blue-violet
                 2: (1, 0, 1),  # violet
                 3: (1, 0, 0.5),  # red-violet,
                 4: (1, 0, 0),  # red
                 5: (1, 0.5, 0),  # red-oranage,
                 6: (1, 1, 0),  # orange
                 7: (0.5, 1, 0),  # orange-yellow
                 8: (0, 1, 0),  # yellow,
                 9: (0, 1, 0.5),  # yellow-green,
                 10: (0, 1, 1),  # green,
                 11: (0, 0.5, 1)  # green-blue,
                 }
    return color_map[angle_round]

def time_sync(self, signal, start_est, stop_est):
        self.symbol_found = 0
        self.symbol_start = 0
        self.log.debug('time_sync()')      
        corr_res = np.zeros( self.sym_len[self.mode], dtype=complex)
        for s in range(start_est, self.sym_len[self.mode]-1):
            corr_res[s] = self.gi_correlation(self.gi_len[self.mode], self.sym_len[self.mode], signal[s: s+self.sym_len[self.mode]] )
        corr_max = np.argmax( np.abs(corr_res) )     

        self.symbol_found = 1
        self.symbol_start = corr_max + 256 - 20
        self.generated_samples = self.sym_len[self.mode] - self.gi_len[self.mode]

        # Tracking
        self.fine_freq_off = np.angle(corr_res[corr_max])/(2*np.pi*(self.sym_len[self.mode] - self.gi_len[self.mode])*(1.0/self.fs))

        return

def compute_by_noise_pow(self, signal, n_pow):
        s_spec = np.fft.fftpack.fft(signal * self._window)
        s_amp = np.absolute(s_spec)
        s_phase = np.angle(s_spec)
        gamma = self._calc_aposteriori_snr(s_amp, n_pow)
        xi = self._calc_apriori_snr(gamma)
        self._prevGamma = gamma
        nu = gamma * xi / (1.0 + xi)
        self._G = (self._gamma15 * np.sqrt(nu) / gamma) * np.exp(-nu / 2.0) *\
                  ((1.0 + nu) * spc.i0(nu / 2.0) + nu * spc.i1(nu / 2.0))
        idx = np.less(s_amp ** 2.0, n_pow)
        self._G[idx] = self._constant
        idx = np.isnan(self._G) + np.isinf(self._G)
        self._G[idx] = xi[idx] / (xi[idx] + 1.0)
        idx = np.isnan(self._G) + np.isinf(self._G)
        self._G[idx] = self._constant
        self._G = np.maximum(self._G, 0.0)
        amp = self._G * s_amp
        amp = np.maximum(amp, 0.0)
        amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
        self._prevAmp = amp
        spec = amp2 * np.exp(s_phase * 1j)
        return np.real(np.fft.fftpack.ifft(spec))

def compute_by_noise_pow(self, signal, n_pow):
        s_spec = np.fft.fftpack.fft(signal * self._window)
        s_amp = np.absolute(s_spec)
        s_phase = np.angle(s_spec)
        gamma = self._calc_aposteriori_snr(s_amp, n_pow)
        xi = self._calc_apriori_snr(gamma)
        # xi = self._calc_apriori_snr2(gamma,n_pow)
        self._prevGamma = gamma
        nu = gamma * xi / (1.0 + xi)
        self._G = xi / (1.0 + xi) * np.exp(0.5 * spc.exp1(nu))
        idx = np.less(s_amp ** 2.0, n_pow)
        self._G[idx] = self._constant
        idx = np.isnan(self._G) + np.isinf(self._G)
        self._G[idx] = xi[idx] / (xi[idx] + 1.0)
        idx = np.isnan(self._G) + np.isinf(self._G)
        self._G[idx] = self._constant
        self._G = np.maximum(self._G, 0.0)
        amp = self._G * s_amp
        amp = np.maximum(amp, 0.0)
        amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
        self._prevAmp = amp
        spec = amp2 * np.exp(s_phase * 1j)
        return np.real(np.fft.fftpack.ifft(spec))

def compute_by_noise_pow(self, signal, n_pow):
        s_spec = np.fft.fftpack.fft(signal * self._window)
        s_amp = np.absolute(s_spec)
        s_phase = np.angle(s_spec)
        gamma = self._calc_aposteriori_snr(s_amp, n_pow)
        # xi = self._calc_apriori_snr2(gamma,n_pow)
        xi = self._calc_apriori_snr(gamma)
        self._prevGamma = gamma
        u = 0.5 - self._mu / (4.0 * np.sqrt(gamma * xi))
        self._G = u + np.sqrt(u ** 2.0 + self._tau / (gamma * 2.0))
        idx = np.less(s_amp ** 2.0, n_pow)
        self._G[idx] = self._constant
        idx = np.isnan(self._G) + np.isinf(self._G)
        self._G[idx] = xi[idx] / (xi[idx] + 1.0)
        idx = np.isnan(self._G) + np.isinf(self._G)
        self._G[idx] = self._constant
        self._G = np.maximum(self._G, 0.0)
        amp = self._G * s_amp
        amp = np.maximum(amp, 0.0)
        amp2 = self._ratio * amp + (1.0 - self._ratio) * s_amp
        self._prevAmp = amp
        spec = amp2 * np.exp(s_phase * 1j)
        return np.real(np.fft.fftpack.ifft(spec))

def plot_eigen_function(ax, se, n, psi1l, psi1r):

    # plt.figure(figsize=(8, 8))
    for x in range(se.info['L']):
        for y in range(se.info['W']):
            i = x + y * se.info['L']

            w = np.sqrt(10) * np.abs(psi1r[i, n])
            arg = np.angle(psi1r[i, n])

            circle = plt.Circle((x, y), w, color='black', zorder=10)
            ax.add_artist(circle)

            ax.plot([x, x + w * np.cos(arg)], [y, y + w * np.sin(arg)], color='white', lw=.8, zorder=12)

        ax.set_xlim([-.5, se.info['L'] - .5])
        ax.set_ylim([-.5, se.info['W'] - .5])
    # plt.show()

def single_spectrogram(inseq,fs,wlen,h,imag=False):
    """
        imag: Return Imaginary Data of the STFT on True 
    """
    NFFT = int(2**(np.ceil(np.log2(wlen)))) 
    K = np.sum(hamming(wlen, False))/wlen
    raw_data = inseq.astype('float32')
    raw_data = raw_data/np.amax(np.absolute(raw_data))
    stft_data,_,_ = STFT(raw_data,wlen,h,NFFT,fs)
    s = np.absolute(stft_data)/wlen/K;
    if np.fmod(NFFT,2):
        s[1:,:] *=2
    else:
        s[1:-2] *=2        
    real_data = np.transpose(20*np.log10(s + 10**-6)).astype(np.float32)
    if imag:
        imag_data = np.angle(stft_data).astype(np.float32)
        return real_data,imag_data 
    return real_data

def complex2rgbalin(s):
   """
   Displays complex image with intensity corresponding to the MODULUS and color (hsv) correponging to PHASE.
   From: pyVincent/ptycho.py
   """    
   ph=np.angle(s)
   t=np.pi/3
   nx,ny=s.shape
   rgba=np.zeros((nx,ny,4))
   rgba[:,:,0]=(ph<t)*(ph>-t) + (ph>t)*(ph<2*t)*(2*t-ph)/t + (ph>-2*t)*(ph<-t)*(ph+2*t)/t
   rgba[:,:,1]=(ph>t)         + (ph<-2*t)      *(-2*t-ph)/t+ (ph>0)*(ph<t)    *ph/t
   rgba[:,:,2]=(ph<-t)        + (ph>-t)*(ph<0) *(-ph)/t + (ph>2*t)         *(ph-2*t)/t
   a=np.abs(s)
   a/=a.max()
   rgba[:,:,3]=a
   return rgba

def autocorr(trace):
    """This function takes an obspy trace object and performs a phase autocorrelation
    of the trace with itself. First, the hilbert transform is taken to obtain the 
    analytic signal and hence the instantaneous phase. This is then passed to a 
    fortran script where phase correlation is performed after Schimmel et al., 2011.

    This function relies on the shared object file phasecorr.so, which is the file
    containing the fortran subroutine for phase correlation.
    """

    import numpy as np

    from scipy.signal import hilbert
    from phasecorr import phasecorr
    # Hilbert transform to obtain the analytic signal
    htrans = hilbert(trace)
    # Perform phase autocorrelation with instantaneous phase
    pac = phasecorr(np.angle(htrans),np.angle(htrans),len(htrans))

    return pac

def crosscorr(tr1,tr2):
    """This function takes 2 obspy traces and performs a phase crosscorrelation
    of the traces. First, the hilbert transform is taken to obtain the 
    analytic signal and hence the instantaneous phase. This is then passed to a 
    fortran script where phase correlation is performed after Schimmel et al., 2011.

    This function relies on the shared object file phasecorr.so, which is the file
    containing the fortran subroutine for phase correlation.
    """

    import numpy as np

    from scipy.signal import hilbert
    from phasecorr import phasecorr
    # Hilbert transform to obtain the analytic signal
    htrans1 = hilbert(tr1)
    htrans2 = hilbert(tr2)
    # Perform phase autocorrelation with instantaneous phase
    pcc = phasecorr(np.angle(htrans1),np.angle(htrans2),len(htrans1))

    return pcc

def __init__(self, gain):
        self.gain = gain * np.pi  # pi term puts angle output to pi range

        # components
        self.conjugate = Conjugate()
        self.complex_mult = ComplexMultiply()
        self.angle = Angle()
        self.out = Sfix()

        # specify component delay
        self._delay = self.conjugate._delay + \
                      self.complex_mult._delay + \
                      self.angle._delay + 1

        # constants
        self.gain_sfix = Const(Sfix(self.gain, 3, -14))

def NMF(stft, n_sources):
    """
    Sound source separation using NMF
    :param stft: the short-time Fourier Transform of the signal
    :param n_sources: the number of sources                                                                         
    :returns:                                                                                                         
      - Ys: sources
    """
    print "Computing approximations"
    W, H = librosa.decompose.decompose(np.abs(stft), n_components=n_sources, sort=True)
    print "Reconstructing signals"
    Ys = list(scipy.outer(W[:,i], H[i])*np.exp(1j*np.angle(stft)) for i in xrange(n_sources))
    print "Saving sound arrays"
    ys = [librosa.istft(Y) for Y in Ys]
    del Ys
    return ys

def undo_melspect(spect, sample_rate=SAMPLE_RATE, fps=FPS, frame_len=FRAME_LEN, min_freq=MIN_FREQ, max_freq=MAX_FREQ, invert_melbank_method='transpose', phases='random', normalize=False):
    """
    Resynthesizes a mel spectrogram into a numpy array of samples.
    """
    # undo logarithmic scaling
    spect = undo_logarithmize(spect)
    # undo Mel filterbank
    spect = undo_melfilter(spect, sample_rate, frame_len, min_freq, max_freq, invert_melbank_method)
    # randomize or reuse phases
    if phases == 'random':
        spect = spect * np.exp(np.pi*2.j*np.random.random(spect.shape))
    elif phases is not None:
        spect = spect * np.exp(1.j * np.angle(phases))
    # undo STFT
    hop_size = sample_rate / fps
    samples = undo_stft(spect, hop_size, frame_len)
    # normalize if needed
    if normalize:
        samples -= samples.mean()
        samples /= np.abs(samples).max()
    return samples.astype(np.float32)

def __new__(cls, realpart, imagpart=None):
        """Create a new EMArray."""

        # Create complex obj
        if np.any(imagpart):
            obj = np.real(realpart) + 1j*np.real(imagpart)
        else:
            obj = np.asarray(realpart, dtype=complex)

        # Ensure its at least a 1D-Array, view cls
        obj = np.atleast_1d(obj).view(cls)

        # Store amplitude
        obj.amp = np.abs(obj)

        # Calculate phase, unwrap it, transform to degrees
        obj.pha = np.rad2deg(np.unwrap(np.angle(obj.real + 1j*obj.imag)))

        return obj


# 2. Input parameter checks for modelling

# 2.a <Check>s (alphabetically)

def ph_dec(m_phs, m_phc, mode='angle'):  

    if mode == 'sign':    
        m_bs = np.arcsin(m_phs)
        m_bc = np.arccos(m_phc)   
        m_ph = np.sign(m_bs) * np.abs(m_bc)   

    elif mode == 'angle':
        m_ph = np.angle(m_phc + m_phs * 1j)

    return m_ph 


#==============================================================================
# From 'analysis_with_del_comp_and_ph_encoding_from_files'
# f0_type: 'f0', 'lf0'

def phase_string(sig):
    """Take the angle and create a string as \pi if close to some \pi values"""
    if isinstance(sig, np.ndarray):
        sig = sig.ravel()[0]
    if isinstance(sig, np.complex):
        angle = np.angle(sig)
    else:
        angle = sig

    fractions = [Fraction(i, 12) for i in np.arange(-12, 12 + 1)]
    pi_multiples = np.array([np.pi * frac_to_float(f) for f in fractions])
    strings = [frac_to_str(f) for f in fractions]

    if np.min(np.abs(angle - pi_multiples)) < 1e-3:
        return strings[np.argmin(np.abs(angle - pi_multiples))]
    else:
        return '%.2f' % angle

def cimshow(f_impulse):
    plt.figure(figsize=(7,5))
    plt.subplot(2,2,1)
    plt.imshow(np.real(f_impulse))
    plt.colorbar()
    plt.title('Re{}')

    plt.subplot(2,2,2)
    plt.imshow(np.imag(f_impulse))
    plt.colorbar()
    plt.title('Img{}')

    plt.subplot(2,2,2+1)
    plt.imshow(np.abs(f_impulse))
    plt.colorbar()
    plt.title('Magnitude')

    plt.subplot(2,2,2+2)
    plt.imshow(np.angle(f_impulse))
    plt.colorbar()
    plt.title('Phase')

def cimshow(f_impulse):
    plt.figure(figsize=(7,5))
    plt.subplot(2,2,1)
    plt.imshow(np.real(f_impulse))
    plt.colorbar()
    plt.title('Re{}')

    plt.subplot(2,2,2)
    plt.imshow(np.imag(f_impulse))
    plt.colorbar()
    plt.title('Img{}')

    plt.subplot(2,2,2+1)
    plt.imshow(np.abs(f_impulse))
    plt.colorbar()
    plt.title('Magnitude')

    plt.subplot(2,2,2+2)
    plt.imshow(np.angle(f_impulse))
    plt.colorbar()
    plt.title('Phase')

def cimshow(f_impulse):
    plt.figure(figsize=(7,5))
    plt.subplot(2,2,1)
    plt.imshow(np.real(f_impulse))
    plt.colorbar()
    plt.title('Re{}')

    plt.subplot(2,2,2)
    plt.imshow(np.imag(f_impulse))
    plt.colorbar()
    plt.title('Img{}')

    plt.subplot(2,2,2+1)
    plt.imshow(np.abs(f_impulse))
    plt.colorbar()
    plt.title('Magnitude')

    plt.subplot(2,2,2+2)
    plt.imshow(np.angle(f_impulse))
    plt.colorbar()
    plt.title('Phase')

def update_recon_py(Recon1, support):
    err1 = 1.0
    Constraint = np.ones(Recon1.shape)
    # cdef int R1y, R1x
    Recon1_abs = np.abs(Recon1)
    Recon1_pwr2 = np.power(Recon1_abs, 2)
    R1y, R1x = Recon1.shape

    for p in range(R1y):
        for q in range(R1x):
            if support[p, q] == 1:
                Constraint[p, q] = Recon1_abs[p, q]
                err1 += Recon1_pwr2[p, q]
            if Recon1_abs[p, q] > 1:
                Constraint[p, q] = 1

    Recon1_update = Constraint * np.exp(1j * np.angle(Recon1))

    return Recon1_update, err1

def imshow_GfpGbp(Gfp, Gbp):
    plt.subplot(2, 2, 1)
    plt.imshow(np.abs(Gfp), cmap='Greys_r')
    plt.title('abs(Gfp)')

    plt.subplot(2, 2, 2)
    plt.imshow(np.angle(Gfp), cmap='Greys_r')
    plt.title('angle(Gfp)')

    plt.subplot(2, 2, 1 + 2)
    plt.imshow(np.abs(Gbp), cmap='Greys_r')
    plt.title('abs(Gbp)')

    plt.subplot(2, 2, 2 + 2)
    plt.imshow(np.angle(Gbp), cmap='Greys_r')
    plt.title('angle(Gbp)')