Spaces:
Running
Running
File size: 4,145 Bytes
8c92a11 |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 |
# Copyright (c) 2023 Amphion.
#
# This source code is licensed under the MIT license found in the
# LICENSE file in the root directory of this source tree.
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)
) # get periodogram, parametrized like in matlab
fundBin = np.argmax(
ps
) # estimate fundamental at maximum amplitude, get the bin number
fundIndizes = getIndizesAroundPeak(
ps, fundBin
) # get bin numbers around fundamental peak
fundFrequency = faxis[fundBin] # frequency of fundamental
nHarmonics = 18
harmonicFs = getHarmonics(
fundFrequency, sr, nHarmonics=nHarmonics, aliased=True
) # get harmonic frequencies
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]]) # get power of fundamental
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)
|