Python signal.blackmanharris函数代码示例

def THDN(signal, sample_rate, filename):

    signal = signal - mean(signal)  # this is so that the DC offset is not the peak in the case of PDM
    windowed = signal * blackmanharris(len(signal)) 

    

    # Find the peak of the frequency spectrum (fundamental frequency), and filter 
    # the signal by throwing away values between the nearest local minima
    f = rfft(windowed)

    #limit the bandwidth
    if len(f) > 24000:
       f[24000:len(f)] = 0

    bandwidth_limited = irfft(f)
    total_rms = rms_flat(bandwidth_limited)

    #for i in range(6000):
       #print abs(f[i])
    i = argmax(abs(f))

    #This will be exact if the frequency under test falls into the middle of an fft bin 
    print 'Frequency: %f Hz' % (sample_rate * (i / len(windowed)))
    lowermin, uppermin = find_range(abs(f), i)


    #notch out the signal
    f[lowermin: uppermin] = 0



    noise = irfft(f)
    THDN = rms_flat(noise) / total_rms
    print "THD+N:     %.12f%% or %.12f dB" % (THDN * 100, 20 * log10(THDN))

def freq_from_HPS(sig, fs):
    """
    Estimate frequency using harmonic product spectrum (HPS)
    
    """
    windowed = sig * blackmanharris(len(sig))

    from pylab import subplot, plot, log, copy, show

    #harmonic product spectrum:
    c = abs(rfft(windowed))
    maxharms = 8
    subplot(maxharms,1,1)
    plot(log(c))
    for x in range(2,maxharms):
        a = copy(c[::x]) #Should average or maximum instead of decimating
        # max(c[::x],c[1::x],c[2::x],...)
        c = c[:len(a)]
        i = argmax(abs(c))
        true_i = parabolic(abs(c), i)[0]
        print 'Pass %d: %f Hz' % (x, fs * true_i / len(windowed))
        c *= a
        subplot(maxharms,1,x)
        plot(log(c))
    show()

def window(f,start,stop,type='blackman'):
    """
    runs the data through a hamming window.
    @param f: The data matrix
    @param start: The start index of the hamming window.
    @param stop: The end index of the hamming window.
    """
    h=numpy.zeros(f.shape,dtype=float)

    if len(h.shape)==1:
        if type=='hamming':
            h[start:stop]=signal.hamming(stop-start)
        elif type=='blackman':
            h[start:stop]=signal.blackman(stop-start)
        elif type=='hann':
            h[start:stop]=signal.hann(stop-start)
        elif type=='blackmanharris':
            h[start:stop]=signal.blackmanharris(stop-start)
        elif type=='rectangular' or type=='rect' or type=='boxcar':
            h[start:stop]=signal.boxcar(stop-start)
    else:
        if type=='hamming':
            h[:,start:stop]=signal.hamming(stop-start)
        elif type=='blackman':
            h[:,start:stop]=signal.blackman(stop-start)
        elif type=='hann':
            h[:,start:stop]=signal.hann(stop-start)
        elif type=='blackmanharris':
            h[:,start:stop]=signal.blackmanharris(stop-start)
        elif type=='rectangular' or type=='rect' or type=='boxcar':
            h[:,start:stop]=signal.boxcar(stop-start)
    return numpy.multiply(f,h)

def freq_from_fft(signal, fs):
    # Compute Fourier transform of windowed signal
    windowed = signal * blackmanharris(len(signal))
    f = rfft(windowed)
    i = argmax(abs(f))
    
    # Convert to equivalent frequency
    return fs*i/len(windowed)

def delayTransformVisibilities(visList,df,kernel=None,wind='Blackman-Harris'):
    '''
     ARGS:
     visList: nfreqxnvis complex matrix of visibilities
     df: float indicating channel width
     RETURNS:
     ndelayxnvis grid of delay-transformed visibilities
    '''
    nf=visList.shape[0]
    if kernel is None:
        kernel=np.ones(nf)

    nv=visList.shape[1]
    windowGrid=np.zeros_like(visList)
    if wind=='Blackman-Harris':
        window=signal.blackmanharris(nf)
    elif wind=='Nithya':
        window=signal.blackmanharris(nf/2)
        window=signal.fftconvolve(window,window,mode='full')
        window=np.append(window,window[-1])
    window/=window.max()
    window/=np.sqrt(np.mean(np.abs(kernel*window)**2.)) 
    for vNum in range(nv):
        windowGrid[:,vNum]=window*kernel
    return df*fft.fftshift(fft.fft(fft.fftshift(visList*windowGrid,axes=[0]),axis=0),axes=[0])

 def test_basic(self):
     assert_allclose(signal.blackmanharris(6, False),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645])
     assert_allclose(signal.blackmanharris(6),
                     [6.0e-05, 0.1030114893456638, 0.7938335106543362,
                      0.7938335106543364, 0.1030114893456638, 6.0e-05])
     assert_allclose(signal.blackmanharris(7, sym=True),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645,
                      6.0e-05])

 def test_basic(self):
     assert_allclose(signal.blackmanharris(6, False),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645])
     assert_allclose(signal.blackmanharris(7, sym=False),
                     [6.0e-05, 0.03339172347815117, 0.332833504298565,
                      0.8893697722232837, 0.8893697722232838,
                      0.3328335042985652, 0.03339172347815122])
     assert_allclose(signal.blackmanharris(6),
                     [6.0e-05, 0.1030114893456638, 0.7938335106543362,
                      0.7938335106543364, 0.1030114893456638, 6.0e-05])
     assert_allclose(signal.blackmanharris(7, sym=True),
                     [6.0e-05, 0.055645, 0.520575, 1.0, 0.520575, 0.055645,
                      6.0e-05])

def sineModelSynth(tfreq, tmag, tphase, N, H, fs):
    """
    Synthesis of a sound using the sinusoidal model
    tfreq, tmag, tphase: frequencies, magnitudes and phases of sinusoids
    N: synthesis FFT size, H: hop size, fs: sampling rate
    returns y: output array sound
    """

    hN = N / 2 # half of FFT size for synthesis
    L = tfreq.shape[0] # number of frames
    pout = 0
    ysize = H * (L+3)
    y = np.zeros(ysize)
    # synthesis window
    sw = np.zeros(N)
    ow = triang(2 * H)
    sw[hN-H:hN+H] = ow
    bh = blackmanharris(N)
    bh = bh / sum(bh)
    sw[hN-H:hN+H] = sw[hN-H:hN+H] / bh[hN-H:hN+H]
    lastytfreq = tfreq[0,:]
    ytphase = 2*np.pi*np.random.rand(tfreq[0,:].size)
    for l in range(L): # iterate over all frames
        if (tphase.size > 0):
            ytphase = tphase[l, :]

def freq_power_from_fft(sig, fs):
	# Compute Fourier transform of windowed signal
	windowed = sig * blackmanharris(CHUNK)
	fftData = abs(rfft(windowed))
	# powerData = 20*log10(fftData)
	# Find the index of the peak and interpolate to get a more accurate peak
	i = argmax(fftData) # Just use this for less-accurate, naive version
	# make sure i is not an endpoint of the bin
	if i != len(fftData)-1 and i != 0:
		#print("data: " + str(parabolic(log(abs(fftData)), i)))
		true_i,logmag = parabolic(log(abs(fftData)), i)
		# Convert to equivalent frequency
		freq=  fs * true_i / len(windowed)
		#print("frequency="+ str(freq) + " log of magnitude=" + str(logmag))
		if logmag < 0:
			logmag = 0	
		return freq,logmag
	
	else:
		freq = fs * i / len(windowed)
		logmag = log(abs(fftData))[i]
		if logmag < 0:
			logmag = 0
		#print("frequency="+ str(freq) + "log of magnitude not interp=" + str(logmag))
		return freq,logmag

    def run(self, tf, dt, K, c1, c2):
        """ Run a simulation """
        n, m, k = self.n, self.m, self.k

        # Total simulation time
        simTime = int(tf / dt)

        # Returns the three synaptic connections kernels
        W12, W21, W22, delays = self.build_kernels(K)

        # Compute delays by dividing distances by axonal velocity
        delays12 = np.floor(delays[0] / c2)
        delays21 = np.floor(delays[1] / c1)
        delays22 = np.floor(delays[2] / c2)
        maxDelay = int(max(delays12[0].max(), delays21[0].max(), delays22[0].max()))

        # Set the initial conditions and the history
        self.initial_conditions(simTime)

        # Initialize the cortical and striatal inputs
        Cx = 0.5
        Str = 0.4

        # Presynaptic activities
        pre12, pre21, pre22 = np.empty((m,)), np.empty((m,)), np.empty((m,))

        # Simulation
        for i in range(maxDelay, simTime):
            # Take into account the history of rate for each neuron according
            # to its axonal delay
            for idxi, ii in enumerate(range(m)):
                mysum = 0.0
                for jj in range(k, n):
                    mysum += (W12[ii, jj] * self.X2[i - delays12[ii, jj], jj]) * self.dx
                pre12[idxi] = mysum

            for idxi, ii in enumerate(range(k, n)):
                mysum = 0.0
                for jj in range(0, m):
                    mysum += (W21[ii, jj] * self.X1[i - delays21[ii, jj], jj]) * self.dx
                pre21[idxi] = mysum

            for idxi, ii in enumerate(range(k, n)):
                mysum = 0.0
                for jj in range(k, n):
                    mysum += (W22[ii, jj] * self.X2[i - delays22[ii, jj], jj]) * self.dx
                pre22[idxi] = mysum

            # Forward Euler step
            self.X1[i, :m] = self.X1[i - 1, :m] + (-self.X1[i - 1, :m] + self.S1(-pre12 + Cx)) * dt / self.tau1
            self.X2[i, k:] = self.X2[i - 1, k:] + (-self.X2[i - 1, k:] + self.S2(pre21 - pre22 - Str)) * dt / self.tau2
        dx = 1.0 / float(m)
        fr = self.X1.sum(axis=1) * dx / 1.0

        signal = detrend(fr)
        windowed = signal * blackmanharris(len(signal))
        f = rfft(windowed)
        i = np.argmax(np.abs(f))
        # true_i = parabolic(np.log(np.abs(f)), i)[0]
        return i

def _BLACKMANHARRIS_Window(signal):
    """ blackman-harris window process
        Use Windows to smooth transient data in the begin/end
    """
    windowed = signal * blackmanharris(len(signal))  # TODO Kaiser?

    return windowed

def stochasticResidualAnal(x, N, H, sfreq, smag, sphase, fs, stocf):
	"""
	Subtract sinusoids from a sound and approximate the residual with an envelope
	x: input sound, N: fft size, H: hop-size
	sfreq, smag, sphase: sinusoidal frequencies, magnitudes and phases
	fs: sampling rate; stocf: stochastic factor, used in the approximation
	returns stocEnv: stochastic approximation of residual
	"""

	hN = N/2                                              # half of fft size
	x = np.append(np.zeros(hN),x)                         # add zeros at beginning to center first window at sample 0
	x = np.append(x,np.zeros(hN))                         # add zeros at the end to analyze last sample
	bh = blackmanharris(N)                                # synthesis window
	w = bh/ sum(bh)                                       # normalize synthesis window
	L = sfreq.shape[0]                                    # number of frames, this works if no sines
	pin = 0
	for l in range(L):
		xw = x[pin:pin+N] * w                               # window the input sound
		X = fft(fftshift(xw))                               # compute FFT
		Yh = UF_C.genSpecSines(N*sfreq[l,:]/fs, smag[l,:], sphase[l,:], N)   # generate spec sines
		Xr = X-Yh                                           # subtract sines from original spectrum
		mXr = 20*np.log10(abs(Xr[:hN]))                     # magnitude spectrum of residual
		mXrenv = resample(np.maximum(-200, mXr), mXr.size*stocf)  # decimate the mag spectrum
		if l == 0:                                          # if first frame
			stocEnv = np.array([mXrenv])
		else:                                               # rest of frames
			stocEnv = np.vstack((stocEnv, np.array([mXrenv])))
		pin += H                                            # advance sound pointer
	return stocEnv

def sineSubtraction(x, N, H, sfreq, smag, sphase, fs):
	"""
	Subtract sinusoids from a sound
	x: input sound, N: fft-size, H: hop-size
	sfreq, smag, sphase: sinusoidal frequencies, magnitudes and phases
	returns xr: residual sound
	"""

	hN = N/2                                           # half of fft size
	x = np.append(np.zeros(hN),x)                      # add zeros at beginning to center first window at sample 0
	x = np.append(x,np.zeros(hN))                      # add zeros at the end to analyze last sample
	bh = blackmanharris(N)                             # blackman harris window
	w = bh/ sum(bh)                                    # normalize window
	sw = np.zeros(N)                                   # initialize synthesis window
	sw[hN-H:hN+H] = triang(2*H) / w[hN-H:hN+H]         # synthesis window
	L = sfreq.shape[0]                                 # number of frames, this works if no sines
	xr = np.zeros(x.size)                              # initialize output array
	pin = 0
	for l in range(L):
		xw = x[pin:pin+N]*w                              # window the input sound
		X = fft(fftshift(xw))                            # compute FFT
		Yh = UF_C.genSpecSines(N*sfreq[l,:]/fs, smag[l,:], sphase[l,:], N)   # generate spec sines
		Xr = X-Yh                                        # subtract sines from original spectrum
		xrw = np.real(fftshift(ifft(Xr)))                # inverse FFT
		xr[pin:pin+N] += xrw*sw                          # overlap-add
		pin += H                                         # advance sound pointer
	xr = np.delete(xr, range(hN))                      # delete half of first window which was added in stftAnal
	xr = np.delete(xr, range(xr.size-hN, xr.size))     # delete half of last window which was added in stftAnal
	return xr

def sinewaveSynth(freq, mag, N, H, fs):
  # Synthesis of a time-varying sinusoid
  # freq,mag, phase: frequency, magnitude and phase of sinusoid,
  # N: synthesis FFT size, H: hop size, fs: sampling rate
  # returns y: output array sound
  hN = N/2                                                # half of FFT size for synthesis
  L = freq.size                                           # number of frames
  pout = 0                                                # initialize output sound pointer         
  ysize = H*(L+3)                                         # output sound size
  y = np.zeros(ysize)                                     # initialize output array
  sw = np.zeros(N)                                        # initialize synthesis window
  ow = triang(2*H);                                       # triangular window
  sw[hN-H:hN+H] = ow                                      # add triangular window
  bh = blackmanharris(N)                                  # blackmanharris window
  bh = bh / sum(bh)                                       # normalized blackmanharris window
  sw[hN-H:hN+H] = sw[hN-H:hN+H]/bh[hN-H:hN+H]             # normalized synthesis window
  lastfreq = freq[0]                                      # initialize synthesis frequencies
  phase = 0                                               # initialize synthesis phases 
  for l in range(L):                                      # iterate over all frames
    phase += (np.pi*(lastfreq+freq[l])/fs)*H            # propagate phases
    Y = UF.genSpecSines(freq[l], mag[l], phase, N, fs)    # generate sines in the spectrum         
    lastfreq = freq[l]                                    # save frequency for phase propagation
    yw = np.real(fftshift(ifft(Y)))                       # compute inverse FFT
    y[pout:pout+N] += sw*yw                               # overlap-add and apply a synthesis window
    pout += H                                             # advance sound pointer
  y = np.delete(y, range(hN))                             # delete half of first window
  y = np.delete(y, range(y.size-hN, y.size))              # delete half of the last window 
  return y

def plot_specgram(ax, data, fs, nfft=256, noverlap=128, window='hann',
                  cmap='jet', interpolation='bilinear', rasterized=True):

    if window not in SPECGRAM_WINDOWS:
        raise ValueError("Window not supported")

    elif window == "boxcar":
        mwindow = signal.boxcar(nfft)
    elif window == "hamming":
        mwindow = signal.hamming(nfft)
    elif window == "hann":
        mwindow = signal.hann(nfft)
    elif window == "bartlett":
        mwindow = signal.bartlett(nfft)
    elif window == "blackman":
        mwindow = signal.blackman(nfft)
    elif window == "blackmanharris":
        mwindow = signal.blackmanharris(nfft)

    specgram, freqs, time = mlab.specgram(data, NFFT=nfft, Fs=fs,
                                          window=mwindow,
                                          noverlap=noverlap)
    specgram = 10 * np.log10(specgram[1:, :])
    specgram = np.flipud(specgram)

    freqs = freqs[1:]
    halfbin_time = (time[1] - time[0]) / 2.0
    halfbin_freq = (freqs[1] - freqs[0]) / 2.0
    extent = (time[0] - halfbin_time, time[-1] + halfbin_time,
              freqs[0] - halfbin_freq, freqs[-1] + halfbin_freq)

    ax.imshow(specgram, cmap=cmap, interpolation=interpolation,
                            extent=extent, rasterized=rasterized)
    ax.axis('tight')