Python numpy.kaiser函数代码示例

def GenerateBLI(sampling_rate = 48000.0,
                cutoff = 15000.0,
                length = 4,
                ppiv = 2700,
                beta = 8.3,
                apodization_factor = 0.5,
                apodization_beta = 0.5):
    """
    Generates a bandlimited impulse
    """
    kEpsilon = 1e-7
    points_count = ppiv * length
    xaxis = numpy.linspace(0.0, points_count, points_count)
    x2 = length * 2.0 * (xaxis - (points_count) / 2.0 + kEpsilon) / (points_count)
    x3 = numpy.pi * cutoff / sampling_rate * x2

    brickwall_impulse = numpy.sin(x3) / x3

    kaiser_window = numpy.kaiser(points_count, beta)
    bli_data = numpy.transpose(brickwall_impulse) * kaiser_window

    # Apodization
    apodization_window = (1.0 - apodization_factor) \
                            * numpy.kaiser(points_count,apodization_beta)
    bli_data *= apodization_window

    return bli_data

def plot_fourier(t, y, subplotnums=[[2,1,1],[2,1,2]], savedata="none", kaiserwindowparameter=0):
   ''' Function for plotting fourier series of some given arrays t and y
       :param t         time variable (plotted along the x-axis)
       :param y         variable to be plotted along the y-axis
       :subplotnums     subplotnum for plotting multiple plots in one window
   '''
   if len(subplotnums) - 1 != len(np.atleast_1d(kaiserwindowparameter)):
      print "Bad subplotnum vs kaiserwindowparameter"
      return
   # First check the t array whether it has a constant dt
   dt = t[1] - t[0]
   for i in xrange(len(t)-1):
      if dt != t[i+1] - t[i]:
         print "Gave bad timestep to plot_fourier, the time step in array t must be constant (for now)"
   # Use kaiser window on y
   _y = y*np.kaiser(len(y), np.atleast_1d(kaiserwindowparameter)[0])
   # Plot the raw data as well as the fitted data (fourier series)
   pl.subplot(subplotnums[0][0],subplotnums[0][1],subplotnums[0][2])
   pl.plot(t,_y,'.',color='r')
   # Do FFT on the data
   fourier=np.fft.fft(_y)*(1/(float)(len(t)))
   # Get frequencies of the fourier
   freq=np.fft.fftfreq(len(fourier), d=dt)
   # Declare t2 (Note: This is the same as t but we want the steps to be thicker so the data looks smoother
   dt2=dt*0.01
   t2=np.arange(len(t)*100)*dt2
   # Declare y2
   y2=np.array([np.sum(fourier*np.exp(complex(0,1)*2*np.pi*freq*T)) for T in t2])
   pl.plot(t2,y2,'-',color='b')
   pl.legend(["data", "fourier_fit"])
   # Plot the frequency spectrums
   j = 0
   for i in subplotnums[1:]:
      pl.subplot(i[0],i[1],i[2])
      # Use window function on the data
      _y = y*np.kaiser(len(y), np.atleast_1d(kaiserwindowparameter)[j])
      # Do FFT on the data
      fourier=np.fft.fft(_y)*(1/(float)(len(t)))
      # Get frequencies of the fourier
      freq=np.fft.fftfreq(len(fourier), d=dt)
      # Get the indexes:
      toIndex = (int)((len(freq)/2)/2.0 + 1)
      #maxValue = np.max(np.abs(fourier[1:len(fourier)/2]))
      #while True:
      #   if toIndex <= 2 or np.abs(fourier[toIndex-1]) > 0.02*maxValue:
      #      break
      #   toIndex = toIndex - 1
      #pl.plot(freq[1:len(freq)/2],2*np.abs(fourier[1:len(fourier)/2]), marker='.', linestyle='-', linewidth=0.5)
      if savedata != "none":
         saveplotdata( freq[1:toIndex], np.log(2*np.abs(fourier[1:toIndex])), savedata + "_" + str(j) + ".npy" )
      pl.plot(freq[1:toIndex],2*np.abs(fourier[1:toIndex]), marker='.', linestyle='-', linewidth=0.5)
      pl.legend(["frequency_spectrum"])
      pl.ylim([0,1.05*max(2*np.abs(fourier[1:len(fourier)/2]))])
      #pl.xlim([0, 0.13])
      #xTicks = np.arange(14)/100.0
      #pl.xticks(xTicks)
      # Put interactive mode on and show the plot:
      #pl.ion()
      #pl.show()
      j = j + 1

 def make_sense(self,u0):
     st=self.st
     L=numpy.shape(u0)[-1]
     u0dims= numpy.ndim(u0)
     
     rows=numpy.shape(u0)[0]
     cols=numpy.shape(u0)[1]
     # dpss rely on ctypes which are not suitable for cx_freeze 
     #    [dpss_rows, eigens] = spectrum.mtm.dpss(rows, 50, 1)
     #    [dpss_cols, eigens] = spectrum.mtm.dpss(cols, 50, 1)
     dpss_rows = numpy.kaiser(rows, 100)
     dpss_cols = numpy.kaiser(cols, 100)
     
     dpss_rows = numpy.fft.fftshift(dpss_rows)
     dpss_cols = numpy.fft.fftshift(dpss_cols)
     
     dpss_rows = appendmat(dpss_rows,cols)
     dpss_cols  = appendmat(dpss_cols,rows)
     
     
     dpss_fil=dpss_rows*dpss_cols.T # low pass filter
     
     rms=numpy.sqrt(numpy.mean(u0*u0.conj(),-1)) # Root of sum square
     st['sensemap']=numpy.ones(numpy.shape(u0),dtype=numpy.complex64)
     
     #    print('L',L)
     #    print('rms',numpy.shape(rms))
     for ll in numpy.arange(0,L):
         st['sensemap'][:,:,ll]=(u0[:,:,ll]+1e-15)/(rms+1e-15)
         st['sensemap'][:,:,ll]=scipy.fftpack.ifft2(scipy.fftpack.fft2(st['sensemap'][:,:,ll])*dpss_fil)
     return st

    def _kaiser_window(self, width, beta):
        """
        Generates a kaiser window

        beta Window shape
        0    Rectangular
        5    Similar to a Hamming
        6    Similar to a Hann
        8.6  Similar to a Blackman
        """
        return np.kaiser(width[0], beta).reshape((-1,1,1)) * np.kaiser(width[1], beta).reshape((-1,1)) * np.kaiser(width[2], beta)

 def signal_filter(self, channel, freqs=(1, 15)):
     """Filter signal from specific channel. A FFT is performed in
     the signal, and the result is multiplied by a window function
     (kaiser function), which nulls the undesired beating
     frequencies, that are outside of 'freqs'. The signal is then
     recreated with an IFFT."""
     # zero padding size:
     zer_pad_filter = 4
     # alias for sweep_size
     N = int(self.sweep_size)
     # FFT with zero padding
     fft = np.fft.rfft(self.single_sweep_data[channel], zer_pad_filter * N)
     # bp=fft[:] #used for other plot functions
     fmin, fmax = freqs
     # creates the beating frequency axis, and finds the position of
     # the frequency limits in the axis
     fft_freq = np.linspace(0, self.rate  / 2.,
                            num=zer_pad_filter * N / 2)
     cmin = (abs(fft_freq - fmin)).argmin()
     cmax = (abs(fft_freq - fmax)).argmin()
     # creates window function for filter. Second argument of kaiser
     # function must be float
     window_func = np.concatenate((np.zeros(cmin + 1),
                                   np.kaiser(cmax - cmin, 2.), np.zeros(zer_pad_filter * N / 2 - cmax)))
     # multiply the window by the signal's FFT.
     bp = np.multiply(fft, window_func)
     # Since the signal is REAL, we use IRFFT, and takes only the
     # first half, since the other half is symmetric.
     newsig = np.fft.irfft(bp)[:N]
     return newsig

def main():
    import sys
    from scipy import misc
    softness = 2
    r_scale = 3
    windowmaker = lambda x: np.kaiser(x, softness)
    for infilename in sys.argv[1:]:
        print("File: %s" % infilename)
        array = misc.imread(infilename)
        circles = get_circles(array, 5)
        l_pnoise, l_order, c_order, r_order, r_pnoise = circles

        window = get_holed_window(windowmaker, l_order[2] * r_scale, 10)
        mask = get_mask(array.shape, window, l_order[1])
        window = windowmaker(l_pnoise[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, l_pnoise[1])
        window = windowmaker(c_order[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, c_order[1])
#        showimage(mask * 255)
#        showimage(logscale(mask * array))

        window = get_holed_window(windowmaker, r_order[2] * r_scale, 10)
        mask = get_mask(array.shape, window, r_order[1])
        window = windowmaker(r_pnoise[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, r_pnoise[1])
        window = windowmaker(c_order[2] * r_scale)
        mask *= 1 - get_mask(array.shape, window, c_order[1])

def periodogram(x,y,z,t,name):
    vel=x*x+y*y
    components=[vel,z*z]
    F=[]
    A=[]

    for i in range(0,len(components)):
        window=np.kaiser(components[i].shape[-1],5)
        vel=components[i]*window
        f,a=DO_FFT(vel,20)
        F.append(f)
        A.append(a)

    
    plt.figure(1,figsize=(11,7))
    plt.title("Spectral Gap - LiCor Fontanella1 Precampagna")
    plt.subplot(121)
    plt.title("FFT horizontal velocity")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Power Spectrum [dB]")
    plt.loglog(F[0],A[0],'k', markersize=0.1)

    plt.subplot(122)
    plt.title("FFT vertical velocity")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Power Spectrum [dB]")
    plt.loglog(F[1],A[1],'r',markersize=0.1)

    try:
        plt.savefig("graph/"+name+"_FFT.png", figuresize=(8,6), dpi=320, format="png")
        print_ok("Graph saved in: "+"graph/"+name+"_FFT.png")
    except IOError as IoE:
        print_fail("I/O Error! Erro number = {0}; {1}".format(IoE.errno,IoE.strerror))
        exit(IoE.errno)

def kaiser(active):
    """Kaiser tapering window."""
    idx = _windowidx(active)
    # compute coefficients
    window = np.zeros(active.shape)
    window[idx] = np.kaiser(len(idx), 2)
    return window

def hps(x, fs=44100, lf=255, harmonics=3, precision=2, window=lambda l:_np.kaiser(l, 7.14285)):
    """ Estimates the pitch (fundamental frequency) of the given sample array by a standard HPS implementation. """
    x -= _np.mean(x)
    N = x.size
    w = x*window(N)

    # Append zeros to the end of the window so that each bin has at least the desired precision.
    if fs/N > precision:
        delta = int(fs/precision) - N
        w = _np.append(w, _np.zeros(delta))
        N = w.size

    X = _np.log(_np.abs(_np.fft.rfft(w)))

    # Sequentially decimate 'X' 'harmonics' times and add it to itself.
    # 'precision < fs/N' must hold, lest the decimation loses all the precision we'd gain.
    hps = _np.copy(X)
    for h in range(2, 2 + harmonics):
        dec = _sig.decimate(X, h)
        hps[:dec.size] += dec*(0.8**h)

    # Find the bin corresponding to the lowest detectable frequency.
    lb = lf*N/fs

    # And then the bin with the highest spectral content.
    arg_peak = lb + _np.argmax(hps[lb:dec.size])

    # TODO: Return the full array? A ranked list of identified notes?
    return fs*arg_peak/N

def single_taper_spectrum(data, delta, taper_name=None):
    """
    Returns the spectrum and the corresponding frequencies for data with the
    given taper.
    """
    length = len(data)
    good_length = length // 2 + 1
    # Create the frequencies.
    # XXX: This might be some kind of hack
    freq = abs(np.fft.fftfreq(length, delta)[:good_length])
    # Create the tapers.
    if taper_name == 'bartlett':
        taper = np.bartlett(length)
    elif taper_name == 'blackman':
        taper = np.blackman(length)
    elif taper_name == 'boxcar':
        taper = np.ones(length)
    elif taper_name == 'hamming':
        taper = np.hamming(length)
    elif taper_name == 'hanning':
        taper = np.hanning(length)
    elif 'kaiser' in taper_name:
        taper = np.kaiser(length, beta=14)
    # Should never happen.
    else:
        msg = 'Something went wrong.'
        raise Exception(msg)
    # Detrend the data.
    data = detrend(data)
    # Apply the taper.
    data *= taper
    spec = abs(np.fft.rfft(data)) ** 2
    return spec, freq

def preprocess(X, y=None, box_width=256, overlap=32, pad_width=0):
    X = pre_flow(X)
    
    out_X = []
    out_y = []
    m = int(X.shape[0] / box_width) * box_width
    
    w = np.kaiser(box_width + pad_width, 14)

    for start in range(0, m, overlap):
        end = start + box_width
        if end >= len(X):
            break
        
        x = X[start:end,]

        if pad_width > 0:
            x = np.vstack([x, np.zeros((pad_width, x.shape[1]))])

        x = (x.T * w).T

        psd = abs(np.fft.rfft(x, axis=0))

        #freqs, psd = signal.welch(X[start:end,], axis=0)
    
        out_X.append(psd.flatten())
        if y != None:
            out_y.append(stats.mode(y[start:end])[0][0])

    if y != None:
        return np.array(out_X), np.array(out_y)
    else:
        return np.array(out_X)

def simple_segmentation( wav_data, sample_rate, db_diff = 17, max_silence = 0.8 ): 

    # filter data
    filtered_data = filter_chain(wav_data, sample_rate, 1000, 0.6)

    # create figure and axis
    f, ax = plt.subplots()

    # get spectograma nd spectogram plot (i.e. STFT)
    NFFT = 1024
    noverlap = 256
    window = np.kaiser(NFFT, 8)
    spectogram, freqs, time_array, im = ax.specgram(filtered_data, NFFT=NFFT,
                                                    Fs=sample_rate, noverlap=noverlap,
                                                    window=window, cmap = 'Greys')
    ax.set_xlim(0, len(wav_data)/float(sample_rate)) # remove empty plot region
    ax.set_ylim(0, sample_rate/2.)                   # remove empty plot region

    syllables = find_syllables(spectogram, time_array, db_diff)

    segments = join_in_segments(syllables, max_silence) 

    # add segments to plot as red shadowed areas
    hspans = []
    for segment in segments:
        hspans.append(ax.axvspan(segment[0], segment[1], facecolor='r', alpha = 0.2))
    
    return segments, f, ax

def plot_Temperature(fig_size=None):
	#{{{
	N = 256
	hammingWindow = np.hamming(N)
	hanningWindow = np.hanning(N)
	kaiserWindow = np.kaiser(N,5)
	fig = plt.figure()
	if fig_size != None:
		fig.set_size_inches(fig_size[0], fig_size[1])
	for i in range(0,10):
		Time = Array[i-1][:,0]
		PlData1 = Array[i-1][:, 5]
		PlData2 = Array[i-1][:, 6]
		PlData3 = Array[i-1][:, 7]
		PlData4 = Array[i-1][:, 8]
		PlData5 = Array[i-1][:, 9]
		PlData6 = Array[i-1][:,10]
		PlData7 = Array[i-1][:,16]
		ax = fig.add_subplot(5,2,i)
		#plt.legend((NameIndex[5],NameIndex[6],NameIndex[7],NameIndex[8],NameIndex[9],NameIndex[10]),'upper left')
		#plt.plot(Time,PlData1,label=' 5-OriTemp')
		#plt.plot(Time,PlData2,label=' 6-InjTemp')
		#plt.plot(Time,PlData3,label=' 7-NUWTemp')
		plt.plot(Time,PlData4,label=' 8-NUCTemp')
		plt.plot(Time,PlData5,label=' 9-NDCTemp')
		plt.plot(Time,PlData6,label='10-FrnTemp')
		plt.plot(Time,PlData7,label='16-MixBoil')
		plt.legend(fontsize=15,loc = 'upper right')
		plt.xlabel('Time[s]')
		plt.ylabel('Temperature[K]')
		plt.xlim(0,10)
		plt.ylim(-200,1400)

def kaiser(x, width=None, weights=None, do_fit_edges=False):
    """Smooth the values in `x` with the Kaiser windowed filter.

    See: https://en.wikipedia.org/wiki/Kaiser_window

    Parameters
    ----------
    x : array-like
        1-dimensional numeric data set.
    width : float
        Fraction of x's total length to include in the rolling window (i.e. the
        proportional window width), or the integer size of the window.
    """
    if len(x) < 2:
        return x
    if width is None:
        width = guess_window_size(x, weights)
    x, wing, signal = check_inputs(x, width, False)
    # Apply signal smoothing
    window = np.kaiser(2 * wing + 1, 14)
    if weights is not None:
        y, _w = convolve_weighted(window, signal, weights)
    else:
        y = convolve_unweighted(window, signal, wing)
    if do_fit_edges:
        _fit_edges(x, y, wing)  # In-place
    return y

def gen_map(vals, n=nbins, hue=False):
	# saturation
	fvals = vals.flatten()
	yh, xh, patches = plt.hist(fvals, bins=n, range=(0,1), normed=False , cumulative=False , histtype='step')
	if hue:
		# apply window
		M = 9
		win = np.kaiser(M, 3.0)
		yh = np.insert(yh, 0, np.zeros(M/2))
		yh = np.append(yh, np.zeros(M/2))
		yh = rolling_window(yh.T, M)
		yh = np.dot(yh, win)
	yh /= sum(yh)
	if hue:
		# adapted norm
		#yh = np.minimum(yh, hcut)
		yh[yh<=hcut] = 0
		yh /= sum(yh)
	yh = np.cumsum(yh)

	xhi = np.linspace(0,1,256)
	yhi = np.interp(xhi, yh, xh[1:])
	yhinv = np.interp(xhi, xh[1:], yh)
	#plt.plot(xhi, yhi)
	return (yhi, yhinv)