% function [delta, psi, qt] = viterbiDecodePCG_Springer(observation_sequence, pi_vector, b_matrix, total_obs_distribution, heartrate, systolic_time, Fs, figures) % % This function calculates the delta, psi and qt matrices associated with % the Viterbi decoding algorithm from: % L. R. Rabiner, "A tutorial on hidden Markov models and selected % applications in speech recognition," Proc. IEEE, vol. 77, no. 2, pp. % 257-286, Feb. 1989. % using equations 32a - 35, and equations 68 - 69 to include duration % dependancy of the states. % % This decoding is performed after the observation probabilities have been % derived from the logistic regression model of Springer et al: % D. Springer et al., "Logistic Regression-HSMM-based Heart Sound % Segmentation," IEEE Trans. Biomed. Eng., In Press, 2015. % % Further, this function is extended to allow the duration distributions to extend % past the beginning and end of the sequence. Without this, the label % sequence has to start and stop with an "entire" state duration being % fulfilled. This extension takes away that requirement, by allowing the % duration distributions to extend past the beginning and end, but only % considering the observations within the sequence for emission probability % estimation. More detail can be found in the publication by Springer et % al., mentioned above. % %% Inputs: % observation_sequence: The observed features % pi_vector: the array of initial state probabilities, dervived from % "trainSpringerSegmentationAlgorithm". % b_matrix: the observation probabilities, dervived from % "trainSpringerSegmentationAlgorithm". % heartrate: the heart rate of the PCG, extracted using % "getHeartRateSchmidt" % systolic_time: the duration of systole, extracted using % "getHeartRateSchmidt" % Fs: the sampling frequency of the observation_sequence % figures: optional boolean variable to show figures % %% Outputs: % logistic_regression_B_matrix: % pi_vector: % total_obs_distribution: % As Springer et al's algorithm is a duration dependant HMM, there is no % need to calculate the A_matrix, as the transition between states is only % dependant on the state durations. % %% Copyright (C) 2016 David Springer % dave.springer@gmail.com % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by % the Free Software Foundation, either version 3 of the License, or % any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU General Public License for more details. % % You should have received a copy of the GNU General Public License % along with this program. If not, see . function [delta, psi, qt] = viterbiDecodePCG_Springer(observation_sequence, pi_vector, b_matrix, total_obs_distribution, heartrate, systolic_time, Fs,figures) if nargin < 8 figures = false; end %% Preliminary springer_options = default_Springer_HSMM_options; T = length(observation_sequence); N = 4; % Number of states % Setting the maximum duration of a single state. This is set to an entire % heart cycle: max_duration_D = round((1*(60/heartrate))*Fs); %Initialising the variables that are needed to find the optimal state path along %the observation sequence. %delta_t(j), as defined on page 264 of Rabiner, is the best score (highest %probability) along a single path, at time t, which accounts for the first %t observations and ends in State s_j. In this case, the length of the %matrix is extended by max_duration_D samples, in order to allow the use %of the extended Viterbi algortithm: delta = ones(T+ max_duration_D-1,N)*-inf; %The argument that maximises the transition between states (this is %basically the previous state that had the highest transition probability %to the current state) is tracked using the psi variable. psi = zeros(T+ max_duration_D-1,N); %An additional variable, that is not included on page 264 or Rabiner, is %the state duration that maximises the delta variable. This is essential %for the duration dependant HMM. psi_duration =zeros(T + max_duration_D-1,N); %% Setting up observation probs observation_probs = zeros(T,N); for n = 1:N %MLR gives P(state|obs) %Therefore, need Bayes to get P(o|state) %P(o|state) = P(state|obs) * P(obs) / P(states) %Where p(obs) is derived from a MVN distribution from all %obserbations, and p(states) is taken from the pi_vector: pihat = mnrval(cell2mat(b_matrix(n)),observation_sequence(:,:)); for t = 1:T Po_correction = mvnpdf(observation_sequence(t,:),cell2mat(total_obs_distribution(1)),cell2mat(total_obs_distribution(2))); %When saving the coefficients from the logistic %regression, it orders them P(class 1) then P(class 2). When %training, I label the classes as 0 and 1, so the %correct probability would be pihat(2). observation_probs(t,n) = (pihat(t,2)*Po_correction)/pi_vector(n); end end %% Setting up state duration probabilities, using Gaussian distributions: [d_distributions, max_S1, min_S1, max_S2, min_S2, max_systole, min_systole, max_diastole, min_diastole] = get_duration_distributions(heartrate,systolic_time); duration_probs = zeros(N,3*Fs); duration_sum = zeros(N,1); for state_j = 1:N for d = 1:max_duration_D if(state_j == 1) duration_probs(state_j,d) = mvnpdf(d,cell2mat(d_distributions(state_j,1)),cell2mat(d_distributions(state_j,2))); if(d < min_S1 || d > max_S1) duration_probs(state_j,d)= realmin; end elseif(state_j==3) duration_probs(state_j,d) = mvnpdf(d,cell2mat(d_distributions(state_j,1)),cell2mat(d_distributions(state_j,2))); if(d < min_S2 || d > max_S2) duration_probs(state_j,d)= realmin; end elseif(state_j==2) duration_probs(state_j,d) = mvnpdf(d,cell2mat(d_distributions(state_j,1)),cell2mat(d_distributions(state_j,2))); if(d < min_systole|| d > max_systole) duration_probs(state_j,d)= realmin; end elseif (state_j==4) duration_probs(state_j,d) = mvnpdf(d,cell2mat(d_distributions(state_j,1)),cell2mat(d_distributions(state_j,2))); if(d < min_diastole ||d > max_diastole) duration_probs(state_j,d)= realmin; end end end duration_sum(state_j) = sum(duration_probs(state_j,:)); end if(length(duration_probs)>3*Fs) duration_probs(:,(3*Fs+1):end) = []; end if(figures) figure('Name', 'Duration probabilities'); plot(duration_probs(1,:)./ duration_sum(1),'Linewidth',2); hold on; plot(duration_probs(2,:)./ duration_sum(2),'r','Linewidth',2); hold on; plot(duration_probs(3,:)./ duration_sum(3),'g','Linewidth',2); hold on; plot(duration_probs(4,:)./ duration_sum(4),'k','Linewidth',2); hold on; legend('S1 Duration','Systolic Duration','S2 Duration','Diastolic Duration'); pause(); end %% Perform the actual Viterbi Recursion: qt = zeros(1,length(delta)); %% Initialisation Step %Equation 32a and 69a, but leave out the probability of being in %state i for only 1 sample, as the state could have started before time t = %0. delta(1,:) = log(pi_vector) + log(observation_probs(1,:)); %first value is the probability of intially being in each state * probability of observation 1 coming from each state %Equation 32b psi(1,:) = -1; % The state duration probabilities are now used. %Change the a_matrix to have zeros along the diagonal, therefore, only %relying on the duration probabilities and observation probabilities to %influence change in states: %This would only be valid in sequences where the transition between states %follows a distinct order. a_matrix = [0,1,0,0;0 0 1 0; 0 0 0 1;1 0 0 0]; %% Run the core Viterbi algorith if(springer_options.use_mex) %% Run Mex code % Ensure you have run the mex viterbi_PhysChallenge.c code on the % native machine before running this: [delta, psi, psi_duration] = viterbi_Springer(N,T,a_matrix,max_duration_D,delta,observation_probs,duration_probs,psi, duration_sum); else %% Recursion %% The Extended Viterbi algorithm: %Equations 33a and 33b and 69a, b, c etc: %again, ommitting the p(d), as state could have started before t = 1 % This implementation extends the standard implementation of the % duration-dependant Viterbi algorithm by allowing the durations to % extend beyond the start and end of the time series, thereby allowing % states to "start" and "stop" outside of the recorded signal. This % addresses the issue of partial states at the beginning and end of the % signal being labelled as the incorrect state. For instance, a % short-duration diastole at the beginning of a signal looks a lot like % systole, and can lead to labelling errors. % t spans input 2 to T + max_duration_D: for t = 2:T+ max_duration_D-1 for j = 1:N for d = 1:1:max_duration_D %The start of the analysis window, which is the current time %step, minus d (the time horizon we are currently looking back), %plus 1. The analysis window can be seen to be starting one %step back each time the variable d is increased. % This is clamped to 1 if extending past the start of the % record, and T-1 is extending past the end of the record: start_t = t - d; if(start_t<1) start_t = 1; end if(start_t > T-1) start_t = T-1; end %The end of the analysis window, which is the current time %step, unless the time has gone past T, the end of the record, in %which case it is truncated to T. This allows the analysis %window to extend past the end of the record, so that the %timing durations of the states do not have to "end" at the end %of the record. end_t = t; if(t>T) end_t = T; end %Find the max_delta and index of that from the previous step %and the transition to the current step: %This is the first half of the expression of equation 33a from %Rabiner: [max_delta, max_index] = max(delta(start_t,:)+log(a_matrix(:,j))'); %Find the normalised probabilities of the observations over the %analysis window: probs = prod(observation_probs(start_t:end_t,j)); %Find the normalised probabilities of the observations at only %the time point at the start of the time window: if(probs ==0) probs = realmin; end emission_probs = log(probs); %Keep a running total of the emmission probabilities as the %start point of the time window is moved back one step at a %time. This is the probability of seeing all the observations %in the analysis window in state j: if(emission_probs == 0 || isnan(emission_probs)) emission_probs =realmin; end %Find the total probability of transitioning from the last %state to this one, with the observations and being in the same %state for the analysis window. This is the duration-dependant %variation of equation 33a from Rabiner: % fprintf('log((duration_probs(j,d)./duration_sum(j))):%d\n',log((duration_probs(j,d)./duration_sum(j)))); delta_temp = max_delta + (emission_probs)+ log((duration_probs(j,d)./duration_sum(j))); %Unlike equation 33a from Rabiner, the maximum delta could come %from multiple d values, or from multiple size of the analysis %window. Therefore, only keep the maximum delta value over the %entire analysis window: %If this probability is greater than the last greatest, %update the delta matrix and the time duration variable: if(delta_temp>delta(t,j)) delta(t,j) = delta_temp; psi(t,j) = max_index; psi_duration(t,j) = d; end end end end end %% Termination % For the extended case, need to find max prob after end of actual % sequence: % Find just the delta after the end of the actual signal temp_delta = delta(T+1:end,:); %Find the maximum value in this section, and which state it is in: [~, idx] = max(temp_delta(:)); [pos, ~] = ind2sub(size(temp_delta), idx); % Change this position to the real position in delta matrix: pos = pos+T; %1) Find the last most probable state %2) From the psi matrix, find the most likely preceding state %3) Find the duration of the last state from the psi_duration matrix %4) From the onset to the offset of this state, set to the most likely state %5) Repeat steps 2 - 5 until reached the beginning of the signal %The initial steps 1-4 are equation 34b in Rabiner. 1) finds P*, the most %likely last state in the sequence, 2) finds the state that precedes the %last most likely state, 3) finds the onset in time of the last state %(included due to the duration-dependancy) and 4) sets the most likely last %state to the q_t variable. %1) [~, state] = max(delta(pos,:),[],2); %2) offset = pos; preceding_state = psi(offset,state); %3) % state_duration = psi_duration(offset, state); onset = offset - psi_duration(offset,state)+1; %4) qt(onset:offset) = state; %The state is then updated to the preceding state, found above, which must %end when the last most likely state started in the observation sequence: state = preceding_state; count = 0; %While the onset of the state is larger than the maximum duration %specified: while(onset > 2) %2) offset = onset-1; % offset_array(offset,1) = inf; preceding_state = psi(offset,state); % offset_array(offset,2) = preceding_state; %3) % state_duration = psi_duration(offset, state); onset = offset - psi_duration(offset,state)+1; %4) % offset_array(onset:offset,3) = state; if(onset<2) onset = 1; end qt(onset:offset) = state; state = preceding_state; count = count +1; if(count> 1000) break; end end qt = qt(1:T);