various commercial applications provide cardiologists with tools that deliver diagnostic parameters , including heart rate variability ( hrv ) [ 13 ] , heart rate turbulence ( hrt ) , t - wave alternans ( twa ) [ 57 ] and many more . the above - mentioned commercial tools , however , often lack the ability to automatically compare the results with other parameters or to verify the applied methods , thus presenting doctors with a final result without being able to examine or adjust the employed algorithms , and thereby precluding further studies . in order to perform diagnoses based on various parameters , some of these tools require advanced programming skills such as c - language programming to take advantage of the offered open source procedures . as a staring point for description of the program reported in this work we provide a short overview of the available toolkits for hrv , hrt and twa   analyses . the heart rate variability ( hrv ) all listed packages offer both time and spectral analyses of hrv , but only some of them implement methods derived from non - linear dynamic techniques . only the finnish kubios hrv analysis software is equipped with a graphical user interface ( gui ) . the toolkit comes with an extensive user guide and provides user support via e - mail . heart rate turbulence ( hrt ) denotes transient fluctuations in the heart rate that directly precede and follow a ventricular premature beat ( vpb ) . the term hrt was coined in 1999 by schmidt et al , who defined the method in a seminal paper published in lancet . currently , many clinical studies have confirmed the significance of hrt as an important predictor of sudden cardiac death risk in patients with myocardial infarction . the computations involved in measurement of hrt parameters are quite simple and rely on measuring geometric relations of the hr signal in the vicinity of the vpb beat . recently , guidelines for hrt measurement standards and physiologic interpretation in clinical studies were published under the supervision of the international society for holter and noninvasive electrophysiology ( ishne ) . in the case of t - wave alternans ( twa ) , to the authors best knowledge , there is only one twa analyser toolkit made available as open source software . the toolkit that can be accessed from the physionet web - page   is a suite of matlab procedures for t - wave analysis , which implement the spectral method ( sm ) and the modified moving average method ( mma ) . a possible explanation for such limited availability of twa analysis software tools is the fact that this cardiologic technique requires special ecg measurement protocols and apparatus to detect microvolt t - wave alternans . however , a considerable number of twa analysis methods have been proposed . for an excellent and comprehensive review of 12 methods see martinez and olmos . in 2010 the u.s . food and drug administration gave market clearance for the microvolt twa spectral method as a viable measure for risk of sudden cardiac death ( scd ) . the indicated indices that proved important in risk stratification of patients ( e.g. , post- infarction patients ) , unfortunately , are not computed in a single software environment or tool . this makes an overall evaluation of a patient s cardiac state a time consuming task . it is also difficult to find the relationship between the parameters since the examination registrations are not simultaneous and the conditions may vary significantly . in this paper a computer program featuring a unified graphical user interface ( gui ) for analysis of hrv , hrt and twa is presented . all of the parameters are computed from the ecg recordings or , alternatively in the case of hrv , from previously computed rr series . at its current development stage , the software is designed predominantly for the purpose of scientific research rather than clinical diagnosis . different file formats of ecg recordings can be loaded into the program . the user is presented with a graphical interface that enables plotting of ecgs and visual verification of intermediate steps of the analyses . the test version of the program can be downloaded without charge from http://www.eletel.p.lodz.pl / kudrynski / matlab/. in order to run the program it is necessary to install matlab compiler runtime ( mcr ) first . the recorded ecg signals or pre - computed rr series can be loaded from binary , text or system - specific format files ( e.g. digital ecg medea , gliwice , poland , and finometer , finapres medical systems , amsterdam , the netherlands ) . suitable filtering procedures are available for denoising the recorded ecgs and beat - by - beat synchronous averaging . baseline wandering removal is implemented by detecting signal samples corresponding to pq segments followed by cubic spline interpolation [ 2022 ] . a haar filter of 20 order is a bandpass linear - phase filter whose passband is centered at f=18.5 hz ( for sampling frequency fs=500 hz ) , while its first and second zeros correspond to 0 hz and 50 hz , respectively . therefore , such a filter retains signal frequency components of qrs complexes and suppresses other signal components such as base - line wander and 50 hz interference coming from the electric powerline network . the reference pq - interval points are found 66msec prior to the q - wave s maximum downslope . the software enables plotting ecg records , visualizing qrs detections , computing and plotting hrv series . the current version of the software enables evaluation of parameters that are claimed to be important for scd risk stratification , mortality after myocardial infarction and for assessment of autonomic system dysfunction in patients with diabetic neuropathy or hypertension . however , further studies are required to evaluate the true validity of these parameters in clinical practice . t - wave alternans is a phenomenon of beat - to - beat alternation of t - wave amplitude and duration occurring at an increased heart rate ( up to 110 bpm ) . numerous studies have established that twa analysis is a useful tool for prediction of arrhythmias and scd risk stratification . the dominant mechanism is dispersion of depolarisation , which can be induced either functionally during rapid heart pacing or when metabolic impairment of ion pumps causes regional changes in ionic gradients as a result of acute myocardial ischemia . therefore , being invisible to the unaided eye , microvolt twa must be analyzed using computer algorithms . in the designed program , twa analysis is performed both in frequency   and in time domain by applying poincare maps . the t - wave period for analysis the selection is marked graphically on the plot of ecg cycles superposed and aligned to previously determined qrs complexes . the amplitude of the t - wave is measured for 7 equally spaced time instances  the first corresponding to the start , and the last to the end of the t - wave selected by the user as shown in figure 2 . in the next step , for each of the seven t - wave phases , a series xn of successive samples is constructed . in the frequency domain analysis method , the level of the amplitude spectrum at frequency equal to 0.5 cycles per beat indicates the level of beat - to - beat alternation of the t - wave amplitude at a given phase . the alternans voltage and ratio indices are calculated according to the formulae 1 and 2 , respectively : valt alternans voltage , s0.5  level of spectrum at frequency 0.5 cycles per beat , snb and nb mean level and standard deviation of spectrum at the noise band defined as 0.440.49 cycles per beat . for the time domain analysis the scatter graph of the poincare mapping the scatter graph is a plot of xn+1 versus xn , where x is the first order differential of series xn . in case of ecg registration of a patient unaffected by twa , the points on the plot form an unstructured cluster otherwise , the points are situated in the second and fourth quarter of the scatter graph ( figure 3b ) . the dtwa index , defined in equation 3 , can be calculated as the distance between the two clusters formed by even and odd scatter graph point pairs . values of both spectral valt and temporal dtwa are computed for all six precordial leads . by using the poincar mapping as well as spectral analysis methods , we have constructed spatio - temporal imaging of twa content across v1v6 then , in order to generate a smooth image of such a spatio - temporal distribution of twa , a biharmonic spline interpolation technique is employed for spanning an interpolation surface supported by 42 knots at space - time locations of the precomputed dtwa indices . example spatio - temporal maps for twa+ patients are shown in figure 3c , d . hot colors indicate regions of high twa levels that can be used for both spatial ( in terms of precordial lead location ) and temporal ( t - wave section number ) identification of heart repolarization abnormalities . heart rate turbulence is a response of the sinus node to the ectopic beat , characterized by rapid acceleration prior to ventricular premature beat ( vpb ) , followed by deceleration of the heart rate . the most probable mechanism underlying hrt is the baro - reflex autoregulation process , which responds to the deviation in blood pressure caused by ectopic beat occurrence . missing or weakened response of such a nature reflects a malfunctioning of the autonomic control system and is considered a powerful risk indicator in patients after myocardial infarction . in the designed program , 4 ) as the mean percentage difference between the average interbeat interval ( ibi ) after and before a ventricular premature beat . where e[. ] denotes the expectation operator and the turbulence slope ( ts ) is the maximum slope of the linear regression line for each sequence of 5 consecutive normal intervals taken from the 15 ibis following the vpb , providing they are free from arrhythmia and artifacts . because the hrt pattern can be masked by heart rate variability , these parameters are usually calculated as an averaging responses to > 5 ventricular premature contractions . heart rate variability ( hrv ) is a measure of beat - to - beat variations in the length of the rr interval , and is a marker of autonomic nervous system modulation of the sinus node activity . limited studies have linked short - term hrv to sudden death , and its use for risk stratification is currently not recommended . the designed program enables computation of statistical parameters in time domain ( mnn , sdnn , rmssd , pnn50 ) and a number of parameters in the frequency domain ( lf and hf power , lf / hr ratio ) . the mnn and sdnn are mean and standard deviation of normal - to - normal rr intervals , rmssd is the root mean square of the differences between consecutive rr intervals , and pnn50 is the percentage of a number of consecutive rr intervals differing by more than 50 ms . spectral analysis based on the fourier transform ( ft ) requires preprocessing of the rr interval sequence . artefacts , including ectopic beats , and missed and false detections , are located in the series by applying thresholding to the second derivative of the hrv series and then removed . in order to obtain uniformly sampled series , the rr samples are interpolated and the result is re - sampled with an appropriate frequency . a frequency of 2hz is assumed as the maximum possible frequency of occurrence of the heart beat event . interpolation is performed using the cubic spline method , which is based on piecewise polynomial functions of the third order . the frequency characteristic is determined using the parametric autoregressive moving average ( arma ) parametric model and nonparametric ( periodogram ) model , ( i.e. the squared value of the fourier transform of the analysed hrv ) . in the implemented welch periodogram method , the signal is divided into a number of overlapping intervals and multiplied by a window function prior to ft analysis . the periodograms are computed for each interval and averaged . by applying such a procedure , the obtained spectrum is smoother than in the case of classic periodogram . by adjusting the degree of overlap , the number of intervals , and type of the window , it is possible to improve the resolution of the spectral characteristic at the cost of its smoothness . in non - parametric methods the frequency characteristics are computed directly from the examined signal . in parametric methods , on the other hand since the heart s natural pacemaker follows complex periodic dynamics , many simplifications need to be introduced into the model . nevertheless , a good approximation of the underlying mechanisms controlling heart rate variability can be achieved with the use of such models . in the moving average modeling ( ma ) approach , the hrv series is modelled by a weighted linear combination of s(k ) , k=1, the modeling task is to determine the vector of weight coefficients bk that would generate the best fit to the analyzed hrv series y(n ) according to the formula : one of the model design parameters is the order of the model given by q. too low an order of the model would omit important hrv details , whereas too high an order can yield a noisy series that would be difficult for meaningful interpretation . in the autoregressive modeling ( ar ) approach , the hrv series is modeled by a recursive process in which the current model sample is generated from a predetermined number of its previous samples ( hence the model name  autoregressive  ) multiplied by appropriate coefficients ak , k=1 ,  , p. the corresponding equation is the following : again , the model order p in an important parameter that has to be decided in order to obtain a suitable ar model of the hrv series . the ma model is driven by gaussian noise and omits inherent dynamics of the modelled object . on the other hand , the ar model does not take into account any stochastic mechanisms that influence the hrv . hence , a modeling approach that combines the two models is a frequently adopted technique in hrv parametric modeling . the resulting model is termed the autoregressive moving average model ( arma ) and the resulting model equation assumes the form : here , the goal of the modelling process is to find vectors of coefficients ak and bk so that the output y(n ) of the system would be as close as possible to the examined hrv series . this is a challenging task because in practice 1020 parameters need to be simultaneously determined . there have been , however , special algorithms developed to solve this multicriterion optimisation problem . another , indicated problem is the selection of the orders p , q of the model . assuming the same orders for all the patients or even for the same patient at various conditions should be considered too far a simplification . the algorithms implemented in the presented software provide means for automatic computation of the model order by optimising model complexity versus the modelling error . the implemented order selection algorithm is based on the so - called akaike s information criterion , which indicates the goodness of fit of the model to the experimentally measured hrv series . the interested reader may probe further into ar , ma and arma modelling by studying classic texts on the subject [ 26 , 28 ] . t - wave alternans is a phenomenon of beat - to - beat alternation of t - wave amplitude and duration occurring at an increased heart rate ( up to 110 bpm ) . numerous studies have established that twa analysis is a useful tool for prediction of arrhythmias and scd risk stratification . the dominant mechanism is dispersion of depolarisation , which can be induced either functionally during rapid heart pacing or when metabolic impairment of ion pumps causes regional changes in ionic gradients as a result of acute myocardial ischemia . therefore , being invisible to the unaided eye , microvolt twa must be analyzed using computer algorithms . in the designed program , twa analysis is performed both in frequency   and in time domain by applying poincare maps . the selection is marked graphically on the plot of ecg cycles superposed and aligned to previously determined qrs complexes . the amplitude of the t - wave is measured for 7 equally spaced time instances  the first corresponding to the start , and the last to the end of the t - wave selected by the user as shown in figure 2 . in the next step , for each of the seven t - wave phases , a series xn of successive samples is constructed . in the frequency domain analysis method , the level of the amplitude spectrum at frequency equal to 0.5 cycles per beat indicates the level of beat - to - beat alternation of the t - wave amplitude at a given phase . the alternans voltage and ratio indices are calculated according to the formulae 1 and 2 , respectively : valt alternans voltage , s0.5  level of spectrum at frequency 0.5 cycles per beat , snb and nb mean level and standard deviation of spectrum at the noise band defined as 0.440.49 cycles per beat . for the time domain analysis the scatter graph of the poincare mapping the scatter graph is a plot of xn+1 versus xn , where x is the first order differential of series xn . in case of ecg registration of a patient unaffected by twa , the points on the plot form an unstructured cluster . otherwise , the points are situated in the second and fourth quarter of the scatter graph ( figure 3b ) . the dtwa index , defined in equation 3 , can be calculated as the distance between the two clusters formed by even and odd scatter graph point pairs . values of both spectral valt and temporal dtwa are computed for all six precordial leads . by using the poincar mapping as well as spectral analysis methods , we have constructed spatio - temporal imaging of twa content across v1v6 then , in order to generate a smooth image of such a spatio - temporal distribution of twa , a biharmonic spline interpolation technique is employed for spanning an interpolation surface supported by 42 knots at space - time locations of the precomputed dtwa indices . example spatio - temporal maps for twa+ patients are shown in figure 3c , d . hot colors indicate regions of high twa levels that can be used for both spatial ( in terms of precordial lead location ) and temporal ( t - wave section number ) identification of heart repolarization abnormalities . heart rate turbulence is a response of the sinus node to the ectopic beat , characterized by rapid acceleration prior to ventricular premature beat ( vpb ) , followed by deceleration of the heart rate . the most probable mechanism underlying hrt is the baro - reflex autoregulation process , which responds to the deviation in blood pressure caused by ectopic beat occurrence . missing or weakened response of such a nature reflects a malfunctioning of the autonomic control system and is considered a powerful risk indicator in patients after myocardial infarction . in the designed program , 4 ) as the mean percentage difference between the average interbeat interval ( ibi ) after and before a ventricular premature beat . where e[. ] denotes the expectation operator and the turbulence slope ( ts ) is the maximum slope of the linear regression line for each sequence of 5 consecutive normal intervals taken from the 15 ibis following the vpb , providing they are free from arrhythmia and artifacts . because the hrt pattern can be masked by heart rate variability , these parameters are usually calculated as an averaging responses to > 5 ventricular premature contractions . heart rate variability ( hrv ) is a measure of beat - to - beat variations in the length of the rr interval , and is a marker of autonomic nervous system modulation of the sinus node activity . limited studies have linked short - term hrv to sudden death , and its use for risk stratification is currently not recommended . the designed program enables computation of statistical parameters in time domain ( mnn , sdnn , rmssd , pnn50 ) and a number of parameters in the frequency domain ( lf and hf power , lf / hr ratio ) . the mnn and sdnn are mean and standard deviation of normal - to - normal rr intervals , rmssd is the root mean square of the differences between consecutive rr intervals , and pnn50 is the percentage of a number of consecutive rr intervals differing by more than 50 ms . spectral analysis based on the fourier transform ( ft ) requires preprocessing of the rr interval sequence . artefacts , including ectopic beats , and missed and false detections , are located in the series by applying thresholding to the second derivative of the hrv series and then removed . in order to obtain uniformly sampled series , the rr samples are interpolated and the result is re - sampled with an appropriate frequency . a frequency of 2hz is assumed as the maximum possible frequency of occurrence of the heart beat event . interpolation is performed using the cubic spline method , which is based on piecewise polynomial functions of the third order . the frequency characteristic is determined using the parametric autoregressive moving average ( arma ) parametric model and nonparametric ( periodogram ) model , ( i.e. the squared value of the fourier transform of the analysed hrv ) . in the implemented welch periodogram method , the signal is divided into a number of overlapping intervals and multiplied by a window function prior to ft analysis . the periodograms are computed for each interval and averaged . by applying such a procedure , the obtained spectrum is smoother than in the case of classic periodogram . by adjusting the degree of overlap , the number of intervals , and type of the window , it is possible to improve the resolution of the spectral characteristic at the cost of its smoothness . in non - parametric methods the frequency characteristics are computed directly from the examined signal . in parametric methods , on the other hand since the heart s natural pacemaker follows complex periodic dynamics , many simplifications need to be introduced into the model . nevertheless , a good approximation of the underlying mechanisms controlling heart rate variability can be achieved with the use of such models . in the moving average modeling ( ma ) approach , the hrv series is modelled by a weighted linear combination of s(k ) , k=1, the modeling task is to determine the vector of weight coefficients bk that would generate the best fit to the analyzed hrv series y(n ) according to the formula : one of the model design parameters is the order of the model given by q. too low an order of the model would omit important hrv details , whereas too high an order can yield a noisy series that would be difficult for meaningful interpretation . in the autoregressive modeling ( ar ) approach , the hrv series is modeled by a recursive process in which the current model sample is generated from a predetermined number of its previous samples ( hence the model name  autoregressive  ) multiplied by appropriate coefficients ak , k=1 ,  , p. the corresponding equation is the following : again , the model order p in an important parameter that has to be decided in order to obtain a suitable ar model of the hrv series . the ma model is driven by gaussian noise and omits inherent dynamics of the modelled object . on the other hand , the ar model does not take into account any stochastic mechanisms that influence the hrv . hence , a modeling approach that combines the two models is a frequently adopted technique in hrv parametric modeling . the resulting model is termed the autoregressive moving average model ( arma ) and the resulting model equation assumes the form : here , the goal of the modelling process is to find vectors of coefficients ak and bk so that the output y(n ) of the system would be as close as possible to the examined hrv series . this is a challenging task because in practice 1020 parameters need to be simultaneously determined . there have been , however , special algorithms developed to solve this multicriterion optimisation problem . another , indicated problem is the selection of the orders p , q of the model . assuming the same orders for all the patients or even for the same patient at various conditions should be considered too far a simplification . the algorithms implemented in the presented software provide means for automatic computation of the model order by optimising model complexity versus the modelling error . the implemented order selection algorithm is based on the so - called akaike s information criterion , which indicates the goodness of fit of the model to the experimentally measured hrv series . the interested reader may probe further into ar , ma and arma modelling by studying classic texts on the subject [ 26 , 28 ] . a user - friendly software package for quantitative evaluation of hrv , hrt and twa in ecgs was developed . the software tool offers the following main features for ecg processing and analysis : ecg preprocessing tools and file handling tools : loading ecg recordings from text or binary files . plotting the loaded ecg recordings of every lead as well as the result of qrs detections accompanied by the preview of rr series . the preview plots of the currently loaded ecgs are built - in within the program s gui ( figure 1 ) . detection of qrs complexes of the san ( sinoatrial node ) rhythm as well as the rhythm stimulated by artificial pacemakers ( e.g. atrial sensed , ventricular paced rhythm ) ( figure 5 ) . denoising ecgs with synchronized averaging method as shown in figure 8 ( averaging a number of ecg cycles synchronized by the qrs complex ) . such an approach enables reduction of random noise signal components and filtering out ecg characteristic features repeating at every averaged cycle ( e.g. hidden p wave ) . exporting the loaded signals ( either source ecg recordings or preprocessed ecgs ) to .txt and .mat files . the recommended number of 15 artefact - free cycles following vpb can be selected by the user . the values of to and ts parameters are accompanied by the plots presenting the algorithm s detection of vpb beats fulfilling the conditions described in the section about hrt analysis . the user is also provided with the ecg plot where the mentioned detections are marked . after graphical selection of t - wave start and end points ( figure 2 ) the user is provided with both spectral and poincare analysis plots for each of the precordial leads at every time section . the sensitivity of artefact correction can be set by the user , which may be found useful in case of highly contaminated recordings or in situations were certain natural rhythm changes are wrongly interpreted by the algorithm as artefacts . the parameters , including type of analysis , model type , and model order are also user - defined . clinical recordings that have already been examined by this software package were provided by the department of electrocardiology of the medical university in lodz . programs used in this software package were used to assess twa in postinfarction patients with low ejection fraction ( ef% ) and preserved left ventricular function , hrv in postinfarction patients with impaired left ventricular function , and hrv in patients with chronic heart failure and implanted cardiac resynchronization devices . each of the supported analyses can be performed using different methods with various parameters , which makes it possible to cross - validate the results obtained . the derived parameters ( given in time , frequency and spatio - temporal domain ) are visualized for easier interpretation . preview of the intermediate steps of the analysis ( qrs detection , hrv artefact correction , etc . ) is a good and easily available source of information on possible errors such as beat misdetections . thanks to this previewing functionality , the parts of registrations prone to errors ( eg , fragments with very frequent occurrence of ectopic beats ) can be excluded from further analyses . spectra of the series of t - wave amplitudes for recordings from the precordial leads are plotted in figure 3a . results of the temporal analysis obtained by applying the poincare maps are shown in figure 3b  d illustrate the interpolation maps of twa voltage obtained from spectral and spatio - temporal analysis , respectively . similar results can be obtained using the two methods , while the computational cost of the latter is much lower . to achieve interpretable results of the spectral analysis , a relatively long analysis period is required . in spectral analysis of hrv the choice of the method modelling the bands of the smooth spectra is inefficient with the ar method . on the other hand , the moving average ( ma ) arma modelling being the combination of both autoregressive and moving average models is capable of estimating spectra partially rich with peaks and partially smooth . the cardiologist analyzing the results of the hrv modelling should be aware of the properties of these algorithms . the shape of the hrv spectra is only partially dependent on the nature of autonomous nervous system behaviour ; it also depends on the applied model type and model order ( figure 9 ) . in the designed program , an automatic model order selection algorithm is implemented and the results of different methods are available for comparison and cross - validation . although isolated studies have demonstrated the utility of many different variables for sudden cardiac death ( scd ) risk stratification ( imaging - based , ecg - based , exercise test - based , holter - based and others ) , the ideal tool to identify patients who might experience life- threatening arrhythmias has not been found . among the risk factors assessed non - invasively , low left ventricular ejection fraction ( lvef% ) is the most powerful predictor of serious events in ischemic and non - ischemic cardiomyopathy . in fact , it is the sole variable determining the recommendation for implantable cardioverter - defibrillator ( icd ) . however , the majority of scds occur in patients with more preserved or normal ef , and the evaluation of risk in this large group of patients is justified by means of the secondary factors computed within the proposed software tool . assessment of risk factors derived from the most frequently performed examinations , such as supine and exercise ecg , can help to resolve the problem of risk stratification for scd in patients without left ventricular impairment . the developed software package is a flexible and comprehensive tool for twa , hrt and hrv quantitative analysis from ecgs . it offers the possibility to choose from various modelling methods and to modify parameters of the applied algorithms . the interface is user - friendly and enables interactive selection of ecg signal fragments for analysis . the authors believe that the proposed software tool for ecg analyses can offer a basis for a better common understanding of engineers and clinicians working together in the difficult field of cardiology .
summarybackgroundthis paper presents a software package for quantitative evaluation of heart rate variability ( hrv ) , heart rate turbulence ( hrt ) , and t - wave alternans ( twa ) from ecg recordings . 
 the software has been developed for the purpose of scientific research rather than clinical diagnosis.material/methodsthe software is written in matlab mathematical language . 
 procedures for evaluation of hrv , hrt and twa were implemented . 
 hrv analysis was carried out by applying statistical and spectral parametric and nonparametric methods . 
 hrt parameters were derived using the schmidt algorithm . 
 twa analysis was performed both in spectral and in time domain by applying poincare mapping . 
 a flexibility of choosing from a number of classical modelling approaches and their modifications was foreseen and implemented . 
 the software underwent preliminary verification tests both on ecgs from the physionet online ecg signal repository and recordings taken at the department of electrocardiology of the medical university hospital in lodz.resultsthe result of the research is a program enabling simultaneous analysis of a number of parameters computed from ecg recordings with the use of the indicated analysis methods 
 . the program offers options to preview the intermediate results and to alter the preprocessing steps.conclusionsby offering the possibility to cross - validate the results of analyses obtained by several methods and to preview the intermediate analysis steps , the program can serve as a helpful aid for clinicians in comprehensive research studies . 
 the software tool can also be utilized in training programs for students and medical personnel .