Skip to content

Commit

Permalink
update the inference model for semi-markov events
Browse files Browse the repository at this point in the history
  • Loading branch information
@charlesneimog
charlesneimog committed Sep 2, 2024
1 parent 337804c commit 910dcac
Show file tree
Hide file tree
Showing 7 changed files with 222 additions and 186 deletions.
131 changes: 76 additions & 55 deletions Sources/OScofo/mdp.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include "mdp.hpp"
#include "log.hpp"
#include "states.hpp"
#include <algorithm>

#include <boost/math/special_functions/bessel.hpp>

Expand All @@ -10,13 +11,14 @@
OScofoMDP::OScofoMDP(float Sr, float WindowSize, float HopSize) {

// Check if the outer map contains the key

m_HopSize = HopSize;
m_WindowSize = WindowSize;
m_Sr = Sr;
m_AccumulationFactor = 0.5;
m_CouplingStrength = 0.5;

m_BlockDur = 1 / m_Sr * m_HopSize;

SetTunning(440);
CreateKappaTable();
}
Expand All @@ -34,27 +36,37 @@ void OScofoMDP::SetScoreStates(States States) {
m_PsiN = 60.0f / m_States[0].BPMExpected;
m_PsiN1 = 60.0f / m_States[0].BPMExpected;
m_LastPsiN = 60.0f / m_States[0].BPMExpected;
m_BeatsAhead = m_States[0].BPMExpected / 60 * 2;
m_BeatsAhead = m_States[0].BPMExpected / 60 * 4;

CreateMarkovMatrix();

int size = m_States.size();

for (int i = 0; i < size; i++) {
m_EventObservations[i].resize(m_MaxHistory + 1, 1);
m_PitchObservations[i].resize(m_MaxHistory + 1, 1);
}
}

// ─────────────────────────────────────
void OScofoMDP::CreateMarkovMatrix() {
int size = m_States.size();
// TODO: REVER

int i = -1;
for (int i = -1; i < size; i++) {
double Dur = 0;
for (int j = i + 1; j < m_States.size(); j++) {
double DurProb = 1 - m_States[j].Duration / m_BeatsAhead;
double DurProb = exp(-5 * (Dur / m_BeatsAhead));
double LastDur = Dur / m_BeatsAhead;
if (i == j - 1) {
m_TransitionsProb[i][j] = 1;
} else {
m_TransitionsProb[i][j] = 0;
}
// printf("%d -> %d: %f\n", i, j, m_TransitionsProb[i][j]);
m_TransitionsProb[i][j] = DurProb;
// if (i - 1 == j) {
// m_TransitionsProb[i][j] = 1;
// } else {
// m_TransitionsProb[i][j] = 0;
// }

Dur += m_States[j].Duration;
}
double Sum = 0;
int innerSize = m_TransitionsProb[i].size();
Expand Down Expand Up @@ -92,9 +104,6 @@ void OScofoMDP::UpdatePitchTemplate() {
double num = std::pow(i - (RootBinFreq * (j + 1)), 2);
double den = 2 * M_PI * m_PitchTemplateSigma * m_PitchTemplateSigma;
m_PitchTemplates[RootBinFreq][i] += amp * std::exp(-(num / den));
if (m_PitchTemplateHigherBin < RootBinFreq) {
m_PitchTemplateHigherBin = RootBinFreq;
}
}
}
}
Expand Down Expand Up @@ -355,9 +364,6 @@ double OScofoMDP::GetPitchSimilarity(State &State, Description &Desc) {
}

for (size_t i = 0; i < m_WindowSize / 2; i++) {
if (i > m_PitchTemplateHigherBin) {
break;
}
double P = PitchTemplate[i] * Desc.MaxAmp;
double Q = Desc.NormSpectralPower[i];
if (P > 0 && Q > 0) {
Expand All @@ -372,63 +378,75 @@ double OScofoMDP::GetPitchSimilarity(State &State, Description &Desc) {
}

// ─────────────────────────────────────
double OScofoMDP::GetSojournTime(State &State, Description &Desc) {
double T = m_Tn + m_TimeInThisEvent;
double OScofoMDP::GetSojournTime(State &State, int u) {
double T = m_Tn + (m_BlockDur * u);
double Duration = State.Duration;
double Sojourn = std::exp(-((T - m_Tn) / (m_PsiN1 * Duration)));
double Sojourn = std::exp(-(T - m_Tn) / (m_PsiN1 * Duration));
return Sojourn;
}

// ─────────────────────────────────────
double OScofoMDP::SemiMarkovState(Description &Desc, int j) {
if (m_T == 0) {
return m_LastDecodedPosition;
}
int OScofoMDP::MaxBlocksInState(State &State) {
int MaxBlocksInState = State.Duration / m_BlockDur;
return MaxBlocksInState + 1;
}

unsigned int t = m_T;
if (t > 20) {
t = 20;
// ─────────────────────────────────────
double OScofoMDP::SemiMarkovState(Description &Desc, int j) {
int t = m_T;
if (m_T > m_MaxHistory) {
t = m_MaxHistory;
}

// Observation probability b_j(x_t)
double KLDiv = GetPitchSimilarity(m_States[j], Desc);
m_PitchKL[j][m_T] = KLDiv;

// Maximum value initialization
double maxAlphaJ = -INFINITY;

// Loop over possible time steps u
for (int u = 1; u <= t; u++) {
// Product term ∏_{v=1}^{u-1} b_j(x_{t-v})
double ProdTerm = 1;
// max
double MaxOverU = 0;
for (int u = 1; u < std::min(t, MaxBlocksInState(m_States[j])); u++) {
double Obs = 1;
for (int v = 1; v < u; v++) {
ProdTerm *= m_PitchKL[j][m_T - v];
Obs *= m_PitchObservations[j][t - v];
}

// Survival function d_j(u) (this should be defined based on your model)
double SojournTime = GetSojournTime(m_States[j], Desc);
// value 2
double Sojourn = GetSojournTime(m_States[j], u);

// Transition probabilities max_{i ≠ j} (p_{ij} a_i(t-u))
double MaxProb = -INFINITY;
// value 3
double Prob = 1;
for (int i = m_CurrentEvent; i < m_StateWindow; i++) {
if (i != j) {
double EventProb = m_TransitionsProb[i][j] * m_PrevObserved[i][m_T - u];
if (MaxProb < EventProb) {
MaxProb = EventProb;
double EventIObs;
if (i == -1) {
EventIObs = m_TransitionsProb[i][j];
} else {
EventIObs = m_EventObservations[i][t - u];
}
Prob = fmax(Prob, m_TransitionsProb[i][j] * EventIObs);
}
}
double AlphaJ = ProdTerm * SojournTime; //* MaxProb;
if (AlphaJ > maxAlphaJ) {
maxAlphaJ = AlphaJ;
}
MaxOverU = fmax(MaxOverU, Obs * Sojourn * Prob);
}

double Pitch = GetPitchSimilarity(m_States[j], Desc);
double AlphaJ = Pitch * MaxOverU;

// Rotate the observations history
if (m_T < m_MaxHistory) {
m_EventObservations[j][m_T] = MaxOverU;
m_PitchObservations[j][m_T] = Pitch;
} else {
std::copy(m_EventObservations[j].begin() + 1, m_EventObservations[j].end(),
m_EventObservations[j].begin());
std::copy(m_PitchObservations[j].begin() + 1, m_PitchObservations[j].end(),
m_PitchObservations[j].begin());
m_EventObservations[j][m_MaxHistory] = MaxOverU;
m_PitchObservations[j][m_MaxHistory] = MaxOverU;
}
return KLDiv * maxAlphaJ;

return AlphaJ;
}

// ─────────────────────────────────────
double OScofoMDP::MarkovState(Description &Desc, int j) {
m_PitchKL[j][m_T] = GetPitchSimilarity(m_States[j], Desc);
// m_PitchKL[j][m_T] = GetPitchSimilarity(m_States[j], Desc);
// m_PrevObserved[j][m_T] = AlphaJ;

return 0;
}
Expand All @@ -440,16 +458,18 @@ double OScofoMDP::GetBestEventIndex(Description &Desc) {

double BestAlphaJ = 0;
for (int j = m_LastDecodedPosition; j < m_StateWindow + 1; j++) {
double AlphaJ;

if (j < 0) {
// TODO: Rethink about first event
j++;
}
double AlphaJ;

if (m_States[j].MarkovType == SEMIMARKOV) {
AlphaJ = SemiMarkovState(Desc, j);
} else if (m_States[j].MarkovType == MARKOV) {
AlphaJ = MarkovState(Desc, j);
}
m_PrevObserved[j][m_T] = AlphaJ;

if (AlphaJ > BestAlphaJ) {
BestAlphaJ = AlphaJ;
Expand All @@ -467,13 +487,14 @@ double OScofoMDP::GetBestEventIndex(Description &Desc) {

// ─────────────────────────────────────
int OScofoMDP::GetEvent(Description &Desc) {
double BlockDur = 1 / m_Sr;
m_TimeInThisEvent += BlockDur * m_HopSize;
m_TimeInThisEvent += m_BlockDur;
m_T += 1;
m_AudioTick += 1;

if (!Desc.PassTreshold || m_CurrentEvent == m_States.size()) {
return m_CurrentEvent;
}

if (m_CurrentEvent == -1) {
m_BPM = m_States[0].BPMExpected;
m_PsiN = 60 / m_BPM;
Expand Down
18 changes: 12 additions & 6 deletions Sources/OScofo/mdp.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,9 @@
#endif

using PitchTemplateArray = std::vector<double>;

using MatrixMap = std::unordered_map<unsigned int, std::unordered_map<unsigned int, double>>;
using Array2DMap = std::unordered_map<unsigned int, std::vector<double>>;

// ╭─────────────────────────────────────╮
// │ Markov Description Process │
Expand Down Expand Up @@ -45,8 +47,9 @@ class OScofoMDP {
double m_Sr;
double m_WindowSize;
double m_HopSize;
int m_AudioTick = 0;
double m_Harmonics = 10;
double m_PitchTemplateHigherBin = 0;
int m_PitchTemplateHigherBin = 0;
double m_dBTreshold = -55;

// Events
Expand All @@ -70,17 +73,21 @@ class OScofoMDP {
double ForwardRecursion(Description &Desc);
double SemiMarkovState(Description &Desc, int j);
double MarkovState(Description &Desc, int j);
double SemiMarkovStateTest(Description &Desc, int j);
double GetSojournTime(State &PossibleState, int u);
int MaxBlocksInState(State &State);

int m_StateWindow = 0;
int m_T = 0;
int m_BlocksInThisEvent = 0;
int m_LastDecodedPosition = -1;
MatrixMap m_Believes;
MatrixMap m_PitchKL;
int m_MaxHistory = 80;

MatrixMap m_TransitionsProb;
MatrixMap m_Observed;
MatrixMap m_PrevObserved;
Array2DMap m_EventObservations;
Array2DMap m_PitchObservations;

double m_BlockDur = 0;
double m_BeatsAhead = 8;
double m_TimeInThisEvent = 0;
double m_LastTn = 0;
Expand Down Expand Up @@ -117,7 +124,6 @@ class OScofoMDP {
std::vector<State> m_States;
double GetReward(States &NextPossibleStates, Description &Desc);
double GetPitchSimilarity(State &NextPossibleState, Description &Desc);
double GetSojournTime(State &PossibleState, Description &Desc);

// Viterbi | SemiMarkov
double ComputeViterbiRecursion(State &JState, Description &Desc);
Expand Down
2 changes: 1 addition & 1 deletion Sources/OScofo/mir.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ class OScofoMIR {
std::vector<double> m_LastPhase2;

// Env
double m_dBTreshold = -80;
double m_dBTreshold = -60;
void GetRMS(const std::vector<double> &In, Description &Desc);

// Onset Detection
Expand Down
1 change: 1 addition & 0 deletions Sources/OScofo/states.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ enum HMMType { SEMIMARKOV = 0, MARKOV };
// ─────────────────────────────────────
struct State {
int Index;
int ScoreEventNumber;
EventType Type; // NOTE, CHORD, PIZZ, SILENCE
HMMType MarkovType; // NOTE, CHORD, PIZZ, SILENCE
std::vector<double> Freqs;
Expand Down
2 changes: 1 addition & 1 deletion Tests/Test0.txt
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
BPM 30
BPM 40
NOTE 82 0.75
NOTE 81 0.125
NOTE 83 0.125
Expand Down
Loading

0 comments on commit 910dcac

Please sign in to comment.