Python numpy.blackman函数代码示例

def noiseFilter(data_in):

	N = int(np.ceil((4 / b)))
	if not N % 2: N += 1  # Make sure that N is odd.
	n = np.arange(N)
	 
	# Compute a low-pass filter with cutoff frequency fH.
	hlpf = np.sinc(2 * fH * (n - (N - 1) / 2.))
	hlpf *= np.blackman(N)
	hlpf = hlpf / np.sum(hlpf)
	 
	# Compute a high-pass filter with cutoff frequency fL.
	hhpf = np.sinc(2 * fL * (n - (N - 1) / 2.))
	hhpf *= np.blackman(N)
	hhpf = hhpf / np.sum(hhpf)
	hhpf = -hhpf
	hhpf[(N - 1) / 2] += 1
	 
	# Convolve both filters.
	h = np.convolve(hlpf, hhpf)
	s = np.convolve(data_in, hlpf)

	fig, ax = plt.subplots()
	ax.plot(data_in)
	plt.show()
	fig1, ax1 = plt.subplots()
	ax1.plot(s)
	plt.show()
	return s

 def applyWindow(self, window="hanning", ww=0, cf=0):
     '''
     Apply window function to frequency domain data
     cf: the frequency the window is centered over [Hz]
     ww: the window width [Hz], if ww equals 0 the window covers the full range
     '''
     self.info("Applying %s window ..." % window)
     if window == "hanning":
         if ww == 0:
             w = np.hanning(self.numfreq)
         else:
             pos = int((cf - self.lowF) / self.deltaF)
             halfwidth = int(ww / (2.0 * self.deltaF))
             w = np.zeros(self.numfreq)
             w[pos - halfwidth:pos + halfwidth] = np.hanning(2 * halfwidth)
     elif window == "hamming":
         if ww == 0:
             w = np.hamming(self.numfreq)
         else:
             pos = int((cf - self.lowF) / self.deltaF)
             halfwidth = int(ww / (2.0 * self.deltaF))
             w = np.zeros(self.numfreq)
             w[pos - halfwidth:pos + halfwidth] = np.hamming(2 * halfwidth)
     elif window == "blackman":
         if ww == 0:
             w = np.blackman(self.numfreq)
         else:
             pos = int((cf - self.lowF) / self.deltaF)
             halfwidth = int(ww / (2.0 * self.deltaF))
             w = np.zeros(self.numfreq)
             w[pos - halfwidth:pos + halfwidth] = np.blackman(2 * halfwidth)
     self.data = self.data * w
     self.done()

    def test_calculate_lanczos_kernel(self):
        """
        Tests the kernels implemented in C against their numpy counterpart.
        """
        x = np.linspace(-5, 5, 11)

        values = calculate_lanczos_kernel(x, 5, "hanning")
        np.testing.assert_allclose(
            values["only_sinc"], np.sinc(x), atol=1E-9)
        np.testing.assert_allclose(
            values["only_taper"], np.hanning(len(x)), atol=1E-9)
        np.testing.assert_allclose(
            values["full_kernel"], np.sinc(x) * np.hanning(len(x)),
            atol=1E-9)

        values = calculate_lanczos_kernel(x, 5, "blackman")
        np.testing.assert_allclose(
            values["only_sinc"], np.sinc(x), atol=1E-9)
        np.testing.assert_allclose(
            values["only_taper"], np.blackman(len(x)), atol=1E-9)
        np.testing.assert_allclose(
            values["full_kernel"], np.sinc(x) * np.blackman(len(x)),
            atol=1E-9)

        values = calculate_lanczos_kernel(x, 5, "lanczos")
        np.testing.assert_allclose(
            values["only_sinc"], np.sinc(x), atol=1E-9)
        np.testing.assert_allclose(
            values["only_taper"], np.sinc(x / 5.0), atol=1E-9)
        np.testing.assert_allclose(
            values["full_kernel"], np.sinc(x) * np.sinc(x / 5.0),
            atol=1E-9)

 def __init__(self, fl=1200, fs=240, fftl=2048):
     self.fl = fl
     self.fs = fs
     self.fftl = fftl
     winpower = np.sqrt(np.sum(np.square(np.blackman(self.fl).astype(np.float32))))
     self.window = {'ana' : np.blackman(self.fl).astype(np.float32) / winpower,
                    'syn' : winpower / 0.42 / (self.fl / self.fs)}

    def _update_mode(self):
        if self.allocation.width <= 0:
            return

        if self.fft_show:
            max_freq = (self.freq_div * self.get_ticks())
            wanted_step = 1.0 / max_freq / 2 * self._input_freq
            self.input_step = max(floor(wanted_step), 1)

            self.draw_interval = 5.0

            self.set_max_samples(
                ceil(self.allocation.width / \
                        float(self.draw_interval) * 2) * self.input_step)

            # Create the (blackman) window
            self.fft_window = blackman(
                ceil(self.allocation.width / float(self.draw_interval) * 2))

            self.draw_interval *= wanted_step / self.input_step
        else:
            # Factor is just for triggering:
            time = (self.time_div * self.get_ticks())
            if time == 0:
                return
            samples = time * self._input_freq
            self.set_max_samples(samples * self.max_samples_fact)

            self.input_step = max(ceil(samples\
                                           / (self.allocation.width / 3.0)), 1)
            self.draw_interval = self.allocation.width\
                                           / (float(samples) / self.input_step)

            self.fft_window = None

def hpf():

    """
    Ref:
        https://tomroelandts.com/articles/how-to-create-a-simple-high-pass-filter

    Usage:
        h = hpf()
        s = np.convolve(s, h)
    """

    fc = 200.0 / 16000
    b = 150.0 / 16000
    N = int(np.ceil((4 / b)))
    if not N % 2: N += 1 # N, odd
    n = np.arange(N)

    # compute a low-pass filter
    h = np.sinc(2 * fc * (n - (N - 1) / 2.))
    w = np.blackman(N)
    h = h * w
    h = h / np.sum(h)

    # spectral inversion, make high-pass filter
    h = -h
    h[(N - 1) / 2] += 1

    return h

    def __init__(self, sample_rate, fft_peak=7.0, sample_format=None, search_size=1, verbose=False, signal_width=40e3):
        self._sample_rate = sample_rate
        self._fft_size=int(math.pow(2, 1+int(math.log(self._sample_rate/1000,2)))) # fft is approx 1ms long
        self._bin_size = float(self._fft_size)/self._sample_rate * 1000 # How many ms is one fft now?
        self._verbose = verbose
        self._search_size = search_size
        self._fft_peak = fft_peak

        if sample_format == "rtl":
            self._struct_elem = numpy.uint8
            self._struct_len = numpy.dtype(self._struct_elem).itemsize * self._fft_size *2
        elif sample_format == "hackrf":
            self._struct_elem = numpy.int8
            self._struct_len = numpy.dtype(self._struct_elem).itemsize * self._fft_size *2
        elif sample_format == "float":
            self._struct_elem = numpy.complex64
            self._struct_len = numpy.dtype(self._struct_elem).itemsize * self._fft_size
        else:
            raise Exception("No sample format given")

        self._window = numpy.blackman(self._fft_size)
        self._fft_histlen=500 # How many items to keep for moving average. 5 times our signal length
        self._data_histlen=self._search_size
        self._data_postlen=5
        self._signal_maxlen=1+int(30/self._bin_size) # ~ 30 ms
        self._fft_freq = numpy.fft.fftfreq(self._fft_size)
        self._signal_width=signal_width/(self._sample_rate/self._fft_size) # Area to ignore around an already found signal in Hz
        
        if self._verbose:
            print "fft_size=%d (=> %f ms)"%(self._fft_size,self._bin_size)
            print "calculate fft once every %d block(s)"%(self._search_size)
            print "require %.1f dB"%(10*math.log(self._fft_peak,10))
            print "signal_width: %d (= %.1f Hz)"%(self._signal_width,self._signal_width*self._sample_rate/self._fft_size)

 def apply_filter(self, array):
     if len(array.shape)!=1:
         return False # need 2D data
     nx = array.shape[0]
     # Apodization function
     if self.type == 1:
         ft = np.bartlett(nx)
     elif self.type == 2:
         ft = np.blackman(nx)
     elif self.type == 3:
         ft = np.hamming(nx)
     elif self.type == 4:
         ft = np.hanning(nx)
     elif self.type == 5:
         # Blackman-Harris window
         a0 = 0.35875
         a1 = 0.48829
         a2 = 0.14128
         a3 = 0.01168
         x = np.linspace(0,1,nx)
         ft = a0 - a1*np.cos(2*np.pi*x) + a2*np.cos(4*np.pi*x) + a3*np.cos(6*np.pi*x) 
     elif self.type == 6:
         # Lanczos window
         x = np.linspace(-1,1,nx)
         ft = np.sinc(x)
     elif self.type == 7:
         x = np.linspace(-1,1,nx)
         exec('x='+self.custom)
         ft = x
     else:
         ft = np.ones(nx)
     
     array.data = array.data*ft
     return True

def smavg():
    # 权重设置
    weights_exp=np.exp(np.linspace(-1,0,N))
    weights_exp/=weights_exp.sum()
    weights_sim=np.ones(N)
    weights_sim/=weights_sim.sum()
    weights_blw=np.blackman(N)
    weights_blw/=weights_blw.sum()

    # 平滑操作
    sm_exp=np.convolve(weights_exp,data)[N-1:-N+1]
    sm_sim=np.convolve(weights_sim ,data)[N-1:-N+1]
    sm_blw=np.convolve(weights_blw ,data)[N-1:-N+1]

    #np.savetxt('result',expa)
    exp_mean=1.01*np.average(sm_exp)
    sim_mean=1.02*np.average(sm_sim)

    #sub mean
    #expa=expa-exp_mean
    #sima=sima-sim_mean

    # 绘图
    #plt.plot(t,data[N-1:],lw=1.0)
    #plt.figure()
    plt.plot(t,data[N-1:],'k',label='origin')
    #plt.plot(t,sm_exp,'r',label='exp avg') #指数移动平均
    #plt.plot(t,sm_sim,'g',label='sim avg') #简单移动平均
    plt.plot(t,sm_blw,'y',lw=2.0,label='blackman win avg')  #布莱克曼窗函数
    plt.legend(loc='upper center')

    #plt.axhline(y=exp_mean,lw=2,color='y')
    #plt.axhline(y=sim_mean,lw=2,color='b')
    plt.show()

	def blackman_window(self, show=0):
		apod = N.blackman(self.data_points)
		for i in range(2):
			self.the_result.y[i] = self.the_result.y[i]*apod
		if show == 1:
			return self.the_result
		return self

    def Execute(self,data):

       #compute power in both components
       pReal = (data['real']/2.0**self.bits)
       pImag = (data['imag']/2.0**self.bits)

       # create empty complex array to combine pReal and pImag
       data = np.empty(pReal.shape, dtype=complex)
       data.real = pReal;
       data.imag = pImag;
       
       # subtract average dc from each range
       dc = np.mean(data,axis=0)
       data = np.subtract(data[:,],dc)
       data = np.clip(data,1e-15,np.max(data))

     
       coeff = np.blackman(64)
       for i in range(0,data.shape[1]):
          data[:,i] = np.convolve(data[:,i],coeff,'same')


       # compute complex 256-point for each range cell
       dMap = np.fft.fft(data,256, axis=0)
       dMap = np.fft.fftshift(dMap,axes=(0,))

       # convert values to dB
       dMap = 20*np.log10(np.abs(dMap));
    
       return dMap

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 get_freq2(sample, rate):
    sample_len = 22100
    frequencies = (np.arange(0, sample_len/2+1)).astype('float')*rate/sample_len
    
    
    avg_sample = np.convolve(sample, np.ones((5,)))
    avg_sample = avg_sample[2:-2]
    window = np.blackman(len(sample))
    fft_data4 = np.log(abs(np.fft.rfft(avg_sample*window, sample_len))**2)

    w = 50
    t = 3

    peaks = []    
    for j in range(0,len(fft_data4)-w):
        mx_i = np.argmax(fft_data4[j:j+w])
        strength = fft_data4[j+mx_i]-np.median(fft_data4[j:j+w])
        if mx_i == w/2 and strength > t:
           d_x = frequencies[j+mx_i]
           y0,y1,y2 = fft_data4[j+mx_i-1:j+mx_i+2]
           x1 = (y2 - y0) * .5 / (2 * y1 - y2 - y0)
           peaks.append([d_x, x1, strength])
           
    peaks = np.array(peaks)
    (freq, new_peaks) = freq_from_harmonics(peaks)
    return get_note(freq)  

    def get_wave_data(self, file_path):
        with wave.open(file_path, 'rb') as wf:
            sample_width = wf.getsampwidth()
            frame_rate = wf.getframerate()
            num_frames = wf.getnframes()

            length_in_seconds = num_frames / (frame_rate * 1.0)

            num_iterations = int(length_in_seconds * 2)
            iteration_size = int(num_frames / num_iterations)

            wave_data = []

            for fragment_num in range(0, num_iterations):
                data = wf.readframes(iteration_size)

                window_size = len(data) / sample_width
                window = np.blackman(window_size)

                fragment_data = wave.struct.unpack("%dh" % window_size, data) * window
                fft_data = abs(np.fft.rfft(fragment_data)) ** 2

                wave_data = np.hstack((wave_data, fft_data[:100]))  # take only 100 frequences

        return wave_data

 def _sweep_harmonics_responses(self, response, N, freqlist, shift=True):
     g = self.input_signal.ptp() / 2
     di = self._fft_convolve(self.sweep_inverse_signal, response / g)  # deconvolved impulse response
     imp_indices = numpy.zeros(
         N + 1, dtype=int
     )  # the indices of each impulse (1st linear, 2nd the first harm. dist. etc.)
     imp_indices[0] = len(self.sweep_inverse_signal) - 1
     for n in range(1, N + 1):
         imp_indices[n] = imp_indices[0] - int(round(math.log(n + 1) * self.sweep_rate))
     imps = []
     for n in range(N):
         dl = (imp_indices[n] - imp_indices[n + 1]) // 2
         lo = imp_indices[n] - dl  # the low limit where to cut
         hi = imp_indices[n] + dl  # the high limit where to cut
         di_w = di[lo:hi]
         win = numpy.blackman(len(di_w))
         if len(di.shape) == 2:
             win = win.reshape((len(win), 1))
         di_w *= win  # smoothed version
         w = freqlist
         if shift:
             w = w * (n + 1)
         h = scipy.signal.freqz(di_w, worN=w)[1]
         off = imp_indices[n] - lo
         h *= numpy.exp(off * 1j * freqlist)
         imps.append(h)
     return imps