|
|
|
|
|
|
|
|
|
|
|
import numpy as np |
|
import scipy.signal as sig |
|
import copy |
|
import librosa |
|
|
|
|
|
def bandpower(ps, mode="time"): |
|
""" |
|
estimate bandpower, see https://de.mathworks.com/help/signal/ref/bandpower.html |
|
""" |
|
if mode == "time": |
|
x = ps |
|
l2norm = np.linalg.norm(x) ** 2.0 / len(x) |
|
return l2norm |
|
elif mode == "psd": |
|
return sum(ps) |
|
|
|
|
|
def getIndizesAroundPeak(arr, peakIndex, searchWidth=1000): |
|
peakBins = [] |
|
magMax = arr[peakIndex] |
|
curVal = magMax |
|
for i in range(searchWidth): |
|
newBin = peakIndex + i |
|
if newBin >= len(arr): |
|
break |
|
newVal = arr[newBin] |
|
if newVal > curVal: |
|
break |
|
else: |
|
peakBins.append(int(newBin)) |
|
curVal = newVal |
|
curVal = magMax |
|
for i in range(searchWidth): |
|
newBin = peakIndex - i |
|
if newBin < 0: |
|
break |
|
newVal = arr[newBin] |
|
if newVal > curVal: |
|
break |
|
else: |
|
peakBins.append(int(newBin)) |
|
curVal = newVal |
|
return np.array(list(set(peakBins))) |
|
|
|
|
|
def freqToBin(fAxis, Freq): |
|
return np.argmin(abs(fAxis - Freq)) |
|
|
|
|
|
def getPeakInArea(psd, faxis, estimation, searchWidthHz=10): |
|
""" |
|
returns bin and frequency of the maximum in an area |
|
""" |
|
binLow = freqToBin(faxis, estimation - searchWidthHz) |
|
binHi = freqToBin(faxis, estimation + searchWidthHz) |
|
peakbin = binLow + np.argmax(psd[binLow : binHi + 1]) |
|
return peakbin, faxis[peakbin] |
|
|
|
|
|
def getHarmonics(fund, sr, nHarmonics=6, aliased=False): |
|
harmonicMultipliers = np.arange(2, nHarmonics + 2) |
|
harmonicFs = fund * harmonicMultipliers |
|
if not aliased: |
|
harmonicFs[harmonicFs > sr / 2] = -1 |
|
harmonicFs = np.delete(harmonicFs, harmonicFs == -1) |
|
else: |
|
nyqZone = np.floor(harmonicFs / (sr / 2)) |
|
oddEvenNyq = nyqZone % 2 |
|
harmonicFs = np.mod(harmonicFs, sr / 2) |
|
harmonicFs[oddEvenNyq == 1] = (sr / 2) - harmonicFs[oddEvenNyq == 1] |
|
return harmonicFs |
|
|
|
|
|
def extract_snr(audio, sr=None): |
|
"""Extract Signal-to-Noise Ratio for a given audio.""" |
|
if sr != None: |
|
audio, _ = librosa.load(audio, sr=sr) |
|
else: |
|
audio, sr = librosa.load(audio, sr=sr) |
|
faxis, ps = sig.periodogram( |
|
audio, fs=sr, window=("kaiser", 38) |
|
) |
|
fundBin = np.argmax( |
|
ps |
|
) |
|
fundIndizes = getIndizesAroundPeak( |
|
ps, fundBin |
|
) |
|
fundFrequency = faxis[fundBin] |
|
|
|
nHarmonics = 18 |
|
harmonicFs = getHarmonics( |
|
fundFrequency, sr, nHarmonics=nHarmonics, aliased=True |
|
) |
|
|
|
harmonicBorders = np.zeros([2, nHarmonics], dtype=np.int16).T |
|
fullHarmonicBins = np.array([], dtype=np.int16) |
|
fullHarmonicBinList = [] |
|
harmPeakFreqs = [] |
|
harmPeaks = [] |
|
for i, harmonic in enumerate(harmonicFs): |
|
searcharea = 0.1 * fundFrequency |
|
estimation = harmonic |
|
|
|
binNum, freq = getPeakInArea(ps, faxis, estimation, searcharea) |
|
harmPeakFreqs.append(freq) |
|
harmPeaks.append(ps[binNum]) |
|
allBins = getIndizesAroundPeak(ps, binNum, searchWidth=1000) |
|
fullHarmonicBins = np.append(fullHarmonicBins, allBins) |
|
fullHarmonicBinList.append(allBins) |
|
harmonicBorders[i, :] = [allBins[0], allBins[-1]] |
|
|
|
fundIndizes.sort() |
|
pFund = bandpower(ps[fundIndizes[0] : fundIndizes[-1]]) |
|
|
|
noisePrepared = copy.copy(ps) |
|
noisePrepared[fundIndizes] = 0 |
|
noisePrepared[fullHarmonicBins] = 0 |
|
noiseMean = np.median(noisePrepared[noisePrepared != 0]) |
|
noisePrepared[fundIndizes] = noiseMean |
|
noisePrepared[fullHarmonicBins] = noiseMean |
|
|
|
noisePower = bandpower(noisePrepared) |
|
|
|
r = 10 * np.log10(pFund / noisePower) |
|
|
|
return r, 10 * np.log10(noisePower) |
|
|