def tone_est_near_index(sdata,index,range,sr):
samples = len(sdata)
fft_size = 2**int(floor(log(samples)/log(2.0)))
freq = fft(sdata[0:fft_size])
pdata = numpy.zeros(fft_size)
for i in xrange(fft_size): pdata[i] = abs(freq[i])
peak = 0
peak_index = 0
if (range == 0):
peak = pdata[index]
peak_index = index;
else:
for i in xrange(2*range):
if (pdata[index+i-range] > peak):
peak = pdata[index+i-range]
peak_index = index+i-range
R = peak*peak;
p = (freq[peak_index+1].real * freq[peak_index].real + freq[peak_index+1].imag * freq[peak_index].imag)/R
g = -p/(1.0-p)
q = (freq[peak_index-1].real * freq[peak_index].real + freq[peak_index-1].imag * freq[peak_index].imag)/R
e = q/(1.0-q)
if ((p>0) and (q>0)):
d = p
else:
d = q
u = peak_index + d
if (debug_estimates):
print "peak is at ",peak_index,"(",u,") and is ",peak
sum_phase = freq[peak_index-1]*c(-1,d) + freq[peak_index]*c(0,d) + freq[peak_index+1]*c(1,d)
sum_c_sq = abs(c(-1,d))*abs(c(-1,d)) + abs(c(0,d))*abs(c(0,d)) + abs(c(1,d))*abs(c(1,d))
amp = (abs(sum_phase)/sum_c_sq)/fft_size
phase_r = cmath.phase(sum_phase)
freq_est = u*sr/fft_size
return (amp,freq_est,phase_r)
python类fft()的实例源码
def zero_crossings(y_axis, window = 11):
"""
Algorithm to find zero crossings. Smoothens the curve and finds the
zero-crossings by looking for a sign change.
keyword arguments:
y_axis -- A list containg the signal over which to find zero-crossings
window -- the dimension of the smoothing window; should be an odd integer
(default: 11)
return -- the index for each zero-crossing
"""
# smooth the curve
length = len(y_axis)
x_axis = np.asarray(range(length), int)
# discard tail of smoothed signal
y_axis = _smooth(y_axis, window)[:length]
zero_crossings = np.where(np.diff(np.sign(y_axis)))[0]
indices = [x_axis[index] for index in zero_crossings]
# check if zero-crossings are valid
diff = np.diff(indices)
if diff.std() / diff.mean() > 0.2:
print(diff.std() / diff.mean())
print(np.diff(indices))
raise ValueError("False zero-crossings found, indicates problem {0} or {1}".format(
"with smoothing window", "problem with offset"))
# check if any zero crossings were found
if len(zero_crossings) < 1:
raise ValueError("No zero crossings found")
return indices
# used this to test the fft function's sensitivity to spectral leakage
#return indices + np.asarray(30 * np.random.randn(len(indices)), int)
############################Frequency calculation#############################
# diff = np.diff(indices)
# time_p_period = diff.mean()
#
# if diff.std() / time_p_period > 0.1:
# raise ValueError,
# "smoothing window too small, false zero-crossing found"
#
# #return frequency
# return 1.0 / time_p_period
##############################################################################
def extrBarkBands(samples):
'''
Extracts the energy of the lower 12 bark bands.
E. Zwicker and E. Terhardt: Analytical expressions for critical band
rate and critical bandwidth as a function of frequency. Journal of the
Acoustic Society of America, vol. 68, 1980, pp 1523-1525
:param samples: audio file samples with fs 44.1kHz
:return: array bb holding the bark band energies for each frame
'''
# PARAMETERS
fs = 44100
wSize = 1024
hSize = 1024
bands = [0.0, 50.0, 100.0, 150.0, 200.0, 300.0, 400.0, 510.0, 630.0, 770.0, 920.0, 1080.0, 1270.0]
fftSize = 1024
# INIT
window = hamming(wSize) # Hanning window
numFrames = int(math.floor(float(len(samples))/float(hSize)))
numBands = len(bands)-1
bb =[]
freqScale = (float(fs)*0.5) / float((wSize-1))
startBin = []
endBin = []
for ii in range(0,numBands):
startBin.append(int(round(bands[ii]/freqScale+0.5)))
endBin.append(int(round(bands[ii+1]/freqScale+0.5)))
# FRAME-WISE BARK BAND EXTRACTION
for i in range(0,numFrames):
frame = samples[i*hSize:i*hSize+wSize]
if len(frame)<wSize:
break
spec = abs(fft(frame*window)) / fftSize
b = []
for ii in range(0,numBands):
b.append(sum(spec[startBin[ii]:endBin[ii]]*spec[startBin[ii]:endBin[ii]]))
if max(b) > 0:
b = b/max(b)
bb.append(b)
return bb
def algoChannelSelection(left, right):
''' Algorithm which automatically selects the channel with dominant vocals from a stereo flamenco recording
based on spectral band energies as described in section 2-A-I of
Kroher, N. & Gomez, E. (2016). Automatic Transcription of Flamenco Singing from Polyphonic Music Recordings.
ACM / IEEE Transactions on Audio, Speech and Language Processing, 24(5), pp. 901-913.
:param left: samples of the left audio channel in 44.1kHz
:param right: samples of the right audio channel in 44.1kHz
:return: index of the dominant vocal channel (0 = left, 1 = right)
'''
# PARAMETERS
fs = 44100 # sample rate
wSize = 2048 # window size in samples
hSize = 2048 # hop size in samples
fftSize = 2048 # FFT size
freqGuitLow = 80.0 # lower bound for guitar band
freqGuitHigh = 400.0 # upper bound for guitar band
freqVocLow = 500.0 # lower bound for vocal band
freqVocHigh = 6000.0 # higher bound for vocal band
# INIT
window = hanning(wSize)
numFrames = int(math.floor(float(len(left))/float(wSize)))
# bin indices corresponding to freqeuncy band limits
indGuitLow = int(round((freqGuitLow/fs)*fftSize))
indGuitHigh = int(round((freqGuitHigh/fs)*fftSize))
indVocLow = int(round((freqVocLow/fs)*fftSize))
indVocHigh = int(round((freqVocHigh/fs)*fftSize))
# frame-wise computation of the spectral band ratio
sbrL = []
sbrR = []
for i in range(0,numFrames-100):
frameL = left[i*hSize:i*hSize+wSize]
specL = fft(frameL*window) / fftSize
specL = abs(specL * conj(specL))
guitMag = sum(specL[indGuitLow:indGuitHigh],0)
vocMag = sum(specL[indVocLow:indVocHigh],0)
sbrL.append(20*math.log10(vocMag/guitMag))
frameR = right[i*hSize:i*wSize+wSize]
specR = fft(frameR*window) / fftSize
specR = abs(specR * conj(specR))
guitMag = sum(specR[indGuitLow:indGuitHigh],0)
vocMag = sum(specR[indVocLow:indVocHigh],0)
sbrR.append(20*math.log10(vocMag/guitMag))
# select channel based on mean SBR
if mean(sbrL)>=mean(sbrR):
ind = 0
else:
ind = 1
return ind
def periodogram(self, signals, hold=False):
"""
Computes the estimation (in dB) for each epoch in a signal
Parameters
----------
signals : array, shape (n_epochs, n_points)
Signals from which one computes the power spectrum
hold : boolean, default = False
If True, the estimation is appended to the list of previous
estimations, else, the list is emptied and only the current
estimation is stored.
Returns
-------
psd : array, shape (n_epochs, n_freq)
Power spectrum estimated with a Welsh method on each epoch
n_freq = fft_length // 2 + 1
"""
fft_length, step = self.check_params()
signals = np.atleast_2d(signals)
n_epochs, n_points = signals.shape
block_length = min(self.block_length, n_points)
window = self.wfunc(block_length)
n_epochs, tmax = signals.shape
n_freq = fft_length // 2 + 1
psd = np.zeros((n_epochs, n_freq))
for i, sig in enumerate(signals):
block = np.arange(block_length)
# iterate on blocks
count = 0
while block[-1] < sig.size:
psd[i] += np.abs(
sp.fft(window * sig[block], fft_length, 0))[:n_freq] ** 2
count = count + 1
block = block + step
if count == 0:
raise IndexError(
'spectrum: first block has %d samples but sig has %d '
'samples' % (block[-1] + 1, sig.size))
# normalize
if self.donorm:
scale = 1.0 / (count * (np.sum(window) ** 2))
else:
scale = 1.0 / count
psd[i] *= scale
if not hold:
self.psd = []
self.psd.append(psd)
return psd
def run(eta, rem_add):
real_RTGF = np.loadtxt("rt_real.txt")
imag_RTGF = np.loadtxt("rt_imag.txt")
time_array = real_RTGF.transpose()[0]
npoints = time_array.shape[0]
delta_t = time_array[1]
frq = np.fft.fftfreq(npoints, delta_t)
frq = np.fft.fftshift(frq)*2.0*np.pi
real_part = real_RTGF.transpose()[1]
imag_part = imag_RTGF.transpose()[1]
if (rem_add == 'rem'):
fftinp = 1j*(real_part + 1j*imag_part)
elif (rem_add == 'add'):
fftinp = 1j*(real_part - 1j*imag_part)
else:
print 'Addition or Removal has not been specified!'
return
for i in range(npoints):
fftinp[i] = fftinp[i]*np.exp(-eta*time_array[i])
Y = fft(fftinp)
Y = np.fft.fftshift(Y)
Y_real = Y.real
Y_real = (Y_real*time_array[-1]/npoints)
Y_imag = Y.imag
Y_imag = (Y_imag*time_array[-1]/npoints)/np.pi
# Plot the results
with open('ldos.out', 'w') as fout:
fout.write('# Omega A(Omega)\n')
for i in range(npoints):
fout.write('%12.8f %12.8f\n' % (frq[i], Y_imag[i]))
with open('real_part.txt', 'w') as fout:
fout.write('# Omega A(Omega)\n')
for i in range(npoints):
fout.write('%12.8f %12.8f\n' % (frq[i], Y_real[i]))
