Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Format python files with black #150

Closed
wants to merge 4 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .git-blame-ignore-revs
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
# format files with black
afe57d2e21ad53c69d67d61a30341bb5a1dc735d

10 changes: 10 additions & 0 deletions .github/workflows/black.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
name: Lint

on: [push, pull_request]

jobs:
lint:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4
- uses: psf/black@stable
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ name = "SMSTools"
authors = [
{ name="Music Technology Group, Universitat Pompeu Fabra", email="mtg-info@upf.edu" },
]
version = "0.99.1"
version = "0.90"
description = "Sound analysis/synthesis tools for music applications"
readme = "README.md"
requires-python = ">=3.9"
Expand Down
50 changes: 31 additions & 19 deletions smstools/models/dftModel.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,15 +11,15 @@

def dftModel(x, w, N):
"""
Analysis/synthesis of a signal using the discrete Fourier transform
x: input signal, w: analysis window, N: FFT size
returns y: output signal
"""
Analysis/synthesis of a signal using the discrete Fourier transform
x: input signal, w: analysis window, N: FFT size
returns y: output signal
"""

if not (UF.isPower2(N)): # raise error if N not a power of two
raise ValueError("FFT size (N) is not a power of 2")

if (w.size > N): # raise error if window size bigger than fft size
if w.size > N: # raise error if window size bigger than fft size
raise ValueError("Window size (M) is bigger than FFT size")

if all(x == 0): # if input array is zeros return empty output
Expand All @@ -35,13 +35,17 @@ def dftModel(x, w, N):
fftbuffer[-hM2:] = xw[:hM2]
X = fft(fftbuffer) # compute FFT
absX = abs(X[:hN]) # compute ansolute value of positive side
absX[absX < np.finfo(float).eps] = np.finfo(float).eps # if zeros add epsilon to handle log
absX[absX < np.finfo(float).eps] = np.finfo(
float
).eps # if zeros add epsilon to handle log
mX = 20 * np.log10(absX) # magnitude spectrum of positive frequencies in dB
pX = np.unwrap(np.angle(X[:hN])) # unwrapped phase spectrum of positive frequencies
# -----synthesis-----
Y = np.zeros(N, dtype=complex) # clean output spectrum
Y[:hN] = 10 ** (mX / 20) * np.exp(1j * pX) # generate positive frequencies
Y[hN:] = 10 ** (mX[-2:0:-1] / 20) * np.exp(-1j * pX[-2:0:-1]) # generate negative frequencies
Y[hN:] = 10 ** (mX[-2:0:-1] / 20) * np.exp(
-1j * pX[-2:0:-1]
) # generate negative frequencies
fftbuffer = np.real(ifft(Y)) # compute inverse FFT
y[:hM2] = fftbuffer[-hM2:] # undo zero-phase window
y[hM2:] = fftbuffer[:hM1]
Expand All @@ -50,10 +54,10 @@ def dftModel(x, w, N):

def dftAnal(x, w, N):
"""
Analysis of a signal using the discrete Fourier transform
x: input signal, w: analysis window, N: FFT size
returns mX, pX: magnitude and phase spectrum
"""
Analysis of a signal using the discrete Fourier transform
x: input signal, w: analysis window, N: FFT size
returns mX, pX: magnitude and phase spectrum
"""

if not (UF.isPower2(N)): # raise error if N not a power of two
raise ValueError("FFT size (N) is not a power of 2")
Expand All @@ -71,20 +75,26 @@ def dftAnal(x, w, N):
fftbuffer[-hM2:] = xw[:hM2]
X = fft(fftbuffer) # compute FFT
absX = abs(X[:hN]) # compute ansolute value of positive side
absX[absX < np.finfo(float).eps] = np.finfo(float).eps # if zeros add epsilon to handle log
absX[absX < np.finfo(float).eps] = np.finfo(
float
).eps # if zeros add epsilon to handle log
mX = 20 * np.log10(absX) # magnitude spectrum of positive frequencies in dB
X[:hN].real[np.abs(X[:hN].real) < tol] = 0.0 # for phase calculation set to 0 the small values
X[:hN].imag[np.abs(X[:hN].imag) < tol] = 0.0 # for phase calculation set to 0 the small values
X[:hN].real[
np.abs(X[:hN].real) < tol
] = 0.0 # for phase calculation set to 0 the small values
X[:hN].imag[
np.abs(X[:hN].imag) < tol
] = 0.0 # for phase calculation set to 0 the small values
pX = np.unwrap(np.angle(X[:hN])) # unwrapped phase spectrum of positive frequencies
return mX, pX


def dftSynth(mX, pX, M):
"""
Synthesis of a signal using the discrete Fourier transform
mX: magnitude spectrum, pX: phase spectrum, M: window size
returns y: output signal
"""
Synthesis of a signal using the discrete Fourier transform
mX: magnitude spectrum, pX: phase spectrum, M: window size
returns y: output signal
"""

hN = mX.size # size of positive spectrum, it includes sample 0
N = (hN - 1) * 2 # FFT size
Expand All @@ -96,7 +106,9 @@ def dftSynth(mX, pX, M):
y = np.zeros(M) # initialize output array
Y = np.zeros(N, dtype=complex) # clean output spectrum
Y[:hN] = 10 ** (mX / 20) * np.exp(1j * pX) # generate positive frequencies
Y[hN:] = 10 ** (mX[-2:0:-1] / 20) * np.exp(-1j * pX[-2:0:-1]) # generate negative frequencies
Y[hN:] = 10 ** (mX[-2:0:-1] / 20) * np.exp(
-1j * pX[-2:0:-1]
) # generate negative frequencies
fftbuffer = np.real(ifft(Y)) # compute inverse FFT
y[:hM2] = fftbuffer[-hM2:] # undo zero-phase window
y[hM2:] = fftbuffer[:hM1]
Expand Down
156 changes: 92 additions & 64 deletions smstools/models/harmonicModel.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,26 +12,28 @@

def f0Detection(x, fs, w, N, H, t, minf0, maxf0, f0et):
"""
Fundamental frequency detection of a sound using twm algorithm
x: input sound; fs: sampling rate; w: analysis window;
N: FFT size; t: threshold in negative dB,
minf0: minimum f0 frequency in Hz, maxf0: maximim f0 frequency in Hz,
f0et: error threshold in the f0 detection (ex: 5),
returns f0: fundamental frequency
"""
if (minf0 < 0): # raise exception if minf0 is smaller than 0
Fundamental frequency detection of a sound using twm algorithm
x: input sound; fs: sampling rate; w: analysis window;
N: FFT size; t: threshold in negative dB,
minf0: minimum f0 frequency in Hz, maxf0: maximim f0 frequency in Hz,
f0et: error threshold in the f0 detection (ex: 5),
returns f0: fundamental frequency
"""
if minf0 < 0: # raise exception if minf0 is smaller than 0
raise ValueError("Minumum fundamental frequency (minf0) smaller than 0")

if (maxf0 >= 10000): # raise exception if maxf0 is bigger than fs/2
if maxf0 >= 10000: # raise exception if maxf0 is bigger than fs/2
raise ValueError("Maximum fundamental frequency (maxf0) bigger than 10000Hz")

if (H <= 0): # raise error if hop size 0 or negative
if H <= 0: # raise error if hop size 0 or negative
raise ValueError("Hop size (H) smaller or equal to 0")

hN = N // 2 # size of positive spectrum
hM1 = int(math.floor((w.size + 1) / 2)) # half analysis window size by rounding
hM2 = int(math.floor(w.size / 2)) # half analysis window size by floor
x = np.append(np.zeros(hM2), x) # add zeros at beginning to center first window at sample 0
x = np.append(
np.zeros(hM2), x
) # add zeros at beginning to center first window at sample 0
x = np.append(x, np.zeros(hM1)) # add zeros at the end to analyze last sample
pin = hM1 # init sound pointer in middle of anal window
pend = x.size - hM1 # last sample to start a frame
Expand All @@ -41,14 +43,15 @@ def f0Detection(x, fs, w, N, H, t, minf0, maxf0, f0et):
f0t = 0 # initialize f0 track
f0stable = 0 # initialize f0 stable
while pin < pend:
x1 = x[pin - hM1:pin + hM2] # select frame
x1 = x[pin - hM1 : pin + hM2] # select frame
mX, pX = DFT.dftAnal(x1, w, N) # compute dft
ploc = UF.peakDetection(mX, t) # detect peak locations
iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc) # refine peak values
ipfreq = fs * iploc / N # convert locations to Hez
f0t = UF.f0Twm(ipfreq, ipmag, f0et, minf0, maxf0, f0stable) # find f0
if ((f0stable == 0) & (f0t > 0)) \
or ((f0stable > 0) & (np.abs(f0stable - f0t) < f0stable / 5.0)):
if ((f0stable == 0) & (f0t > 0)) or (
(f0stable > 0) & (np.abs(f0stable - f0t) < f0stable / 5.0)
):
f0stable = f0t # consider a stable f0 if it is close to the previous one
else:
f0stable = 0
Expand All @@ -59,30 +62,36 @@ def f0Detection(x, fs, w, N, H, t, minf0, maxf0, f0et):

def harmonicDetection(pfreq, pmag, pphase, f0, nH, hfreqp, fs, harmDevSlope=0.01):
"""
Detection of the harmonics of a frame from a set of spectral peaks using f0
to the ideal harmonic series built on top of a fundamental frequency
pfreq, pmag, pphase: peak frequencies, magnitudes and phases
f0: fundamental frequency, nH: number of harmonics,
hfreqp: harmonic frequencies of previous frame,
fs: sampling rate; harmDevSlope: slope of change of the deviation allowed to perfect harmonic
returns hfreq, hmag, hphase: harmonic frequencies, magnitudes, phases
"""

if (f0 <= 0): # if no f0 return no harmonics
Detection of the harmonics of a frame from a set of spectral peaks using f0
to the ideal harmonic series built on top of a fundamental frequency
pfreq, pmag, pphase: peak frequencies, magnitudes and phases
f0: fundamental frequency, nH: number of harmonics,
hfreqp: harmonic frequencies of previous frame,
fs: sampling rate; harmDevSlope: slope of change of the deviation allowed to perfect harmonic
returns hfreq, hmag, hphase: harmonic frequencies, magnitudes, phases
"""

if f0 <= 0: # if no f0 return no harmonics
return np.zeros(nH), np.zeros(nH), np.zeros(nH)
hfreq = np.zeros(nH) # initialize harmonic frequencies
hmag = np.zeros(nH) - 100 # initialize harmonic magnitudes
hphase = np.zeros(nH) # initialize harmonic phases
hf = f0 * np.arange(1, nH + 1) # initialize harmonic frequencies
hi = 0 # initialize harmonic index
if len(hfreqp) == 0: # if no incomming harmonic tracks initialize to harmonic series
if (
len(hfreqp) == 0
): # if no incomming harmonic tracks initialize to harmonic series
hfreqp = hf
while (f0 > 0) and (hi < nH) and (hf[hi] < fs / 2): # find harmonic peaks
pei = np.argmin(abs(pfreq - hf[hi])) # closest peak
dev1 = abs(pfreq[pei] - hf[hi]) # deviation from perfect harmonic
dev2 = (abs(pfreq[pei] - hfreqp[hi]) if hfreqp[hi] > 0 else fs) # deviation from previous frame
dev2 = (
abs(pfreq[pei] - hfreqp[hi]) if hfreqp[hi] > 0 else fs
) # deviation from previous frame
threshold = f0 / 3 + harmDevSlope * pfreq[pei]
if ((dev1 < threshold) or (dev2 < threshold)): # accept peak if deviation is small
if (dev1 < threshold) or (
dev2 < threshold
): # accept peak if deviation is small
hfreq[hi] = pfreq[pei] # harmonic frequencies
hmag[hi] = pmag[pei] # harmonic magnitudes
hphase[hi] = pphase[pei] # harmonic phases
Expand All @@ -92,19 +101,21 @@ def harmonicDetection(pfreq, pmag, pphase, f0, nH, hfreqp, fs, harmDevSlope=0.01

def harmonicModel(x, fs, w, N, t, nH, minf0, maxf0, f0et):
"""
Analysis/synthesis of a sound using the sinusoidal harmonic model
x: input sound, fs: sampling rate, w: analysis window,
N: FFT size (minimum 512), t: threshold in negative dB,
nH: maximum number of harmonics, minf0: minimum f0 frequency in Hz,
maxf0: maximim f0 frequency in Hz,
f0et: error threshold in the f0 detection (ex: 5),
returns y: output array sound
"""
Analysis/synthesis of a sound using the sinusoidal harmonic model
x: input sound, fs: sampling rate, w: analysis window,
N: FFT size (minimum 512), t: threshold in negative dB,
nH: maximum number of harmonics, minf0: minimum f0 frequency in Hz,
maxf0: maximim f0 frequency in Hz,
f0et: error threshold in the f0 detection (ex: 5),
returns y: output array sound
"""

hN = N // 2 # size of positive spectrum
hM1 = int(math.floor((w.size + 1) / 2)) # half analysis window size by rounding
hM2 = int(math.floor(w.size / 2)) # half analysis window size by floor
x = np.append(np.zeros(hM2), x) # add zeros at beginning to center first window at sample 0
x = np.append(
np.zeros(hM2), x
) # add zeros at beginning to center first window at sample 0
x = np.append(x, np.zeros(hM1)) # add zeros at the end to analyze last sample
Ns = 512 # FFT size for synthesis (even)
H = Ns // 4 # Hop size used for analysis and synthesis
Expand All @@ -117,57 +128,70 @@ def harmonicModel(x, fs, w, N, t, nH, minf0, maxf0, f0et):
w = w / sum(w) # normalize analysis window
sw = np.zeros(Ns) # initialize synthesis window
ow = triang(2 * H) # overlapping window
sw[hNs - H:hNs + H] = ow
sw[hNs - H : hNs + H] = ow
bh = blackmanharris(Ns) # synthesis window
bh = bh / sum(bh) # normalize synthesis window
sw[hNs - H:hNs + H] = sw[hNs - H:hNs + H] / bh[hNs - H:hNs + H] # window for overlap-add
sw[hNs - H : hNs + H] = (
sw[hNs - H : hNs + H] / bh[hNs - H : hNs + H]
) # window for overlap-add
hfreqp = []
f0t = 0
f0stable = 0
while pin < pend:
# -----analysis-----
x1 = x[pin - hM1:pin + hM2] # select frame
x1 = x[pin - hM1 : pin + hM2] # select frame
mX, pX = DFT.dftAnal(x1, w, N) # compute dft
ploc = UF.peakDetection(mX, t) # detect peak locations
iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc) # refine peak values
ipfreq = fs * iploc / N
f0t = UF.f0Twm(ipfreq, ipmag, f0et, minf0, maxf0, f0stable) # find f0
if ((f0stable == 0) & (f0t > 0)) \
or ((f0stable > 0) & (np.abs(f0stable - f0t) < f0stable / 5.0)):
if ((f0stable == 0) & (f0t > 0)) or (
(f0stable > 0) & (np.abs(f0stable - f0t) < f0stable / 5.0)
):
f0stable = f0t # consider a stable f0 if it is close to the previous one
else:
f0stable = 0
hfreq, hmag, hphase = harmonicDetection(ipfreq, ipmag, ipphase, f0t, nH, hfreqp, fs) # find harmonics
hfreq, hmag, hphase = harmonicDetection(
ipfreq, ipmag, ipphase, f0t, nH, hfreqp, fs
) # find harmonics
hfreqp = hfreq
# -----synthesis-----
Yh = UF.genSpecSines(hfreq, hmag, hphase, Ns, fs) # generate spec sines
fftbuffer = np.real(ifft(Yh)) # inverse FFT
yh[:hNs - 1] = fftbuffer[hNs + 1:] # undo zero-phase window
yh[hNs - 1:] = fftbuffer[:hNs + 1]
y[pin - hNs:pin + hNs] += sw * yh # overlap-add
yh[: hNs - 1] = fftbuffer[hNs + 1 :] # undo zero-phase window
yh[hNs - 1 :] = fftbuffer[: hNs + 1]
y[pin - hNs : pin + hNs] += sw * yh # overlap-add
pin += H # advance sound pointer
y = np.delete(y, range(hM2)) # delete half of first window which was added in stftAnal
y = np.delete(y, range(y.size - hM1, y.size)) # add zeros at the end to analyze last sample
y = np.delete(
y, range(hM2)
) # delete half of first window which was added in stftAnal
y = np.delete(
y, range(y.size - hM1, y.size)
) # add zeros at the end to analyze last sample
return y


def harmonicModelAnal(x, fs, w, N, H, t, nH, minf0, maxf0, f0et, harmDevSlope=0.01, minSineDur=.02):
def harmonicModelAnal(
x, fs, w, N, H, t, nH, minf0, maxf0, f0et, harmDevSlope=0.01, minSineDur=0.02
):
"""
Analysis of a sound using the sinusoidal harmonic model
x: input sound; fs: sampling rate, w: analysis window; N: FFT size (minimum 512); t: threshold in negative dB,
nH: maximum number of harmonics; minf0: minimum f0 frequency in Hz,
maxf0: maximim f0 frequency in Hz; f0et: error threshold in the f0 detection (ex: 5),
harmDevSlope: slope of harmonic deviation; minSineDur: minimum length of harmonics
returns xhfreq, xhmag, xhphase: harmonic frequencies, magnitudes and phases
"""
Analysis of a sound using the sinusoidal harmonic model
x: input sound; fs: sampling rate, w: analysis window; N: FFT size (minimum 512); t: threshold in negative dB,
nH: maximum number of harmonics; minf0: minimum f0 frequency in Hz,
maxf0: maximim f0 frequency in Hz; f0et: error threshold in the f0 detection (ex: 5),
harmDevSlope: slope of harmonic deviation; minSineDur: minimum length of harmonics
returns xhfreq, xhmag, xhphase: harmonic frequencies, magnitudes and phases
"""

if (minSineDur < 0): # raise exception if minSineDur is smaller than 0

if minSineDur < 0: # raise exception if minSineDur is smaller than 0
raise ValueError("Minimum duration of sine tracks smaller than 0")

hN = N // 2 # size of positive spectrum
hM1 = int(math.floor((w.size + 1) / 2)) # half analysis window size by rounding
hM2 = int(math.floor(w.size / 2)) # half analysis window size by floor
x = np.append(np.zeros(hM2), x) # add zeros at beginning to center first window at sample 0
x = np.append(
np.zeros(hM2), x
) # add zeros at beginning to center first window at sample 0
x = np.append(x, np.zeros(hM2)) # add zeros at the end to analyze last sample
pin = hM1 # init sound pointer in middle of anal window
pend = x.size - hM1 # last sample to start a frame
Expand All @@ -177,19 +201,21 @@ def harmonicModelAnal(x, fs, w, N, H, t, nH, minf0, maxf0, f0et, harmDevSlope=0.
f0t = 0 # initialize f0 track
f0stable = 0 # initialize f0 stable
while pin <= pend:
x1 = x[pin - hM1:pin + hM2] # select frame
x1 = x[pin - hM1 : pin + hM2] # select frame
mX, pX = DFT.dftAnal(x1, w, N) # compute dft
ploc = UF.peakDetection(mX, t) # detect peak locations
iploc, ipmag, ipphase = UF.peakInterp(mX, pX, ploc) # refine peak values
ipfreq = fs * iploc / N # convert locations to Hz
f0t = UF.f0Twm(ipfreq, ipmag, f0et, minf0, maxf0, f0stable) # find f0
if ((f0stable == 0) & (f0t > 0)) \
or ((f0stable > 0) & (np.abs(f0stable - f0t) < f0stable / 5.0)):
if ((f0stable == 0) & (f0t > 0)) or (
(f0stable > 0) & (np.abs(f0stable - f0t) < f0stable / 5.0)
):
f0stable = f0t # consider a stable f0 if it is close to the previous one
else:
f0stable = 0
hfreq, hmag, hphase = harmonicDetection(ipfreq, ipmag, ipphase, f0t, nH, hfreqp, fs,
harmDevSlope) # find harmonics
hfreq, hmag, hphase = harmonicDetection(
ipfreq, ipmag, ipphase, f0t, nH, hfreqp, fs, harmDevSlope
) # find harmonics
hfreqp = hfreq
if pin == hM1: # first frame
xhfreq = np.array([hfreq])
Expand All @@ -200,5 +226,7 @@ def harmonicModelAnal(x, fs, w, N, H, t, nH, minf0, maxf0, f0et, harmDevSlope=0.
xhmag = np.vstack((xhmag, np.array([hmag])))
xhphase = np.vstack((xhphase, np.array([hphase])))
pin += H # advance sound pointer
xhfreq = SM.cleaningSineTracks(xhfreq, round(fs * minSineDur / H)) # delete tracks shorter than minSineDur
xhfreq = SM.cleaningSineTracks(
xhfreq, round(fs * minSineDur / H)
) # delete tracks shorter than minSineDur
return xhfreq, xhmag, xhphase
Loading
Loading