A method for noninvasive beat‐by‐beat visualization of His bundle signals

Abstract Background Invasive recording of His bundle signals (HBS) in electrophysiological study (EPS) is important in determining HV interval, the time taken to activate the ventricles from the His bundle. Noninvasive surface measurements of HBS are attempted by averaging typically 100–200 cardiac cycles of ECG time series in body surface potential mapping (BSPM) and in magnetocardiography (MCG) which records weak cardiac magnetic fields by highly sensitive detectors. However, noninvasive beat‐by‐beat extraction of HBS is challenged by ramp‐like atrial signals and noise in PR segment of the cardiac cycle. Methods By making use of a signal‐averaged trace showing prominent HBS as a guide trace, we developed a method combining interval‐dependent wavelet thresholding (IDWT) and signal space projection (SSP) technique to eliminate artifacts from single beats. The method was applied on MCG recorded on 21 subjects with known HV intervals based on EPS and noninvasive signal‐averaging, including five subjects with BSPM recorded subsequently. The method was also applied on stress‐MCG of a subject featuring autonomic dynamics. Results HBS could be extracted from 19 out of 21 subjects by signal‐averaging whose timing differed from EPS between −8 and 11 ms as tested by 2 observers. HBS in single beats were seen as aligned patterns in inter‐beat contours and were appreciable in stress‐MCG and conspicuous than BSPM. The performance of the method was evaluated on simulated and measured MCG to be adequate if the signal‐to‐noise ratio was at least 20 dB. Conclusions These results suggest the use of this method for noninvasive assessments on HBS.


| INTRODUC TI ON
The measurement of the time difference of the onsets of the spikelike deflection of the His bundle (H) and the QRS complex (V) in the P wave-to-R wave interval (PR interval) of the cardiac cycle (Muresan et al., 2019) is a quantitative index of atrio-ventricular conduction.
Owing to the fact, the His bundle signals (HBS) are significantly weaker than the rest of the features of the cardiac cycle; its recording necessitates an invasive intra-cardiac electrophysiological study (EPS) performed as a catheter-based interventional procedure.
HBS are bi or triphasic deflections ~100 Hz, and the HV interval is determined as the time difference between the onsets of the first wave after the intracardiac recordings of atrial signal and that of the surface recording of the Q wave seen in ECG. The normal range of HV interval is between 35 and 50 ms (Miller et al., 2018). A prolongation in the HV interval exceeding 75 ms hints at an eventual risk of interrupted electrical activation of the ventricles, needing implantation of a permanent artificial pacemaker to avoid complete heart block (Miller et al., 2018). Among the intracardiac signals and their time intervals that could be recorded in an EPS, HBS have a scope to be recorded by noninvasive methods by averaging several beats that are time marked to a chosen fiducial point on the cardiac cycle (Berbari et al., 1973;Flowers et al., 1974;Hombach et al., 1982;Vincent et al., 1978). HV interval is reasonably constant and normally maintains a temporal relationship with the QRS-onset and is relatively unaffected by cardiac dynamics under normal conditions unlike other intervals of the cardiac cycle (Miller et al., 2018;Muresan et al., 2019). Hence, signal averaging procedure has been widely employed in investigations on noninvasive measurements of HBS in voltage-based techniques such as the ECG/body surface potential mapping (BSPM) measurements (Berbari et al., 1973;Flowers et al., 1974;Hombach et al., 1982;Vincent et al., 1978) and in magnetocardiography (MCG), which records subtle cardiac magnetic fields (~100 femto Tesla-100 pico Tesla) (Tavarozzi et al., 2002).
Among different magnetic field detectors and their varying configurations, MCG measured using superconducting quantum interference devices (SQUIDs) operating at 4.2 K temperature inside a magnetically shielded room (MSR) is regarded as a gold standard technique to record clinical MCG (Fenici et al., 1983;Sengottuvel et al., 2017;Senthilnathan et al., 2012;Tavarozzi et al., 2002;Yamada et al., 2003). MCG has been reported to offer complimentary diagnostic information on cardiac electrophysiology which is superior to that offered by voltage-based measurements in a variety of cardiac dysfunctions (Tavarozzi et al., 2002). The list includes noninvasive measurement of highly resolved His-signal features and HV intervals comparable to that of an EPS (Fenici et al., 1983;Sengottuvel et al., 2017;Senthilnathan et al., 2012;Yamada et al., 2003).
Biological artifacts such as very low-frequency ramp-like signals which are believed to be attributed to atrial repolarization in the PR interval undermine HBS in noninvasive measurements (Fenici et al., 1983;Sengottuvel et al., 2017) and further in retrieving them in a beat-by-beat manner. These low-frequency signals are circumvented in invasive EPS measurements using filters with a typical bandwidth setting of 30-300 Hz (Muresan et al., 2019).
Furthermore, HBS are reliably validated by moving the invasive measurement probe closer to the His bundle guided by X-ray fluoroscopy in EPS. Because such maneuvers are not possible in noninvasive surface measurements, they are posed with an additional difficulty for rightly discriminating HBS from distortions and artifacts which are sometimes introduced by the signal conditioning methods itself (Feyter et al., 1981;Sengottuvel et al., 2017;ten Voorde et al., 1988). In order to unambiguously identify His bundle signals in a more reliable way, few researchers used multichannel ECG/MCG to seek additional information on the signal registrations of the His bundle by inspecting the spatial distribution patterns of signals obtained from iso-field contour maps of multichannel data (Farrell et al., 1980;Fenici et al., 1983;ten Voorde et al., 1988). This suggestion had helped to better define HBS in noninvasive surface measurements and in determining the HV intervals comparable to invasive EPS traces . But that study could extract HBS only on signal-averaged beats of the whole cardiac time series for each measurement channel (one representative beat per measurement channel), because of the use of signal averaging which helped to reduce inter-beat noise. Recently, a technique on reducing beat-by-beat noise by treating the cardiac time series, as beat epochs was reported and was demonstrated to be helpful in the identification of MCG cardiac features based on inter-beat similarities and time registries of cardiac events even under extremely noisy measurement conditions . It was intuitively persuading to employ a similar procedure for identifying HBS in every beat, by utilizing the relative constancy of HV interval on cases for whom HBS had already been identified based on multichannel ECG/ MCG. The present work had drawn motivation from similar problems in neuroscience of recording event-related potentials (ERP) in an electroencephalogram (EEG), by repeatedly presenting a specific stimulus to subjects and extracting the brain response pertaining only to the stimulus by trigger-locked averaging (Luck, 2014). Quiroga and Garcia (2003) demonstrated that by using the triggerlocked averaged ERP waveform as a guide trace, single-trial brain responses could be retrieved by interval-dependent wavelet thresholding (IDWT) by masking noise that occur other than the signal regime. Furthermore, the dynamics in neural responses were also visualized with the help of inter-trial contour maps (Rey et al., 2015).
The present work hence attempted to devise a method which would facilitate tracking HBS across beats by the conventional technique of averaging multiple beats to obtain HBS and use it as a template to define the signal regime for suppressing ramps and inter-beat noise in the PR interval using IDWT. In addition, signal space projection (SSP) technique (Kunpeng et al., 2012) which is suitable to cancel common mode signal (i.e., ramp-like variations) was employed in this work. Techniques on extracting HBS in every cardiac beat have not been reported so far in MCG and such studies in BSPM or ECG are also only sparingly discussed in the literature using proprietary hardware and software schemes (Ishijima, 1977;Liu et al., 2015;Wang et al., 2017). Moreover, the reported objective criteria in those ECG measurements involve inferring HBS from multiple micro-wavelets in the PR interval necessitating subjective bias/practice in their identification (Liu et al., 2015;Wang et al., 2017). Hence, a signal treatment scheme is proposed on unambiguously validated HBS (by comparison with EPS), which were obtained by passively averaging individual MCG beats, and thereby, the proposed method has a scope to inspect their inter-beat dynamics which were lost in the process of averaging. Thus, the method is expected to augment the visualization of HBS by tracking its components in a beat-by-beat manner in MCG and BSPM time series. Furthermore, the present work also demonstrates a possible application of the method in tracking HBS in an exercise MCG measurement where the beat-bybeat cardiac dynamics widely varied across the cardiac time series.

| Instrumentation and study protocol
The measurement data sets used in this investigation were recorded as a part of two independent studies (Kesavaraja et al., 2021;Sengottuvel et al., 2017) and were retrospectively analyzed for the present work. Institutional ethics committee approval and informed consent from subjects had been obtained for the original investigation. MCG data of 21 subjects (48 ± 11 years) were analyzed in this work. Out of these subjects (12 male, 9 female), two had ventricular tachycardia, and others had supraventricular tachycardia diagnosed based on ECG and EPS. The measurement of the HV interval formed a part of the usual detailed evaluation in EPS. The EPS was measured following the conventional clinical procedure in a hospital setting. MCG measurements were done in a different institution by an independent research group on different days. All these subjects were hemodynamically stable and were accompanied by a cardiologist when brought for MCG measurements. For five of these subjects, BSPM was recorded subsequently after MCG measurements.
Stress-MCG was measured from a healthy control (M33) involved in physical exertion using a non-magnetic bicycle ergometer attached to the MCG system. MCG was measured using a thirty-seven-channel SQUID system operating inside a low noise environment facilitated by a magnetically shielded room (MSR) (Parasakthi et al., 2012). The room contained two layers of mu-metal and two layers of Aluminium capable of attenuating external electromagnetic interferences with a shielding factor better than 100 dB at 100 Hz, 70 dB at 1 Hz and 50 dB at 0.1 Hz. The MCG measurement system consisted of a hexagonal array of low transition temperature, direct current biased SQUID sensors (Parasakthi et al., 2012). The sensors were housed in a flat-bottomed liquid Helium cryostat made of fiber reinforced plastic. The warm-to-cold distance of the cryostat was 1.5 cm which dictated the minimum most distance of separation between the chest of subjects and the sensor plane. The SQUID electronics was operated with a bandwidth setting of 0-300 Hz. The subjects were positioned in supine posture with the chest facing the bottom tail of the suspended cryostat mounted on a non-magnetic gantry. The chest of the subjects was aligned under the cryostat about anatomical landmarks such as the collar bone, sternal line, fourth intercostal space, and xiphoid process. The average noise floor measured across the 37 channels in the MCG system was 25 ± 9 fT rms /√Hz. For the healthy control subject involved in stress measurement, the MCG measurement was made during a pedaling action in supine posture in the ergometer, creating a moderate stress condition (Kesavaraja et al., 2021). BSPM measurements were made using a sixty-four-channel EEG system (Compumedics Ltd.) which can also function as a generalpurpose bio-potential recorder (Patel et al., 2016). The amplifier was set for a gain of 1000, bandwidth set at 0-1000 Hz (and later bandlimited to 300 Hz). The digitization accuracy of the system was 30 nV per least significant bit, and the peak-peak noise amplitude of the system was <1 μV pp . BSPM measurements were made precisely on the same thirty-seven locations on the chest as measured by the MCG system. Discrete Ag-AgCl electrodes were used to measure the potential difference across the chest locations. The signal processing methods and their sequence were the same for both MCG and BSPM cardiac time series.   such that at least 200 consecutive cardiac cycles (beats) are taken for further analysis. The PR intervals of all 200 cardiac beats were isolated and averaged to get one representative signal-averaged trace. Each of these averaged waveforms showed a P wave followed by ramp-like slowly varying signals, a small bump and again ramps up to the QRS-onset in all the measurement channels but with differing magnitudes and configurations across spatial locations. The objective criteria followed in this work to define a bump (or a deflection) in the signal-averaged PR interval waveforms as the His bundle signal are as follows:

| Presentation of MCG/BSPM time series
a. The chosen deflection in the PR segment was above the noise floor that was observed before the onset of the P wave, where no cardiac activity could be expected.
b. The chosen bump was ensured to be well delineated from the offset of the P wave observing its time registries in all the channels.
c. The chosen bump/deflection was surrounded by ramp-like signals before and after its occurrence. d. The later part of the ramp that followed a chosen bump changed its polarity across the measurement channels, i.e., that which changed from a linearly increasing to decreasing variation and vice-versa.
This physiological sequence of events in the PR interval was wellstudied by researchers (Fenici et al., 1983;ten Voorde et al., 1988) and was utilized to identify HBS in an earlier work . The ramp-like signals were reported to be attributed to atrial repolarization which was known to be dominant in the initial portion of the PR interval and progressively subdued at the later part after His bundle deflection (Fenici et al., 1983). The objective criteria and these descriptions could be comprehended better by plotting an overlay plot of the signal-averaged waveforms as illustrated in Figure  required to test the objective criteria. As could be seen in the magnified view of the PR interval in Figure 2b, a bump-like deflection rides over ramps in the waveforms of all the channels. Two vertical lines (T1 and T2) were manually drawn based on visual inspection such that they were placed before and after a chosen bump in the PR segment to identify the physiological origin of the ramps under discussion. As could be noticed in the figure, the ramps at time instant T2 changes its polarity across measurement channels and thus fulfills the criteria to define the chosen bump as HBS. The identified HBS from signal-averaged traces in MCG and their respective HV intervals which were noninvasively deduced were verified at a later stage with their corresponding EPS report in a blind-folded manner by two independent researchers after MCG measurements and their analysis were completed. Detailed descriptions of this analysis are given in an earlier work . Furthermore, the HBS and the HV intervals of a subset of cases used in this analysis were revalidated in an unbiased manner by two electrophysiologists as independent observers who were unaware of the proposed method. The two observers were requested to identify HBS and to determine HV intervals in the PR interval and the values were compared against the HV intervals determined already by us as a reference set. The morphology of the obtained HBS in MCG varied from a sharp spike ("M" shaped biphasic deflection) to a bump, riding over the ramps. In some cases, the bumps were superposed with highfrequency wiggles which were attributed to myoelectric noise or other random fluctuations in the PR interval (Hombach et al., 1982;Yamada et al., 2003). The signal-averaged MCG trace with verified HBS was then subjected to interval dependent wavelet thresholding and signal space projection technique. The mathematical basics of IDWT and SSP are given elsewhere (Donoho, 1995;Uusitalo & Ilmoniemi, 1997) the following details describe the relevance and the way of implementation of these techniques in the present work to extract HBS in every beat.

| Wavelet thresholding
Wavelet thresholding aimed to reduce the background signals present in the PR interval other than the time segment of HBS (in the HV interval). Wavelet Toolbox of MATLAB software, version R2013a (The MathWorks Inc.) was used for IDWT. The averaged signal was decomposed to a set of approximation and detail coefficients by 'coif 5' wavelet in five levels of decomposition. The choice of this wavelet basis was based on a report which discussed its suitability for analyzing His bundle signal features (Lee, 2004).
The magnitude values of each of the detail coefficients were thresholded by fixed form soft thresholding by keeping the time intervals around the chosen HBS unaffected (few milliseconds of time points on both sides of the HBS bump or deflection) by manually setting the threshold intervals (Donoho, 1995). The number of intervals was fixed based on the discernibility of HBS overlapping with the preceding and succeeding ramps in the PR interval. For the cardiac signals used in the present work, either 4 or 5-time F I G U R E 1 Block diagram of the proposed beat-by-beat analysis method for the visualization of His bundle signals.
intervals were found to be sufficient for IDWT. The chosen time intervals with the set of threshold parameters were stored as a function and applied on all the individual raw epochs for a given subject. Because thresholding technique (Donoho, 1995) could be performed only on the detail coefficients, the attenuation of the low-frequency ramps which fell outside the chosen region of interest (within the HV interval) was taken care of separately by the SSP technique.

| Signal space projection technique
SSP aimed at separating HBS from the PR interval by creating an orthonormal basis function of the atrial components generated from the PR epochs. The wavelet thresholded epochs of the PR intervals were smoothed by applying Savitzky-Golay (SG) filter which acted as a local least square polynomial estimator (Hargittai, 2005).
The order of the SG filter was set as 3 (cubic) with a frame length of 23 to be suitable for over-smoothing operations (Krishnan & Seelamantula, 2013). Since the HBS are high frequency signals, the smoothed epochs were expected to contain only atrial components: the P wave and ramps in the PR interval.
The measured signal in the PR interval is a combination of physiological events of the HBS and the atrial activity (comprising depolarization and repolarization of atria) as represented by Equation (1). It is reasonable and physiologically consistent to assume that both signals are from different sources (Berbari et al., 1973;Flowers et al., 1974;ten Voorde et al., 1988;Vincent et al., 1978).
If b A (t) was characterized by magnetic field/electric potential waveforms of 'm' number of epochs of a chosen measurement channel, then a matrix U n with each column forming an orthonormal basis vector could be constituted by singular value decomposition (Uusitalo & Ilmoniemi, 1997). The vector space spanned by U n was designated as the atrial subspace. The size of U n was less than 'm' comprising only significant basis vectors, which was sufficient to represent b A (t). It was then possible to form a subspace P per that was orthogonal to the atrial subspace and was devoid of atrial components as given by Equation (2). where, I is an identity matrix of the same size as U n .
By projecting the measured signal over this new subspace P per, the His bundle signals could be selectively retrieved from the PR interval epochs as given by Equation (3).
Detailed descriptions of SSP technique and its use for removing common-mode background artifacts are given elsewhere (Sriram et al., 2013).
The extracted HBS were visualized by inter-beat contour maps generated throughout the MCG epochs of 200 consecutive beats.
However, for better visibility of the HBS, the beat contours are illustrated in this work for a time segment of length ~25 beats.

| Simulation study on the proposed method in retrieving HBS
The proposed analysis method was assessed for its ability to extract HBS under varying levels of atrial influence and random noise using a synthetic signal representing a PR interval epoch containing F I G U R E 2 Identification of His bundle signals from overlay plot of signal-averaged MCG waveforms across 37 measurement locations (each colored line representing MCG measurement channel) (a) PR interval of signal-averaged MCG waveforms showing cardiac deflections of interest namely P wave, offset of P wave and QRS-onset. Time instants T1 and T2 are manually marked, respectively, after P wave offset and before QRS-onset on ramp-like signals encompassing a bump (b) Magnified view of (a) showing a bump attributed to HBS (H) and the ramp marked at T2 changes polarity across measurement locations. HV interval computed as the difference between the QRS-onset and T1 is also indicted.
a reference HBS at an interval of 60 ms before the QRS-onset (i.e, HV = 60 ms). Varying magnitudes and configurations of ramps mimicking atrial repolarization activity and fluctuations caused by random noise which corrupt the PR interval were simulated. This simulation involved generating three different types of atrial T wave (PT a ) (atrial repolarization wave) ramps that commonly occur in the PR interval (Farrell et al., 1980;Jayaraman et al., 2016). Atrial repolarization waveform with its magnitude one-third of that of the P wave, but of opposite polarity (designated as Type 1), PT a with magnitude one-half of that of the P wave (designated as Type 2), and PT a with a magnitude exceeding one-half of the magnitude of P wave in the positive direction (designated as Type 3) were used. Additive white Gaussian noise was added to each type of the simulated PR interval waveform with signal-to-noise ratio (SNR) ranging between 5-28 dB (with P wave peak regarded as the signal magnitude). Pearson's correlation coefficient was calculated between a reference HBS (containing only the HBS without random noise and atrial components) and the retrieved HBS using the proposed method in each type of the simulated waveform.

| Quantitative parameters computed from the retrieved HBS in the measured MCG
Two quantitative parameters namely signal-to-noise ratio (SNR) and signal-to-error ratio (SER) (Rajesh et al., 2017), the former in time and the latter in frequency domain were used to evaluate the efficiency of the proposed method in reducing the atrial activity and random noise. SNR values were calculated at P wave peak and that of extracted HBS against the noise present before the onset of the P wave in the PR interval of a signal-averaged beat.
These values were compared against those calculated for one of the beat epochs in which the HBS were extracted using the proposed IDWT-SSP method. SER was computed using the following Equation (4).
SER was expected to be higher at the frequency range of signals signifying the preservation of the content of the signal and was expected to be least when calculated at unwanted noise frequencies indicating their envisaged reduction by the de-noising methods (Rajesh et al., 2017). SER was computed from the power spectral density at the frequency range of P wave (5-30 Hz), HBS (75-100 Hz) and noise (150-200 Hz). These discrete frequency regimes were fixed based on the frequency spectrum of signal-averaged waveforms of PR intervals showing the frequency components of P, HBS and noise for all the subjects discussed in this study. The frequency regime 150-200 Hz always contained noise components and were found to be less likely to contain signal features of interest in this study. Hence, the exclusive contribution of the methods used in this study to the reduction of random noise was evaluated at this noise frequency range.

| Statistical analysis on inter-observer variability of HV intervals
By treating the HV interval values determined already by us as a reference set, the values calculated by the two independent electrophysiologists were compared by computing Spearman's rank correlation coefficient (ρ) which is widely used to assess inter-observer variability in medical research (Patrick et al., 2018) as given by this following equation.
where d is the difference between the two ranks of HV intervals measured by each observer, and n is the number of HV interval measurements.
In addition, Bland-Altman analysis (Giavarina, 2015) was also performed on the differences of HV intervals and those calculated by the two independent observers to evaluate their level of agreements.

| Results on the application of the proposed method on cardiac time series
The stage-by-stage outputs of the proposed analysis method comprising IDWT and SSP are shown in Figure 3. A discernible HBS seen in the averaged PR interval trace was taken as a template for IDWT. HBS in MCG (defined with its onset) was identified based on the criteria followed in an earlier work with respect to multichannel information . The HV interval measured as 110 ms for this subject in MCG was comparable to that determined in intracardiac recordings as 120 ms.   Three representative noisy raw epochs of the PR interval of MCG are shown in Figure 5a, and the application of this method which selectively isolated HBS from the raw epochs, is shown in Figure 5b.
It is to be noted in Figure 5c, the conventional method of averaging these epochs could retrieve HBS (red trace in Figure 5c   Although MCG data of 21 subjects were analyzed in this study, single beat extraction of HBS could be achieved only in 19/21 subjects using the proposed method due to very high level of atrial ramps in the other two subjects. It could be noticed from the Table   that SNR of the as-presented signal-averaged MCG trace was already adequate (~15-20 dB) for HBS, but with significant dominance of atrial activity seen from SNR at P wave. This was drastically reduced by the proposed method even in a single cardiac beat. Similarly, the observations of higher values of SER at HBS frequency and low SER values at P wave and at noise frequencies further endorsed the fact the prominence of HBS alone was facilitated in the PR interval by the proposed IDWT-SSP method. Figure 9 features the application of the proposed method to MCG measured on a healthy subject during a moderate stress condition, characterized by beat-by-beat autonomic variations in the heart rate. Figure 9a shows heart rate variations plotted against the beat numbers. A gradual increase in the heart rate commensurate with physical stress causing fluctuations from a normal resting state ~80 beats per minute (bpm) to a maximum of 120 bpm and back to ~80 bpm is shown. The beat numbers (epochs) were classified into six analysis bins each of 30 beats in length. Signal-averaged epochs (of PR interval) in each analysis bin are shown in Figure 9b in the order from top to bottom (representing bins 1 to 6). The averaged traces show an inverted P wave with short PR segments (the portion of the PR interval after the P wave offset) with bumps attributed to HBS identified based on criteria described earlier with reference to signal registration in multiple channels .

| Results on the statistical analysis on HV interval
The Spearman's rank correlation coefficient (ρ) on the inter-observer variability of HV intervals was computed using Equation (5)  The corresponding invasive EPS recordings featuring HBS spike (indicated by arrows) marked with the measured HV interval of 120 ms.
with p value <.00065 for both the cases. Figure 10 shows the plot on Bland-Altman analysis performed on the differences of HV intervals by the two independent observers. As could be seen from the figure, the mean differences are 1.7 and 1.43 ms, respectively, for observer 1 and 2. The upper and lower limits of agreement at ±1.96 times standard deviation are calculated as −9 to 12 ms for observer 1 and −7 to 10 ms for observer 2. Hence, the overall mean of the limits of agreement is −8 to 11 ms, and these figures indicate that there is no significant difference in HV values across the observers, validating the analysis presented here on further extracting the HBS in a beatby-beat manner.

| DISCUSS ION
The present work reported the surface registration of the HBS using of the His bundle as well [Tavarozzi et al., 2002]). The efficiency of this approach was primarily dependent on the prominence of the His bundle spike or bump in the template beat (averaged beat), which was used for fixing the threshold parameters for wavelet thresholding. Several factors affected the notability of His bundle bump such as, the closeness of the bump to the offset of the P wave, level of atrial repolarization, optimal measurement location on the chest and most importantly, the input SNR of the measurement Senthilnathan et al., 2012;Yamada et al., 2003). Our simulation results prescribed an SNR range ~15 dB and an atrial repolarization presented with a magnitude ~1/3rd of the P wave (Type 1 in simulation) to be optimum for recovering the HBS in every beat using this method. This estimation on the adequacy of SNR was consistent with the quantitative assessment performed on the retrieved MCG-HBS of all the subjects in Table 1.
However, this analysis also suggested a scope for further improvement in the input SNR of the signal by pre-processing MCG data before applying the proposed method. The magnitude limit suggested for PT a was also within the expected range considering independent studies on the measurements of atrial repolarization (Fenici & Brisinda, 2007). The maximum correlation of the extracted HBS that could be achieved using the proposed method was ~70%-80%, even in a high SNR. This indicated a possible distortion in the retrieved signal. Nevertheless, the contribution of these distortions for the purpose of visualizing HBS using contours was insignificant.
The conventional use of filters for HBS like that employed in EPS might require extensive evaluation due to the possible imposition of filtering distortions (Agarwal & Gupta, 2013;An & Stylios, 2020) and might mislead the noninvasive identification of HBS. The present work managed the filtering problem by defining the region of interest and appropriately using wavelet thresholding and SSP as demonstrated. Use of projection-based techniques to selectively remove drift signals whose frequencies overlapped with ramps in the PR segment were reported in literature (An & Stylios, 2020) and it justified the choice of SSP for this purpose of removing atrial components.
It should be admitted that the proposed method has some limitations primarily on its application to cases where significant beat-by- the visualization of HBS in inter-beat contour as shown in this work could still be acceptable for needy cases where continuous tracking of HBS is expected to be informative like for example, in studying the influence of pharmacological interventions (Su et al., 2020) and follow-up measurements where a conventional EPS recording involving fluoroscopic radiation may not be advisable (Muresan et al., 2019). There are ample scopes for refinements on the proposed method with respect to the automation of the IDWT process for fixing thresholds settings and adapting the methodology to encounter arrhythmic beats etc.

| CON CLUS ION
In summary, a new approach has been presented to track and visualize HBS in every beat through inter-beat contour maps with sufficient

ACK N OWLED G M ENTS
We are thankful to all those who had helped us during compilation of this work as an article.

CO N FLI C T O F I NTER E S T S TATEM ENT
There is no conflict of interest for any of the author of this manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this analysis is available with the corresponding author and could be obtained upon reasonable request for academic purposes, subject to the approval of the institutions of the authors.

E TH I C S S TATEM ENT
The research reported in this work was approved by the Institutional