# AUTHORS
# Jose PATINO, EURECOM, Sophia-Antipolis, France, 2019
# http://www.eurecom.fr/en/people/patino-jose
# Contact: patino[at]eurecom[dot]fr, josempatinovillar[at]gmail[dot]com
# DIARIZATION FUNCTIONS
import numpy as np
def extractFeatures(audioFile,framelength,frameshift,nfilters,ncoeff):
import librosa
y, sr = librosa.load(audioFile,sr=None)
frame_length_inSample=framelength*sr
hop = int(frameshift*sr)
NFFT=int(2**np.ceil(np.log2(frame_length_inSample)))
if sr>=16000:
features=librosa.feature.mfcc(y=y, sr=sr,dct_type=2,n_mfcc=ncoeff,n_mels=nfilters,n_fft=NFFT,hop_length=hop,fmin=20,fmax=7600).T
else:
features=librosa.feature.mfcc(y=y, sr=sr,dct_type=2,n_mfcc=ncoeff,n_mels=nfilters,n_fft=NFFT,hop_length=hop).T
return features
def getFeatures(featureFile):
features = HTKFile()
features.load(featureFile)
allData = np.asarray(features.data)
return allData
def readUEMfile(path,filename,ext,nFeatures,frameshift):
uemFile = path + filename + ext
notEvaluatedFramesMask=np.zeros([1, nFeatures])
f=open(uemFile,'r')
C=f.read().splitlines()
f.close()
initTime = []
endTime = []
if C:
if len(C[0])>1:
for idx,x in enumerate(C):
initTime=np.append(initTime,float(x.strip().split(' ')[2]))
endTime=np.append(endTime,float(x.strip().split(' ')[3]))
else:
print('UEM annotations file was empty. The complete file is considered evaluable.')
if type(initTime) is list:
notEvaluatedFramesMask = np.ones([1,nFeatures])
else:
initTime = np.floor(initTime/frameshift)
endTime = np.floor(endTime/frameshift)
for idxT,x in np.ndenumerate(initTime):
it=initTime[idxT]
if endTime[idxT]>nFeatures:
et = nFeatures
else:
et = endTime[idxT]
notEvaluatedFramesMask[0,int(it):int(et)]=1
return notEvaluatedFramesMask
def readSADfile(path,filename,ext,nFeatures,frameshift, format):
sadFile = path + filename + ext
notEvaluatedFramesMask=np.zeros([1, nFeatures])
f=open(sadFile,'r')
C=f.read().splitlines()
f.close()
initTime = []
endTime = []
if C:
if len(C[0])>1:
if format=='LBL':
for idx,x in enumerate(C):
initTime=np.append(initTime,float(x.strip().split(' ')[0]))
endTime=np.append(endTime,float(x.strip().split(' ')[1]))
elif format=='RTTM':
for idx,x in enumerate(C):
initTime=np.append(initTime,float(x.strip().split(' ')[3]))
endTime=np.append(endTime,float(x.strip().split(' ')[3])+float(x.strip().split(' ')[4]))
elif format=='MDTM':
for idx,x in enumerate(C):
initTime=np.append(initTime,float(x.strip().split(' ')[2]))
endTime=np.append(endTime,float(x.strip().split(' ')[2])+float(x.strip().split(' ')[3]))
else:
print('SAD annotations file was empty. The complete file is considered speech.')
if type(initTime) is list:
notEvaluatedFramesMask = np.ones([1,nFeatures])
else:
initTime = np.round(initTime/frameshift)
endTime = np.round(endTime/frameshift)
for idxT,x in np.ndenumerate(initTime):
it=initTime[idxT]
if endTime[idxT]>nFeatures:
et = nFeatures
else:
et = endTime[idxT]
notEvaluatedFramesMask[0,int(it):int(et)]=1
return notEvaluatedFramesMask
def getSADfile(config,filename,nFeatures):
""" The following VAD applying procedure is adapted from the repository
provided as part of the DIHARD II challenge speech enhancement system:
https://github.com/staplesinLA/denoising_DIHARD18/blob/master/main_get_vad.py"""
import os
from librosa import load
data, sr = load(config['PATH']['audio']+filename+config['EXTENSION']['audio'],sr=None)
va_framed = py_webrtcvad(data, fs=sr, fs_vad=sr, hoplength=30, vad_mode=0)
segments = get_py_webrtcvad_segments(va_framed,sr)
if not os.path.isdir(config['PATH']['SAD']):
os.mkdir(config['PATH']['SAD'])
if os.path.isfile(config['PATH']['SAD']+filename+'.lab'):
os.remove(config['PATH']['SAD']+filename+'.lab')
output_file = open(config['PATH']['SAD']+filename+'.lab','w')
for i in range(segments.shape[0]):
start_time = segments[i][0]
end_time = segments[i][1]
output_file.write("%.3f %.3f speech\n" % (start_time, end_time))
output_file.close()
maskSAD = readSADfile(config['PATH']['SAD'],filename,'.lab',nFeatures,config.getfloat('FEATURES','frameshift'),'LBL')
return maskSAD
def py_webrtcvad(data, fs, fs_vad, hoplength=30, vad_mode=0):
import webrtcvad
from librosa.core import resample
from librosa.util import frame
""" Voice activity detection.
This was implementioned for easier use of py-webrtcvad.
Thanks to: https://github.com/wiseman/py-webrtcvad.git
Parameters
----------
data : ndarray
numpy array of mono (1 ch) speech data.
1-d or 2-d, if 2-d, shape must be (1, time_length) or (time_length, 1).
if data type is int, -32768 < data < 32767.
if data type is float, -1 < data < 1.
fs : int
Sampling frequency of data.
fs_vad : int, optional
Sampling frequency for webrtcvad.
fs_vad must be 8000, 16000, 32000 or 48000.
Default is 16000.
hoplength : int, optional
Step size[milli second].
hoplength must be 10, 20, or 30.
Default is 0.1.
vad_mode : int, optional
set vad aggressiveness.
As vad_mode increases, it becomes more aggressive.
vad_mode must be 0, 1, 2 or 3.
Default is 0.
Returns
-------
vact : ndarray
voice activity. time length of vact is same as input data.
If 0, it is unvoiced, 1 is voiced.
"""
# check argument
if fs_vad not in [8000, 16000, 32000, 48000]:
raise ValueError('fs_vad must be 8000, 16000, 32000 or 48000.')
if hoplength not in [10, 20, 30]:
raise ValueError('hoplength must be 10, 20, or 30.')
if vad_mode not in [0, 1, 2, 3]:
raise ValueError('vad_mode must be 0, 1, 2 or 3.')
# check data
if data.dtype.kind == 'i':
if data.max() > 2**15 - 1 or data.min() < -2**15:
raise ValueError(
'when data type is int, data must be -32768 < data < 32767.')
data = data.astype('f')
elif data.dtype.kind == 'f':
if np.abs(data).max() >= 1:
data = data / np.abs(data).max() * 0.9
warnings.warn('input data was rescaled.')
data = (data * 2**15).astype('f')
else:
raise ValueError('data dtype must be int or float.')
data = data.squeeze()
if not data.ndim == 1:
raise ValueError('data must be mono (1 ch).')
# resampling
if fs != fs_vad:
resampled = resample(data, fs, fs_vad)
else:
resampled = data
resampled = resampled.astype('int16')
hop = fs_vad * hoplength // 1000
framelen = resampled.size // hop + 1
padlen = framelen * hop - resampled.size
paded = np.lib.pad(resampled, (0, padlen), 'constant', constant_values=0)
framed = frame(paded, frame_length=hop, hop_length=hop).T
vad = webrtcvad.Vad()
vad.set_mode(vad_mode)
valist = [vad.is_speech(tmp.tobytes(), fs_vad) for tmp in framed]
hop_origin = fs * hoplength // 1000
va_framed = np.zeros([len(valist), hop_origin])
va_framed[valist] = 1
return va_framed.reshape(-1)[:data.size]
def get_py_webrtcvad_segments(vad_info,fs):
vad_index = np.where( vad_info==1.0) # find the speech index
vad_diff = np.diff( vad_index)
vad_temp = np.zeros_like(vad_diff)
vad_temp[ np.where(vad_diff==1) ] = 1
vad_temp = np.column_stack( (np.array([0]), vad_temp, np.array([0]) ))
final_index = np.diff(vad_temp)
starts = np.where( final_index == 1)
ends = np.where( final_index == -1)
sad_info = np.column_stack( (starts[1], ends[1]) )
vad_index = vad_index[0]
segments = np.zeros_like(sad_info,dtype=np.float)
for i in range(sad_info.shape[0]):
segments[i][0] = float( vad_index [ sad_info[i][0] ] ) / fs
segments[i][1] = float( vad_index[ sad_info[i][1] ] +1 ) / fs
return segments # present in seconds
def getSegmentTable(mask, speechMapping, wLength, wIncr, wShift):
changePoints,segBeg,segEnd,nSegs = unravelMask(mask)
segmentTable = np.empty([0,4])
for i in range(nSegs):
begs = np.arange(segBeg[i],segEnd[i],wShift)
bbegs = np.maximum(segBeg[i],begs-wIncr)
ends = np.minimum(begs+wLength-1,segEnd[i])
eends = np.minimum(ends+wIncr,segEnd[i])
segmentTable=np.vstack((segmentTable,np.vstack((bbegs, begs, ends, eends)).T))
return segmentTable
def unravelMask(mask):
changePoints = np.diff(1*mask)
segBeg=np.where(changePoints==1)[0]+1
segEnd=np.where(changePoints==-1)[0]
if mask[0]==1:
segBeg = np.insert(segBeg,0,0)
if mask[-1]==1:
segEnd = np.append(segEnd,np.size(mask)-1)
nSegs = np.size(segBeg)
return changePoints,segBeg,segEnd,nSegs
def trainKBM(data, windowLength, windowRate, kbmSize):
from scipy.stats import multivariate_normal
from scipy.spatial.distance import cdist
# Calculate number of gaussian components in the whole gaussian pool
numberOfComponents = int(np.floor((np.size(data,0)-windowLength)/windowRate))
# Add new array for storing the mvn objects
gmPool = []
likelihoodVector = np.zeros((numberOfComponents, 1))
muVector = np.zeros((numberOfComponents,np.size(data,1)))
sigmaVector = np.zeros((numberOfComponents,np.size(data,1)))
for i in range(numberOfComponents):
mu = np.mean(data[np.arange((i*windowRate),(i*windowRate+windowLength),1,int)],axis=0)
std = np.std(data[np.arange((i*windowRate),(i*windowRate+windowLength),1,int)],axis=0)
muVector[i], sigmaVector[i] = mu, std
mvn = multivariate_normal(mu,std)
gmPool.append(mvn)
likelihoodVector[i] = -np.sum(mvn.logpdf(data[np.arange((i*windowRate),(i*windowRate+windowLength),1,int)]))
# Define the global dissimilarity vector
v_dist = np.inf*np.ones((numberOfComponents,1))
# Create the kbm itself, which is a vector of kbmSize size, and contains the gaussian IDs of the components
kbm = np.zeros((kbmSize,1));
# As the stored likelihoods are negative, get the minimum likelihood
bestGaussianID = np.where(likelihoodVector==np.min(likelihoodVector))[0]
currentGaussianID = bestGaussianID
kbm[0]=currentGaussianID
v_dist[currentGaussianID]=-np.inf
# Compare the current gaussian with the remaining ones
dpairsAll = cdist(muVector,muVector,metric='cosine')
np.fill_diagonal(dpairsAll,-np.inf)
for j in range(1,kbmSize):
dpairs = dpairsAll[currentGaussianID]
v_dist=np.minimum(v_dist,dpairs.T)
# Once all distances are computed, get the position with highest value
# set this position to 1 in the binary KBM ponemos a 1 en el vector kbm
# store the gaussian ID in the KBM
currentGaussianID = np.where(v_dist==np.max(v_dist))[0]
kbm[j]=currentGaussianID
v_dist[currentGaussianID]=-np.inf
return [kbm, gmPool]
def getVgMatrix(data, gmPool, kbm, topGaussiansPerFrame):
print('Calculating log-likelihood table... ',end='')
logLikelihoodTable = getLikelihoodTable(data,gmPool,kbm)
print('done')
print('Calculating Vg matrix... ',end='')
Vg = np.argsort(-logLikelihoodTable)[:,0:topGaussiansPerFrame]
print('done')
return Vg
def getLikelihoodTable(data,gmPool,kbm):
# GETLIKELIHOODTABLE computes the log-likelihood of each feature in DATA
# against all the Gaussians of GMPOOL specified by KBM vector
# Inputs:
# DATA = matrix of feature vectors
# GMPOOL = pool of Gaussians of the kbm model
# KBM = vector of the IDs of the actual Gaussians of the KBM
# Output:
# LOGLIKELIHOODTABLE = NxM matrix storing the log-likelihood of each of
# the N features given each of th M Gaussians in the KBM
kbmSize = np.size(kbm,0)
logLikelihoodTable = np.zeros([np.size(data,0),kbmSize])
for i in range(kbmSize):
logLikelihoodTable[:,i]=gmPool[int(kbm[i])].logpdf(data)
return logLikelihoodTable
def getSegmentBKs(segmentTable, kbmSize, Vg, bitsPerSegmentFactor, speechMapping):
# GETSEGMENTBKS converts each of the segments in SEGMENTTABLE into a binary key
# and/or cumulative vector.
# Inputs:
# SEGMENTTABLE = matrix containing temporal segments returned by 'getSegmentTable' function
# KBMSIZE = number of components in the kbm model
# VG = matrix of the top components per frame returned by 'getVgMatrix' function
# BITSPERSEGMENTFACTOR = proportion of bits that will be set to 1 in the binary keys
# Output:
# SEGMENTBKTABLE = NxKBMSIZE matrix containing N binary keys for each N segments in SEGMENTTABLE
# SEGMENTCVTABLE = NxKBMSIZE matrix containing N cumulative vectors for each N segments in SEGMENTTABLE
numberOfSegments = np.size(segmentTable,0)
segmentBKTable = np.zeros([numberOfSegments,kbmSize])
segmentCVTable = np.zeros([numberOfSegments,kbmSize])
for i in range(numberOfSegments):
# Conform the segment according to the segmentTable matrix
beginningIndex = int(segmentTable[i,0])
endIndex = int(segmentTable[i,3])
# Store indices of features of the segment
# speechMapping is substracted one because 1-indexing is used for this variable
A = np.arange(speechMapping[beginningIndex]-1,speechMapping[endIndex],dtype=int)
segmentBKTable[i], segmentCVTable[i] = binarizeFeatures(kbmSize, Vg[A,:], bitsPerSegmentFactor)
print('done')
return segmentBKTable, segmentCVTable
def binarizeFeatures(binaryKeySize, topComponentIndicesMatrix, bitsPerSegmentFactor):
# BINARIZEMATRIX Extracts a binary key and a cumulative vector from the the
# rows of VG specified by vector A
# Inputs:
# BINARYKEYSIZE = binary key size
# TOPCOMPONENTINDICESMATRIX = matrix of top Gaussians per frame
# BITSPERSEGMENTFACTOR = Proportion of positions of the binary key which will be set to 1
# Output:
# BINARYKEY = 1xBINARYKEYSIZE binary key
# V_F = 1xBINARYKEYSIZE cumulative vector
numberOfElementsBinaryKey = np.floor(binaryKeySize * bitsPerSegmentFactor)
# Declare binaryKey
binaryKey = np.zeros([1, binaryKeySize])
# Declare cumulative vector v_f
v_f = np.zeros([1, binaryKeySize])
unique, counts = np.unique(topComponentIndicesMatrix, return_counts=True)
# Fill CV
v_f[:,unique]=counts
# Fill BK
binaryKey[0,np.argsort(-v_f)[0][0:int(numberOfElementsBinaryKey)]]=1
# CV normalization
v_f = v_f/np.sum(v_f)
return binaryKey, v_f
def performClusteringLinkage(segmentBKTable, segmentCVTable, N_init, linkageCriterion,linkageMetric ):
from scipy.cluster.hierarchy import linkage
from scipy import cluster
if linkageMetric == 'jaccard':
observations = segmentBKTable
elif linkageMetric == 'cosine':
observations = segmentCVTable
else:
observations = segmentCVTable
clusteringTable = np.zeros([np.size(segmentCVTable,0),N_init])
Z = linkage(observations,method=linkageCriterion,metric=linkageMetric)
for i in np.arange(N_init):
clusteringTable[:,i] = cluster.hierarchy.cut_tree(Z,N_init-i).T+1
k=N_init
print('done')
return clusteringTable, k
def performClustering( speechMapping, segmentTable, segmentBKTable, segmentCVTable, Vg, bitsPerSegmentFactor, kbmSize, N_init, initialClustering, clusteringMetric):
numberOfSegments = np.size(segmentTable,0)
clusteringTable = np.zeros([numberOfSegments, N_init])
finalClusteringTable = np.zeros([numberOfSegments, N_init])
activeClusters = np.ones([N_init,1])
clustersBKTable = np.zeros([N_init, kbmSize])
clustersCVTable = np.zeros([N_init, kbmSize])
clustersBKTable, clustersCVTable = calcClusters(clustersCVTable,clustersBKTable,activeClusters,initialClustering,N_init,segmentTable,kbmSize,speechMapping,Vg,bitsPerSegmentFactor)
####### Here the clustering algorithm begins. Steps are:
####### 1. Reassign all data among all existing signatures and retrain them
####### using the new clustering
####### 2. Save the resulting clustering solution
####### 3. Compare all signatures with each other and merge those two with
####### highest similarity, creating a new signature for the resulting
####### cluster
####### 4. Back to 1 if #clusters > 1
for k in range(N_init):
####### 1. Data reassignment. Calculate the similarity between the current segment with all clusters and assign it to the one which maximizes
####### the similarity. Finally re-calculate binaryKeys for all cluster
# before doing anything, check if there are remaining clusters
# if there is only one active cluster, break
if np.sum(activeClusters)==1:
break
clustersStillActive=np.zeros([1,N_init])
segmentToClustersSimilarityMatrix = binaryKeySimilarity_cdist(clusteringMetric,segmentBKTable,segmentCVTable,clustersBKTable,clustersCVTable)
# clusteringTable[:,k] = finalClusteringTable[:,k] = np.argmax(segmentToClustersSimilarityMatrix,axis=1)+1
clusteringTable[:,k] = finalClusteringTable[:,k] = np.nanargmax(segmentToClustersSimilarityMatrix,axis=1)+1
# clustersStillActive[:,np.unique(clusteringTable[:,k]).astype(int)-1] = 1
clustersStillActive[:,np.unique(clusteringTable[:,k]).astype(int)-1] = 1
####### update all binaryKeys for all new clusters
activeClusters = clustersStillActive
clustersBKTable, clustersCVTable = calcClusters(clustersCVTable,clustersBKTable,activeClusters.T,clusteringTable[:,k].astype(int),N_init,segmentTable,kbmSize,speechMapping,Vg,bitsPerSegmentFactor)
####### 2. Compare all signatures with each other and merge those two with highest similarity, creating a new signature for the resulting
clusterSimilarityMatrix = binaryKeySimilarity_cdist(clusteringMetric,clustersBKTable,clustersCVTable,clustersBKTable,clustersCVTable)
np.fill_diagonal(clusterSimilarityMatrix,np.nan)
value = np.nanmax(clusterSimilarityMatrix)
location = np.nanargmax(clusterSimilarityMatrix)
R,C = np.unravel_index(location,(N_init,N_init))
### Then we merge clusters R and C
#print('Merging clusters',R+1,'and',C+1,'with a similarity score of',np.around(value,decimals=4))
print('Merging clusters','%3s'%str(R+1),'and','%3s'%str(C+1),'with a similarity score of',np.around(value,decimals=4))
activeClusters[0,C]=0
### 3. Save the resulting clustering and go back to 1 if the number of clusters >1
mergingClusteringIndices = np.where(clusteringTable[:,k]==C+1)
# update clustering table
clusteringTable[mergingClusteringIndices[0],k]=R+1
# remove binarykey for removed cluster
clustersBKTable[C,:]=np.zeros([1,kbmSize])
clustersCVTable[C,:]=np.nan
# prepare the vector with the indices of the features of thew new cluster and then binarize
segmentsToBinarize = np.where(clusteringTable[:,k]==R+1)[0]
M=[]
for l in np.arange(np.size(segmentsToBinarize,0)):
M = np.append(M,np.arange(int(segmentTable[segmentsToBinarize][:][l,1]),int(segmentTable[segmentsToBinarize][:][l,2])+1))
clustersBKTable[R,:], clustersCVTable[R,:]=binarizeFeatures(kbmSize,Vg[np.array(speechMapping[np.array(M,dtype='int')],dtype='int')-1].T,bitsPerSegmentFactor)
print('done')
return clusteringTable, k
def calcClusters(clustersCVTable,clustersBKTable,activeClusters,clusteringTable,N_init,segmentTable,kbmSize,speechMapping,Vg,bitsPerSegmentFactor):
for i in np.arange(N_init):
if activeClusters[i]==1:
segmentsToBinarize = np.where(clusteringTable==i+1)[0]
M = []
for l in np.arange(np.size(segmentsToBinarize,0)):
M = np.append(M,np.arange(int(segmentTable[segmentsToBinarize][:][l,1]),int(segmentTable[segmentsToBinarize][:][l,2])+1))
clustersBKTable[i], clustersCVTable[i] = binarizeFeatures(kbmSize,Vg[np.array(speechMapping[np.array(M,dtype='int')],dtype='int')-1].T,bitsPerSegmentFactor)
else:
clustersBKTable[i]=np.zeros([1,kbmSize])
clustersCVTable[i]=np.nan
return clustersBKTable, clustersCVTable
def binaryKeySimilarity_cdist(clusteringMetric,bkT1, cvT1, bkT2, cvT2):
from scipy.spatial.distance import cdist
if clusteringMetric == 'cosine':
S = 1 - cdist(cvT1,cvT2,metric=clusteringMetric)
elif clusteringMetric == 'jaccard':
S = 1 - cdist(bkT1,bkT2,metric=clusteringMetric)
else:
print('Clustering metric must be cosine or jaccard')
return S
def getBestClustering(bestClusteringMetric, bkT, cvT, clusteringTable, n, maxNrSpeakers):
from scipy.spatial.distance import cdist
wss = np.zeros([1,n])
overallMean = np.mean(cvT,0)
if bestClusteringMetric == 'cosine':
distances = cdist(np.expand_dims(overallMean,axis=0),cvT,bestClusteringMetric)
elif bestClusteringMetric == 'jaccard':
nBitsTol = np.sum(bkT[0,:])
indices = np.argsort(-overallMean)
overallMean = np.zeros([1,np.size(bkT,1)])
overallMean[0,indices[np.arange(nBitsTol).astype(int)]]=1
distances = cdist(overallMean,bkT,bestClusteringMetric)
distances2 = np.square(distances)
wss[0,n-1] = np.sum(distances2)
for i in np.arange(n-1):
T = clusteringTable[:,i]
clusterIDs = np.unique(T)
variances = np.zeros([np.size(clusterIDs,0),1])
for j in np.arange(np.size(clusterIDs,0)):
clusterIDsIndex = np.where(T==clusterIDs[j])
meanVector = np.mean(cvT[clusterIDsIndex,:],axis=1)
if bestClusteringMetric=='cosine':
distances = cdist(meanVector,cvT[clusterIDsIndex,:][0],bestClusteringMetric)
elif bestClusteringMetric=='jaccard':
indices = np.argsort(-meanVector[0,:])
meanVector = np.zeros([1,np.size(bkT,1)])
meanVector[0,indices[np.arange(nBitsTol).astype(int)]]=1
distances = cdist(meanVector,bkT[clusterIDsIndex,:][0],bestClusteringMetric)
distances2 = np.square(distances)
variances[j] = np.sum(distances2)
wss[0,i]=np.sum(variances)
nPoints = np.size(wss,1)
allCoord = np.vstack((np.arange(1,nPoints+1),wss)).T
firstPoint = allCoord[0,:]
allCoord = allCoord[np.arange(np.where(allCoord[:,1]==np.min(allCoord[:,1]))[0],nPoints),:]
nPoints = np.size(allCoord,0)
lineVec = allCoord[-1,:] - firstPoint
lineVecN = lineVec / np.sqrt(np.sum(np.square(lineVec)))
vecFromFirst = np.subtract(allCoord,firstPoint)
scalarProduct = vecFromFirst*lineVecN
scalarProduct = scalarProduct[:,0]+scalarProduct[:,1]
vecFromFirstParallel = np.expand_dims(scalarProduct,axis=1) * np.expand_dims(lineVecN,0)
vecToLine = vecFromFirst - vecFromFirstParallel
distToLine = np.sqrt(np.sum(np.square(vecToLine),axis=1))
bestClusteringID = allCoord[np.argmax(distToLine)][0]
nrSpeakersPerSolution = np.zeros((clusteringTable.shape[1]))
for k in np.arange(clusteringTable.shape[1]):
nrSpeakersPerSolution[k] = np.size(np.unique(clusteringTable[:,k]))
bestClusteringID = np.maximum(np.where(nrSpeakersPerSolution==np.maximum(np.minimum(maxNrSpeakers,np.max(nrSpeakersPerSolution))-1,1))[0][0],bestClusteringID)
return bestClusteringID
def getSpectralClustering(bestClusteringMetric,clusteringTable,N_init, bkT, cvT, n, sigma, percentile,maxNrSpeakers):
from scipy.ndimage.filters import gaussian_filter
simMatrix = binaryKeySimilarity_cdist(bestClusteringMetric,bkT,cvT,bkT,cvT)
np.fill_diagonal(simMatrix,np.nan)
np.fill_diagonal(simMatrix,np.nanmax(simMatrix,1))
# Similarity matrix smoothed through Gaussian filter
simMatrix_1 = gaussian_filter(simMatrix,sigma)
# Similarity matrix thresholded to the percentile-th components leaving a small amount instead of 0 following Google's implementation
thresholds = np.tile(percentile*np.nanmax(simMatrix,1)/100,(np.size(bkT,0),1)).T
simMatrix_2 = np.copy(simMatrix_1)
mask = simMatrix_2<thresholds
simMatrix_2[np.where(mask==True)] = simMatrix_2[np.where(mask==True)]*0.01
# Similarity matrix made symmetric
simMatrix_3 = np.copy(simMatrix_2)
for k in np.arange(np.size(simMatrix_3,0)):
# for j in np.arange(np.size(simMatrix_3,)):
maximum = np.maximum(simMatrix_3[:,k],simMatrix_3[k,:])
simMatrix_3[:,k] = simMatrix_3[k,:] = maximum
# Similarity matrix diffussion
simMatrix_4 = np.dot(simMatrix_3,simMatrix_3)
# Row-wise max normalization
simMatrix_5 = simMatrix_4 / np.tile(simMatrix_4.max(axis=0),(np.size(simMatrix_4,0),1)).T
# Decomposition in eigenvalues, we don't use eigenvectors for the moment
eigenvalues,eigenvectors = np.linalg.eigh(simMatrix_5)
new_N_init = np.minimum(maxNrSpeakers,N_init)
eigenvalues = np.flip(np.sort(eigenvalues),axis=0)[0:new_N_init]
if new_N_init > 1:
eigengaps = eigenvalues[0:new_N_init-1]/eigenvalues[1:new_N_init]
eigengapssubstract = eigenvalues[0:new_N_init-1] - eigenvalues[1:new_N_init]
else:
eigengaps = eigengapssubstract = eigenvalues
if eigengapssubstract[0] > 340:
kclusters = 1
else:
eigengaps[0]=0
kclusters = np.flip(np.argsort(eigengaps),axis=0)[0]+1
nrElements = np.zeros([N_init,1])
for k in np.arange(N_init):
nrElements[k]=np.size(np.unique(clusteringTable[:,k]),0)
distances = np.abs(nrElements-kclusters)
minidx = np.argmin(distances)
bestClusteringID = minidx
return bestClusteringID
def smooth(a,WSZ):
# a: NumPy 1-D array containing the data to be smoothed
# WSZ: smoothing window size needs, which must be odd number,
# as in the original MATLAB implementation
#From https://stackoverflow.com/a/40443565
out0 = np.convolve(a,np.ones(WSZ,dtype=int),'valid')/WSZ
r = np.arange(1,WSZ-1,2)
start = np.cumsum(a[:WSZ-1])[::2]/r
stop = (np.cumsum(a[:-WSZ:-1])[::2]/r)[::-1]
return np.concatenate((start,out0,stop))
def performResegmentation(data, speechMapping,mask,finalClusteringTable,segmentTable,modelSize,nbIter,smoothWin,numberOfSpeechFeatures):
from sklearn import mixture
np.random.seed(0)
changePoints,segBeg,segEnd,nSegs = unravelMask(mask)
speakerIDs = np.unique(finalClusteringTable)
trainingData = np.empty([2,0])
for i in np.arange(np.size(speakerIDs,0)):
spkID = speakerIDs[i]
speakerFeaturesIndxs = []
idxs = np.where(finalClusteringTable==spkID)[0]
for l in np.arange(np.size(idxs,0)):
speakerFeaturesIndxs = np.append(speakerFeaturesIndxs,np.arange(int(segmentTable[idxs][:][l,1]),int(segmentTable[idxs][:][l,2])+1))
formattedData = np.vstack((np.tile(spkID,(1,np.size(speakerFeaturesIndxs,0))),speakerFeaturesIndxs))
trainingData = np.hstack((trainingData,formattedData))
llkMatrix = np.zeros([np.size(speakerIDs,0),numberOfSpeechFeatures])
for i in np.arange(np.size(speakerIDs,0)):
spkIdxs = np.where(trainingData[0,:]==speakerIDs[i])[0]
spkIdxs = speechMapping[trainingData[1,spkIdxs].astype(int)].astype(int)-1
msize = np.minimum(modelSize,np.size(spkIdxs,0))
w_init = np.ones([msize])/msize
m_init = data[spkIdxs[np.random.randint(np.size(spkIdxs,0), size=(1, msize))[0]],:]
gmm=mixture.GaussianMixture(n_components=msize,covariance_type='diag',weights_init=w_init,means_init=m_init,verbose=0)
gmm.fit(data[spkIdxs,:])
llkSpk = gmm.score_samples(data)
llkSpkSmoothed = np.zeros([1,numberOfSpeechFeatures])
for jx in np.arange(nSegs):
sectionIdx = np.arange(speechMapping[segBeg[jx]]-1,speechMapping[segEnd[jx]]).astype(int)
sectionWin = np.minimum(smoothWin,np.size(sectionIdx))
if sectionWin % 2 ==0:
sectionWin = sectionWin - 1
if sectionWin>=2:
llkSpkSmoothed[0,sectionIdx] = smooth(llkSpk[sectionIdx], sectionWin)
else:
llkSpkSmoothed[0,sectionIdx]=llkSpk[sectionIdx]
llkMatrix[i,:] = llkSpkSmoothed[0].T
segOut = np.argmax(llkMatrix,axis=0)+1
segChangePoints = np.diff(segOut)
changes = np.where(segChangePoints!=0)[0]
relSegEnds = speechMapping[segEnd]
relSegEnds = relSegEnds[0:-1]
changes = np.sort(np.unique(np.hstack((changes,relSegEnds))))
# Create the new segment and clustering tables
currentPoint = 0
finalSegmentTable = np.empty([0,4])
finalClusteringTableResegmentation = np.empty([0,1])
for i in np.arange(np.size(changes,0)):
addedRow = np.hstack((np.tile(np.where(speechMapping==np.maximum(currentPoint,1))[0],(1,2)), np.tile(np.where(speechMapping==np.maximum(1,changes[i].astype(int)))[0],(1,2))))
finalSegmentTable = np.vstack((finalSegmentTable,addedRow[0]))
finalClusteringTableResegmentation = np.vstack((finalClusteringTableResegmentation,segOut[(changes[i]).astype(int)]))
currentPoint = changes[i]+1
addedRow = np.hstack((np.tile(np.where(speechMapping==currentPoint)[0],(1,2)), np.tile(np.where(speechMapping==numberOfSpeechFeatures)[0],(1,2))))
finalSegmentTable = np.vstack((finalSegmentTable,addedRow[0]))
finalClusteringTableResegmentation = np.vstack((finalClusteringTableResegmentation,segOut[(changes[i]+1).astype(int)]))
return finalClusteringTableResegmentation,finalSegmentTable
def getSegmentationFile(format, frameshift,finalSegmentTable, finalClusteringTable, showName, filename, outputPath, outputExt):
numberOfSpeechFeatures = finalSegmentTable[-1,2].astype(int)+1
solutionVector = np.zeros([1,numberOfSpeechFeatures])
for i in np.arange(np.size(finalSegmentTable,0)):
solutionVector[0,np.arange(finalSegmentTable[i,1],finalSegmentTable[i,2]+1).astype(int)]=finalClusteringTable[i]
seg = np.empty([0,3])
solutionDiff = np.diff(solutionVector)[0]
first = 0
for i in np.arange(0,np.size(solutionDiff,0)):
if solutionDiff[i]:
last = i+1
seg1 = (first)*frameshift
seg2 = (last-first)*frameshift
seg3 = solutionVector[0,last-1]
if seg3:
seg = np.vstack((seg,[seg1,seg2,seg3]))
first = i+1
last = np.size(solutionVector,1)
seg1 = (first-1)*frameshift
seg2 = (last-first+1)*frameshift
seg3 = solutionVector[0,last-1]
seg = np.vstack((seg,[seg1,seg2,seg3]))
solution = []
if format=='MDTM':
for i in np.arange(np.size(seg,0)):
solution.append(showName+' 1 '+str(np.around(seg[i,0],decimals=4))+' '+str(np.around(seg[i,1],decimals=4))+' speaker NA unknown speaker'+str(seg[i,2].astype(int)))+'\n'
elif format=='RTTM':
for i in np.arange(np.size(seg,0)):
solution.append('SPEAKER '+showName+' 1 '+str(np.around(seg[i,0],decimals=4))+' '+str(np.around(seg[i,1],decimals=4))+' <NA> <NA> speaker'+str(seg[i,2].astype(int))+' <NA>\n')
else:
print('Output file format must be MDTM or RTTM.')
solution[-1]=solution[-1][0:-1]
outf = open(outputPath+filename+outputExt,"a")
outf.writelines(solution)
outf.write('\n')
outf.close()
class HTKFile:
#Code from: https://github.com/danijel3/PyHTK
""" Class to load binary HTK file.
Details on the format can be found online in HTK Book chapter 5.7.1.
Not everything is implemented 100%, but most features should be supported.
Not implemented:
CRC checking - files can have CRC, but it won't be checked for correctness
VQ - Vector features are not implemented.
"""
data = None
nSamples = 0
nFeatures = 0
sampPeriod = 0
basicKind = None
qualifiers = None
def load(self, filename):
import struct
""" Loads HTK file.
After loading the file you can check the following members:
data (matrix) - data contained in the file
nSamples (int) - number of frames in the file
nFeatures (int) - number if features per frame
sampPeriod (int) - sample period in 100ns units (e.g. fs=16 kHz -> 625)
basicKind (string) - basic feature kind saved in the file
qualifiers (string) - feature options present in the file
"""
with open(filename, "rb") as f:
header = f.read(12)
self.nSamples, self.sampPeriod, sampSize, paramKind = struct.unpack(">iihh", header)
basicParameter = paramKind & 0x3F
if basicParameter is 0:
self.basicKind = "WAVEFORM"
elif basicParameter is 1:
self.basicKind = "LPC"
elif basicParameter is 2:
self.basicKind = "LPREFC"
elif basicParameter is 3:
self.basicKind = "LPCEPSTRA"
elif basicParameter is 4:
self.basicKind = "LPDELCEP"
elif basicParameter is 5:
self.basicKind = "IREFC"
elif basicParameter is 6:
self.basicKind = "MFCC"
elif basicParameter is 7:
self.basicKind = "FBANK"
elif basicParameter is 8:
self.basicKind = "MELSPEC"
elif basicParameter is 9:
self.basicKind = "USER"
elif basicParameter is 10:
self.basicKind = "DISCRETE"
elif basicParameter is 11:
self.basicKind = "PLP"
else:
self.basicKind = "ERROR"
self.qualifiers = []
if (paramKind & 0o100) != 0:
self.qualifiers.append("E")
if (paramKind & 0o200) != 0:
qualifiers.append("N")
if (paramKind & 0o400) != 0:
self.qualifiers.append("D")
if (paramKind & 0o1000) != 0:
self.qualifiers.append("A")
if (paramKind & 0o2000) != 0:
self.qualifiers.append("C")
if (paramKind & 0o4000) != 0:
self.qualifiers.append("Z")
if (paramKind & 0o10000) != 0:
self.qualifiers.append("K")
if (paramKind & 0o20000) != 0:
self.qualifiers.append("0")
if (paramKind & 0o40000) != 0:
self.qualifiers.append("V")
if (paramKind & 0o100000) != 0:
self.qualifiers.append("T")
if "C" in self.qualifiers or "V" in self.qualifiers or self.basicKind is "IREFC" or self.basicKind is "WAVEFORM":
self.nFeatures = sampSize // 2
else:
self.nFeatures = sampSize // 4
if "C" in self.qualifiers:
self.nSamples -= 4
if "V" in self.qualifiers:
raise NotImplementedError("VQ is not implemented")
self.data = []
if self.basicKind is "IREFC" or self.basicKind is "WAVEFORM":
for x in range(self.nSamples):
s = f.read(sampSize)
frame = []
for v in range(self.nFeatures):
val = struct.unpack_from(">h", s, v * 2)[0] / 32767.0
frame.append(val)
self.data.append(frame)
elif "C" in self.qualifiers:
A = []
s = f.read(self.nFeatures * 4)
for x in range(self.nFeatures):
A.append(struct.unpack_from(">f", s, x * 4)[0])
B = []
s = f.read(self.nFeatures * 4)
for x in range(self.nFeatures):
B.append(struct.unpack_from(">f", s, x * 4)[0])
for x in range(self.nSamples):
s = f.read(sampSize)
frame = []
for v in range(self.nFeatures):
frame.append((struct.unpack_from(">h", s, v * 2)[0] + B[v]) / A[v])
self.data.append(frame)
else:
for x in range(self.nSamples):
s = f.read(sampSize)
frame = []
for v in range(self.nFeatures):
val = struct.unpack_from(">f", s, v * 4)
frame.append(val[0])
self.data.append(frame)
if "K" in self.qualifiers:
print("CRC checking not implememented...")
