def test_random_odd(self):
     for n in [33,65,55]:
         f = random((n,))
         af = sum(f,axis=0)/n
         f = f-af
         assert_almost_equal(sum(f,axis=0),0.0)
         assert_array_almost_equal(ihilbert(hilbert(f)),f)
         assert_array_almost_equal(hilbert(ihilbert(f)),f)

 def test_definition(self):
     for n in [16,17,64,127]:
         x = arange(n)*2*pi/n
         y = hilbert(sin(x))
         y1 = direct_hilbert(sin(x))
         assert_array_almost_equal(y,y1)
         assert_array_almost_equal(hilbert(sin(2*x)),
                                   direct_hilbert(sin(2*x)))

    def envelope(self):
        '''
        Computes the envelope of the trace. 
        '''
        from scipy.fftpack import hilbert

        self.values = 20 * np.log10(np.abs(hilbert(self.values)))

def analytic_signal(x):
    """ A short-cut assuming that the incoming signal is reasonable
    e.g. fairly pure sinusoid.
    So far this has no strategy to minimize fourier transform time.
    """ 
    x = x - np.mean(x)
    return(x+1j*hilbert(x))

 def bench_random(self):
     print()
     print(' Hilbert transform of periodic functions')
     print('=========================================')
     print(' size  | optimized |    naive')
     print('-----------------------------------------')
     for size,repeat in [(100,1500),(1000,300),
                         (256,1500),
                         (512,1000),
                         (1024,500),
                         (2048,200),
                         (2048*2,100),
                         (2048*4,50),
                         ]:
         print('%6s' % size, end=' ')
         sys.stdout.flush()
         x = arange(size)*2*pi/size
         if size < 2000:
             f = sin(x)*cos(4*x)+exp(sin(3*x))
         else:
             f = sin(x)*cos(4*x)
         assert_array_almost_equal(hilbert(f),direct_hilbert(f))
         print('| %9.2f' % measure('hilbert(f)',repeat), end=' ')
         sys.stdout.flush()
         print('| %9.2f' % measure('direct_hilbert(f)',repeat), end=' ')
         sys.stdout.flush()
         print(' (secs for %s calls)' % (repeat))

    def get_iprobe(self, leakage=None, t_comp=None):
        """ main purpose is to subtract leakage currents
        Will use the full data, as the t_range is meant to be plasma interval
        returns a copy of the measured courremt, overwritten with the
        corrected iprobe
        """
        # obtain leakage estimate
        if t_comp is None:
            t_comp = self.t_comp
        FFT_size = nice_FFT_size(len(self.imeasfull.timebase), -1)
        self.iprobefull = self.imeasfull.copy()
        self.sweepQ = []  # will keep these for synchronous sampling
        input_leakage = leakage
        for (c, chan) in enumerate(self.imeasfull.channels):
            if self.select is not None and c not in self.select:
                continue
            leakage = input_leakage
            cname = chan.config_name
            sweepV = self.vcorrfull.signal[self.vlookup[self.vassoc[c]]][0:FFT_size]
            sweepQ = hilbert(sweepV)
            self.sweepQ.append(sweepQ)  # save for synchronising segments (it is smoothed)

            # these attempts to make it accept a single channel are only partial
            imeas = self.imeasfull.signal[c] # len(self.imeasfull.channels) >1 else self.imeasfull.signal
            tb = self.imeasfull.timebase

            w_comp = np.where((tb>=t_comp[0]) & (tb<=t_comp[1]))[0]
            if len(w_comp) < 2000:

                raise ValueError('Not enough points {wc} t_comp - try {tt}'
                    .format(tt=np.round([tb[0], tb[0] + t_comp[1]-t_comp[0]],3),
                            wc=len(w_comp)))
            ns = len(w_comp)
            wind = np.blackman(ns)
            offset = np.mean(wind * imeas[w_comp])/np.mean(wind)
            sweepVFT = np.fft.fft(AC(sweepV[w_comp]) * wind)
            imeasFT = np.fft.fft(AC(imeas[w_comp]) * wind)
            ipk = np.argmax(np.abs(sweepVFT)[0:ns//2])  # avoid the upper one
            comp = imeasFT[ipk]/sweepVFT[ipk]

            #print('leakage compensation factor = {r:.2e} + j{i:.2e}'
            #      .format(r=np.real(comp), i=np.imag(comp)))
            print('{u}sing computed leakage comp factor = {m:.2e} e^{p:.2f}j'
                  .format(u = ["Not u", "U"][leakage is None],
                          m=np.abs(comp), p=np.angle(comp)))
            if leakage is None:
                leakage = [np.real(comp), np.imag(comp)]

            # find the common length - assuming they start at the same time????
            comlen = min(len(self.imeasfull.timebase),len(self.vmeasfull.timebase),len(sweepQ))
            # put signals back into rdata (original was copied by reduce_time)
            # overwrite - is this OK?
            self.iprobefull.signal[c] = self.iprobefull.signal[c]*0.  # clear it
            # sweepV has a DC component! beware
            self.iprobefull.signal[c][0:comlen] = self.imeasfull.signal[c][0:comlen]-offset \
                                        - sweepV[0:comlen] * leakage[0] - sweepQ[0:comlen] * leakage[1]
            # remove DC cpt (including that from the compensation sweepV)
            offset = np.mean(wind * self.iprobefull.signal[c][w_comp])/np.mean(wind)
            self.iprobefull.signal[c][0:comlen] -= offset

 def test_tilbert_relation(self):
     for n in [16,17,64,127]:
         x = arange(n)*2*pi/n
         f = sin(x)+cos(2*x)*sin(x)
         y = hilbert(f)
         y1 = direct_hilbert(f)
         assert_array_almost_equal(y,y1)
         y2 = tilbert(f,h=10)
         assert_array_almost_equal(y,y2)

def frequency_estimate_ml(x):
    x_ = scifft.hilbert(x)
    x_ = noiseWhite(x_)
    datalen = len(x_)

    for i in range(0, datalen):
        x_[i] = np.math.log(abs(x_[i]))/(i+1)
    averge_ = np.average(x_)
    result_ = abs(averge_)
    return result_

 def test_random_even(self):
     for n in [32,64,56]:
         f = random((n,))
         af = sum(f,axis=0)/n
         f = f-af
         # zeroing Nyquist mode:
         f = diff(diff(f,1),-1)
         assert_almost_equal(sum(f,axis=0),0.0)
         assert_array_almost_equal(direct_hilbert(direct_ihilbert(f)),f)
         assert_array_almost_equal(hilbert(ihilbert(f)),f)

    def envelope(self):
        '''
        NEEDS TESTING
        
        Computes the envelope of the traces using the Hilbert transform. 
        '''
        from scipy.fftpack import hilbert

        self.info('Applying Hilbert transform ...')
        for trace in range(self.traces):
            self.data[:, trace] = np.abs(hilbert(self.data[:, trace]))
        self.done()

def compressed_sensing_1VD( fid, mask, num_iter=500, factor=0.95, tol = 0.01, maxPeaks=2 ):
    
    sss = numpy.zeros( len(fid))
    sss0 = numpy.zeros( len(fid))
    
    final_iter_value  = 0
    final_tol_value   = 0
    final_numpeaks_value = 0

    fid1 = fid.copy()
    tol_diff = (numpy.sqrt(((abs(fid1)).sum())/32.))
    
    k=0
    kkk = []
    rrr = fftpack.fft(fid1)
    rss0 = []
    rss0.append(abs(rrr).sum())
    
    tol0 = abs(rrr).sum()
    tol_diff = ( tol0 - abs(rrr).sum() )*100.0 / tol0
    while (tol_diff < tol) and (k < num_iter) and numPeaks( sss ) <maxPeaks:
        
        sss0 = 1.0*sss
        
        rrr = fftpack.fft(fid1)
        m1 = max(rrr.real)
        
        sss_max_old = sss.max()
               
        for i,r in enumerate(rrr):
            if r.real > m1* factor:
                sss[i] = sss[i]+rrr[i].real
                rrr[i] = complex(m1*factor)
        sss_max = sss.max()
        
        rrr_iii = fftpack.hilbert( rrr.real )
        
        rrr = rrr.real + 1j * rrr_iii
        
        fid1 = fftpack.ifft(rrr)
        
        fid1 *= mask
        tol_diff = ( tol0 - abs(rrr).sum() )*100.0 / tol0
        k +=1

    final_iter_value = k
    final_numpeaks_value = numPeaks(sss)
    final_tol_value = tol_diff
    
    return( sss0, [final_iter_value, final_numpeaks_value, final_tol_value ] )

 def reconstruct(self, paData):
     if paData.ndim == 2:
         (nSamples, nSteps) = paData.shape
         paData = np.reshape(paData, (nSamples, nSteps, 1))
     (nSamples, nSteps, zSteps) = paData.shape
     reImg1 = np.copy(super().reconstruct(paData))
     # take 90-degree phase shift
     import scipy.fftpack as spfp
     for z in range(zSteps):
         for n in range(nSteps):
             paData[:, n, z] = spfp.hilbert(paData[:, n, z])
     reImg2 = np.copy(super().reconstruct(paData))
     self.reImg = np.sqrt(reImg1 ** 2 + reImg2 ** 2)
     return self.reImg

def analytic_phase(x, t=None, subint=None):
    """ gets the phase from an amazing variety of signals
    http://en.wikipedia.org/wiki/Analytic_signal
    subinterval idea is not debugged and is probably unnecessary
    may shorten data?
    """
    from scipy.fftpack import hilbert
    from pyfusion.utils import fix2pi_skips
    from numpy import zeros, arctan2
    # this subinterval idea does not seem necessary and is not debugged
    if subint != None:
        subint=powerof2(subint)
        nsubints = int(len(x)/subint)
        phs = zeros([nsubints, subint])
        for i in range(nsubints):
            xsub = x[i*subint:(i+1)*subint]
            phs[i,:] = arctan2(xsub, hilbert(xsub))
        phi = phs.flatten()
    else:
        y=hilbert(x)  # use hilbert twice just to remove DC (lazy)
        phi = arctan2(y, hilbert(y))

    return(fix2pi_skips(phi, sign='+'))

def envelope(data):
    """
    Envelope of a function.

    Computes the envelope of the given function. The envelope is determined by
    adding the squared amplitudes of the function and it's Hilbert-Transform
    and then taking the square-root. (See [Kanasewich1981]_)
    The envelope at the start/end should not be taken too seriously.

    :param data: Data to make envelope of, type numpy.ndarray.
    :return: Envelope of input data.
    """
    hilb = hilbert(data)
    data = (data ** 2 + hilb ** 2) ** 0.5
    return data

def Kosambi_Hilbert_phase(X, sampling_rate, passband=None, index=0, moving_window=False):

	y = Kosambi_Hilbert_torsion(X, sampling_rate, passband=passband, index=index, moving_window=moving_window)
	Hy = hilbert(y)

	radius = np.sqrt(y**2+Hy**2)
	phase = np.arctan2(y, Hy)
	phi_u = unmod(phase)

	if phi_u[-1]-phi_u[0] < 0.:	# if phase shrinks, reverse it.
		phase = -phase

	phase = np.mod(phase, pi2)

	return phase, radius

def envelope(x):
    """Computes the envelope of a seismic signal.

    The envelope, e(n), of a signal x(n), is calculated as:

        e(n) = (x(n) ** 2 + h(n) ** 2) ** 0.5

    where h(n) is the Hilbert Transform of x(n)

    Args:
        x: array of data.

    Returns:
        out: Envelope of x, numpy array type.
    """
    x_mean = x.mean()
    x_norm = x - x_mean
    return ((x_norm ** 2 + fftpack.hilbert(x_norm) ** 2) ** 0.5) + x_mean

def _Kosambi_Hilbert_torsion(X, Filter, index=0):
	# X[time, channel] should be X_j(t_i)

	if X.shape[1] == 1: return X[:, 0]

	# declarations and initializations
	X = np.asarray(X, dtype=float)
	channels = range(X.shape[1])
	Y = np.zeros((X.shape[0], 2*X.shape[1]+1), float)			# Y[time,  channel], X, and H(X) with Y[:, 0] as the reference channel.
	Yf = np.zeros(Y.shape, float)						# Filter(ed) version of Y.


	X, reference_amplitude = normalize(X, Filter, index=index)
	Y[:, 0] = X[:, index]							# save the reference channel to vector zero
	channels.pop(index)							# reference channel is treated separately.


	for (c, channel) in enumerate(channels):
		Y[:, 1+2*c] = X[:, channel]
		Y[:, 1+2*c+1] = hilbert(Y[:, channel])


	for i in xrange(Y.shape[1]):
		Yf[:, i] = Filter(Y[:, i])


	pcanode = mdp.nodes.PCANode(svd=True)	# this pcanode is used by the function below. (it's actually some static variable)
	pcanode.execute(Yf)			# get the principle components from Yf
	Proj = pcanode.get_projmatrix()		# ...and their projection matrix.

	if Proj[0, 0] < 0: Proj = -Proj		# ... why do I need to do this?


	KHT_component = np.dot(Y, Proj)[:, 0]		# apply them to Y!!!
	pca_amplitude = np.sqrt(signalAndNoise(KHT_component, Filter)[0])

	return reference_amplitude/pca_amplitude * KHT_component 

# -*- coding: utf-8 -*-
from scipy import fftpack
import numpy as np

x = np.random.rand(16)
y = fftpack.hilbert(x)

X = np.fft.fft(x)
Y = np.fft.fft(y)

print np.imag(Y/X)

def getPhase(x):     
    return numpy.arctan2(x,hilbert(x))

# -*- coding: utf-8 -*-
from scipy import signal
from scipy import fftpack
import numpy as np
import pylab as pl

# 某个均衡滤波器的参数
a = np.array([1.0, -1.947463016918843, 0.9555873701383931])
b = np.array([0.9833716591860479, -1.947463016918843, 0.9722157109523452])

# 44.1kHz， 1秒的频率扫描波
t = np.arange(0, 0.5, 1/44100.0)
x= signal.chirp(t, f0=10, t1 = 0.5, f1=1000.0)

# 直接一次计算滤波器的输出
y = signal.lfilter(b, a, x)

hy = fftpack.hilbert(y)
pl.plot( np.sqrt(y**2 + hy**2),"r", linewidth=2) 
pl.plot(y)
pl.show()