Software algorithm and hardware design for real-time implementation of new spectral estimator

Background Real-time spectral analyzers can be difficult to implement for PC computer-based systems because of the potential for high computational cost, and algorithm complexity. In this work a new spectral estimator (NSE) is developed for real-time analysis, and compared with the discrete Fourier transform (DFT). Method Clinical data in the form of 216 fractionated atrial electrogram sequences were used as inputs. The sample rate for acquisition was 977 Hz, or approximately 1 millisecond between digital samples. Real-time NSE power spectra were generated for 16,384 consecutive data points. The same data sequences were used for spectral calculation using a radix-2 implementation of the DFT. The NSE algorithm was also developed for implementation as a real-time spectral analyzer electronic circuit board. Results The average interval for a single real-time spectral calculation in software was 3.29 μs for NSE versus 504.5 μs for DFT. Thus for real-time spectral analysis, the NSE algorithm is approximately 150× faster than the DFT. Over a 1 millisecond sampling period, the NSE algorithm had the capability to spectrally analyze a maximum of 303 data channels, while the DFT algorithm could only analyze a single channel. Moreover, for the 8 second sequences, the NSE spectral resolution in the 3-12 Hz range was 0.037 Hz while the DFT spectral resolution was only 0.122 Hz. The NSE was also found to be implementable as a standalone spectral analyzer board using approximately 26 integrated circuits at a cost of approximately $500. The software files used for analysis are included as a supplement, please see the Additional files 1 and 2. Conclusions The NSE real-time algorithm has low computational cost and complexity, and is implementable in both software and hardware for 1 millisecond updates of multichannel spectra. The algorithm may be helpful to guide radiofrequency catheter ablation in real time.


Introduction
The discrete Fourier transform (DFT) is used ubiquitously for power spectral estimation of biomedical data. However, it does not readily lend itself to real-time analysis. When successive snapshots of the signal spectrum are desired [1], a real-time spectral estimator is needed. Early work found that it was possible to implement recursive and nonrecursive procedures to update the DFT estimate when an analysis window of N values (i.e., an ordered set, or block) is staggered by a moving size of M sample points, with both N and M being powers of two [1]. Using a complex split-radix implementation, it is possible to improve efficiency by 2-3 × for real-time update as compared with the radix-2 implementation [2]. Other algorithms have been developed that specifically target the implementation of the two-dimensional DFT, with reduction of up to 50% in the computation time being observed [3]. Similarly, a split-vector-radix implementation for two-dimensional DFT analysis was found to reduce the computation time to half or less as compared to the MATLAB implementation [4]. However, this latter algorithm works best for real-time analysis only when the moving size M is large, i.e., there is the need to skip samples, or when the analysis window N is shortened, in order for real-time DFT update to be feasible. Yet increasing M or decreasing N will result, respectively, in a diminished time and frequency resolution of the real-time analysis. Ideally, a sliding window moving size of M = 1 sample point, (i.e., spectral update for every new input sample point) should be used so that the real-time spectral analysis can capture all frequency changes. The optimal window length for atrial fibrillation (AF) frequency analysis was shown to be N = 8192 sample points in a prior study [5].
The development of a real-time spectral analyzer having moving size M = 1, window length N = 8192 discrete sample points, and low computational overhead would potentially be useful for mapping of AF electrograms in order to detect and localize arrhythmogenic regions for radiofrequency ablation. Although current frequency analysis methods hold some promise for this purpose [6,7], due to a lack of resolution, the algorithms are unable to capture transient periodicity, i.e., periodic signal components that are present for only brief intervals, which are likely related to independent, intermittent drivers of the arrhythmia. If an algorithm for real-time spectral update were to be devised without reducing the time or frequency resolution, it could be helpful to detect transient electrophysiologic events, as well as to improve analysis of the detail and trends in sequentially-acquired atrial electrical activation data. In this study both a real-time software algorithm and a hardware design (prototypical hardware electronics board block diagram) for power spectral analysis are described based on a new spectral estimator (NSE) [8][9][10][11] which is applied to fractionated atrial electrogram data.

Clinical data acquisition
Atrial electrograms were recorded from patients referred to the Columbia University Medical Center cardiac electrophysiology laboratory for catheter ablation of the AF substrate. These recordings were obtained prospectively as approved by the Institutional Review Board, and analyzed retrospectively for this study. Nine patients had documented clinical paroxysmal AF. Ten other patients had longstanding persistent AF that did not terminate for several months or years prior to catheter mapping and ablation. The atrial mapping procedure was done using a NaviStar ThermoCool catheter, 7.5 F, with 3.5 mm tip and a 2 mm spacing between the bipoles of the distal ablation electrode (Biosense-Webster Inc, Diamond Bar, CA, USA). The electrogram signals were acquired using a General Electric CardioLab system (GE Healthcare, Waukesha, WI), and filtered at acquisition from 30 to 500 Hz with a bandpass filter (single-pole) to remove baseline drift and high frequency noise. The filtered signals were then sampled at 977 Hz (i.e., approximately 1 millisecond intervals between samples) and the digital data was stored to PC-type computer disk. Although the low pass corner frequency was slightly above the Nyquist value, negligible electrophysiologic signal energy resides above this range [9]. The data was extracted without patient identifiers to a USB drive for subsequent analysis. As in previous studies, to standardize the morphological characteristics, all CFAE were preprocessed to mean zero and unity variance (average signal level = 0 volts, standard deviation and variance = 1).

The NSE algorithm
For computer implementation of the NSE, a software algorithm was developed. The NSE power spectrum is formed from the ensemble averages of signal segments. When successive signal segments of length w are added together, if the segments are correlated, then the sum will reinforce the individual components, which would indicate periodicity at an interval w. The ensemble mean vector e À w of dimension w × 1 is defined as the average of n successive segments of an N × 1 dimensional signal À x N , where À x N is normalized to mean zero and unity variance prior to analysis. Each segment À x w;i of this signal, of dimension w × 1, is used for averaging: where underscore indicates a vector quantity, the first subscript denotes the vector length, and: with x À w,i being signal segments of length w, from i = 1 to n, with the number of summations for averaging, n, given by: and int is the integer function (decimal is rounded down). For each ensemble mean vector, the sum of squares of all elements divided by the vector length is defined as the ensemble power. This can be expressed as: The NSE power spectrum is plotted as: where P wRMS is the root mean square power, and the scaling term √n is used to counter the falloff of the baseline level by 1/√n per number of summations n used for signal averaging [9]. From Eq.'s 1 and 4: where i is the segment number from 1 to n. Thus from Eq.'s 5 and 6a: As in prior studies [8][9][10][11], Eq. 7b was used to implement the offline version of the NSE power spectral calculation.
The real-time NSE power spectral calculation was developed so as to estimate the ensemble mean from a moving average. By using a moving average, only one calculation is needed to update e w for each power spectral calculation, rather than n calculations as in Eq. 1. From Eq.'s 4 and 5, the normalized root mean square power is: where the √N term is a constant, and for the N = 8,192 sample point window size used in this study, it has a value of 90.51, while 1/√N = 0.01105. For real-time power spectral updates, a moving average < e w > is substituted for e w in Eq. 8d, thus: For both the offline and real-time NSE power spectral calculation, √n P wRMS is plotted versus frequency f, which is given by: The ensemble power is computed in the standard electrophysiologic frequency range of 3-12 Hz [12,13]. At the 977 Hz sampling rate, this corresponds to a range of w from 325-81 discrete sample points. Thus in the electrophysiologic frequency range, the NSE is a 245 point power spectrum. The signal length N for analysis does not affect the spectral resolution, but it does affect the number of summations for averaging used to calculate e w and P w . Greater number of summations improves cancellation of random components while forming the ensemble mean; however, the DF may not be stable over long intervals [14,15]. Over short intervals, the NSE estimate will be improved when the signal is composed mostly of periodic components and less random noise. Although the DFT does not depend on a summation of signal components, the DFT frequency resolution is inversely proportional to the window length; thus short intervals have low resolution.

Software implementation
A computer algorithm to implement the offline NSE power spectrum calculation is shown below (called spectral_estimate). In the source code, period w ranges from n1 = 325 to n0 = 81 sample points (3)(4)(5)(6)(7)(8)(9)(10)(11)(12). The calculation window is n2 = 8,192 sample points, which is approximately 8 seconds when sampled at a rate r = 977 Hz. The number of signals (multichannels) analyzed n3 = 216. For calculation, the input data array, with each patient data being received from a separate file, is 'inp' (lines 4-6), the ensemble mean is stored in array 'en' (lines 7-15), and the generated spectrum is stored in array 's' (lines [16][17][18][19]. Note that line 17 in the software code arises from Eq. 7b. For simplicity, the spectrum of a single electrogram was written to disk, arbitrarily taken as channel 4 (line 18). For faster run time, lines 8-12 calculate the ensemble means from w = 163 to w = 325. Then in lines 13-15, the ensemble means from w = 162 to w = 81 are computed by simply adding the two halves of vector 2w that is calculated by lines 8-12, thus: is determined in lines 13-15, which reduces redundancy in calculation. The program spectral_estimate is used for offline (non real-time) computations, i.e., it calculates a single spectrum for each patient record using sample points 1-8192, and it has been used in prior studies [8][9][10][11]. The offline time for computation of a single NSE spectrum for 216 patient records was compared with the offline radix-2 implementation of the Fast Fourier Transform (FFT) [16]. The time for spectral estimation, lines 7-20 was determined using the 'date_and_time' Fortran function, which was inserted between lines 6-7 and between lines 20-21 in the actual source code. Using the same 'date_and_time' function, the run time for spectral estimation using the offline radix-2 FFT implementation was also determined. Since slight temporal changes in processor speed can occur, the mean and standard deviation in runtime over five trials was calculated for both the NSE and FFT methods.
The NSE and FFT spectral estimators were then implemented for real-time analysis. An update was computed once every sample point (moving size M = 1). In Figure 1, a fractionated electrogram is shown. Calculations of the ensemble mean were done using a moving average. Early in the sequence, for example at discrete sample points 2 and 12 (show as dashed vertical lines in the top graph), the moving average was based upon the few available prior data points and was a rough estimate. Later in the sequence, for example at sample point 8192 (dashed vertical line in lower graph), the moving average was based upon many prior data points. The software code for real-time NSE analysis is shown following (spectral_estimate_real_time). The declaration lines are virtually the same as for the offline program (spectral_estimate). In the real-time program, most variables are vectors whose elements represent a range in w from 81 to 325. Declaring length w as an integer was found to considerably reduce computation time as compared to a floating point declaration. Lines 5-8 compute constants c1 and c2 for a moving average, or low pass filter, which is used iteratively for a subsequent calculation in the coding sequence. Lines 9-11 input the data, the same as lines 4-6 in the offline program (spectral_estimate). However the ensemble mean calculation (lines 12-18) is slightly different as compared with the spectral_estimate program. The actual calculations are done in lines 15-17. The index 'ind' is first calculated (line 15), which is used to point to a particular element within each vector that is used for calculation. The index is incremented by 1, and after it reaches a value of w, it wraps around to a value of 1. On line 15 the prior value of the ensemble mean at the index is also stored in a buffer (ee). The ensemble mean is then updated using the moving average (line 16), which includes the constants c1 and c2 calculated in lines 5-8. The realtime ensemble mean calculation can be expressed as: where: with n, and therefore c1 and c2, being dependent upon the value of w (Eq. 3). The power of the ensemble mean vector of length w is then updated in line 17. The ensemble power can be expressed as: where e is the scalar ensemble mean, index i ranges from 1 to w, and the divide by w is accounted for when the spectral point is calculated subsequently in the code. At any particular index value i = 1 to w, the prior squared ensemble mean e i − w ⋅ e i − w , which was calculated w sample points previously in time, is removed, and it is replaced with the newly calculated value e i ⋅ e i . Thus at each index i, the power is updated as: where the index i is incremented by 1 with each new input sample point. For each new sample point, all P w are updated using Eq. 16. All spectral points S w are then calculated as the square roots of P w , scaled by n(w)/√N as in Eq. 8d, with 1/√N symbolized by xn in lines 1 and 17 of the spectral_estimate_real_time software program. It was found that implementation using Eq. 8d rather than 8c, so that there is no actual 'divide by' in the computation, considerably reduces the runtime. One complete spectrum, arbitrarily selected as that for data record 4, is stored for subsequent analysis (lines [19][20][21][22]. The frequency variable, which is the x-axis parameter for graphing, is also stored (r / w, line 21).
Note in the spectral_estimate_real_time program that 16,384 spectra are calculated for each of 216 patient records (lines [12][13]. This program is designed as a real-time implementation, updating the spectrum with each new input sample point. However, it will only be justified as being an actual real-time algorithm if the update calculations for all channels can be done at a speed faster than the discrete sampling interval, which is approximately 1 millisecond. As for the spectral_estimate program, the time for spectral estimation using spectral_estimate_real_time was determined by inserting the 'date_and_time' Fortran library function at the appropriate locations in the code (between lines 11 and 12, and after line 23).
The radix-2 FFT program was also implemented for real time [16]. For ease of comparison, all FFT calculations were repeated for each sliding analysis window rather than using any of the shortcuts suggested by other authors [1,2]. Because the repetition of FFT calculations is time-consuming, it was run for only 512 input sample points. Again inserting the Fortran 'date_and_time' function at appropriate locations in the coding sequence, the time for FFT calculation over 512 sliding windows was determined, which was then scaled by 32 × for comparison with the NSE implementation in which spectra were calculated for 16,384 sliding windows. Since slight temporal changes in processor speed can occur, the mean and standard deviation in spectral computation time over five trials was determined for both the NSE and DFT real-time methods.
As an additional test of the real-time NSE implementation, spectral parameter values at sample point 8,192 were computed for the 216 data files, and the results were compared to values computed using the offline NSE implementation. It was shown previously that significant differences exist in the dominant amplitude (DA), dominant frequency (DF), and mean and standard deviation in spectral profile (MP and SP, respectively) of paroxysmal versus persistent AF fractionated electrogram recordings when using the NSE and DFT offline implementations [17]. The DA is defined as the amplitude and the DF is the frequency of the largest fundamental peak in the power  Figure 1 Windows used for analysis. A moving average was used, which was calculated based upon prior discrete sample points in the sequence.
spectrum over the physiologic range of 3-12 Hz [17]. The MP and SP are, respectively, the mean and standard deviation of the spectral profile after the magnitude has been normalized to a range between 0 and 1. Unlike the regularity and organizational indices, spectral parameters which require guestimates of the width of the DF spectral peak and of its harmonics [18,19], no such guestimates are required to calculate the DA, MP, and SP parameters. These new parameters have also been incorporated in the study of another group [20]. If real-time implementation could produce similar means and standard deviations in these spectral parameters, and similar significant differences, it would provide evidence that the real-time NSE implementation is similar to the offline implementation. The NSE and FFT source code was compiled using an Intel Visual Fortran Compiler (Standard Edition for Windows, Version 9, Intel Corporation, Santa Clara CA, 2005). All computer calculations were done using a PC-type laptop having the following specifications: The software code and data files used for analysis in this study are included as a supplement (please see Additional files 1 and 2).

Hardware implementation
Lines 13-23 of the spectral_estimate_real_time program were then implemented in schematic form in hardware. For ease of implementation, CMOS process integrated circuits were selected to do the calculations shown in the software code. The integrated circuits were also selected for their utility and high performance. The goal was to implement the calculations on a small prototype electronics board. A requirement for implementation was that the electronics board would operate in standalone mode, without any computer control. For simplicity, the delays needed to ensure system stability and proper handshaking were not shown, but can be readily implemented with a 555 timer circuit to provide a delay for latching data. An approximate cost for parts was determined based on the schematic hardware implementation.

Overview
The DFT spectral resolution in this study is given by: Thus in the range 3-12 Hz, the DFT spectrum will consist of 9/0.122 ≈ 74 spectral points. The mean NSE spectral resolution for this study is approximated by [11]: The NSE spectrum in the range 3-12 Hz consists of 245 spectral points, and from Eq. 18, the mean resolution is 0.037 Hz. The NSE spectral resolution is therefore improved by more than 3 × over the DFT spectral resolution for this study.
Examples of power spectra for the offline NSE and DFT implementations are depicted in Figure 2. The spectra were generated from a sequence acquired from the antrum of the right inferior pulmonary vein in a persistent AF patient using the distal ablation catheter bipolar electrode. The DF occurs at approximately 7.2 Hz, and a secondary peak, perhaps caused by an independent generator of electrical activity in the left atrium, is present at approximately 9.7 Hz. The individual spectral points are shown by solid circles. The DFT resolution is fixed at 0.122 Hz and consists of approximately 74 spectral points. The NSE resolution varies with frequency and consists of 245 spectral points. This makes the NSE spectral peaks better defined and sharper.
Graphs of real-time NSE spectral estimates were constructed for a patient with longstanding persistent AF. The spectral estimates were generated from fractionated electrogram records. In Figure 3, instances are shown of the spectral estimate using the original NSE algorithm spectral_estimate (panels A and C) and the real-time NSE algorithm spectral_estimate_real_time (panels B and D). Panels A and B show a recording from the left superior pulmonary vein antrum in a patient with persistent AF, and panels C and D were recorded from the left inferior pulmonary vein antrum in the same patient. For spectral_estimate, the result of analysis of the static window with data points [1][2][3][4][5][6][7][8]192 is shown. For spectral_estimate_real_time, the result of analysis at sample point 8,192 is shown (since the analysis windows are moving averages). There are only subtle differences in the result as can be seen in panels E and F, where the spectra are overlapped after scaling to mean zero and unity variance (offline-black, real-time-gray). The root mean square error difference is 0.23 ± 0.21 mV 2 in panel E and 0.34 ± 0.26 mV 2 in panel F.
The statistical results for paroxysmal versus persistent AF spectral parameters are shown in Table 1. The means and standard deviations for spectral parameters in persistent versus paroxysmal AF patients are noted. For both the offline and real-time algorithms,  the DA parameter is of greater magnitude in persistents versus paroxysmals and the difference is highly significant (p < 0.001). Similarly, for both offline and real-time algorithm, the DF parameter is higher in persistents versus paroxysmals and the difference is highly significant (p < 0.001). Furthermore, for both offline and real-time algorithm, the MP and SP parameters are of lesser magnitude in persistents versus paroxysmals and the difference is highly significant (p < 0.001). The measurements suggest that electrical activity is less organized in paroxysmal AF as compared with persistent AF [17], though ablation of high DF sites has been shown to more likely result in arrhythmia termination for paroxysmal AF [7]. The means and standard deviations of all parameters have similar values for the offline versus real-time algorithm (Table 1). There are no significant differences in the means of offline versus real-time parameters for persistent AF (P ≥ 0.201) and for paroxysmal AF (P ≥ 0.302). Therefore, the real-time implementation provides a similar spectral estimate to the offline implementation at a particular time epoch.
To illustrate how the time-varying spectral content can be detected using the realtime algorithm, Figures 4 and 5 show spectra from successive points in time for data acquired at the posterior left atrial free wall in a patient with persistent AF. The split peaks at the DF change slightly over short time intervals of 50-100 milliseconds between frames as in Figure 4. Over longer intervals of 1000 milliseconds, major changes in the spectra occur ( Figure 5). The split peaks at the DF (k = 12,000 and 13,000 sample points) become single peaks (k = 14,000 and 15,000 sample points). Over successive frames in Figure 5, the background level becomes somewhat reduced in magnitude, particularly at lower frequencies for the spectra acquired at k = 15,000 sample points. Thus temporal changes in spectral details are evident over both short and longer intervals using the real-time implementation of NSE.

Summary software statistics
Computation of single (i.e., offline) 8,192-point spectra for 216 patient records completed in 0.38 ± 0.01 seconds using NSE and in 0.12 ± 0.01 seconds using the radix-2 FFT implementation. Thus the offline FFT calculation is more than 3 × faster than the NSE. This result can be explained by the increased number of calculations required by NSE for a single spectral estimate as compared to FFT. For real-time analysis, 16,384 spectra from 216 patient records (3,538,944 spectra total) were calculated in 11.63 ± 0.14 seconds for NSE, while 512 spectra for 216 patient records (110,592 spectra in total) were calculated in 55.79 ± 0.26 seconds for the radix-2 FFT implementation. Scaling by 32×, 3,538,944 spectra would be calculated in 1785.28 ± 8.32 seconds for FFT. This represents a speed advantage of 153:1 for NSE over FFT when the real-time analysis algorithm is implemented. Dividing each time by 3,538,944, the average time for a single real-time spectral calculation was 3.29 μs for NSE versus 504.5 μs for FFT. During a 1 millisecond sampling period, the NSE algorithm therefore has the capability to spectrally analyze a maximum of 303.95 data channels, although only 216 patient data sequences were actually analyzed for this study. By comparison, the FFT algorithm could only analyze a maximum of 1.98 channels, which for purposes of calculation means that the FFT can only analyze a single recording sequence in real time within the 1 millisecond sampling interval.

Hardware implementation
The schematic for hardware implementation based upon the NSE real-time software algorithm (spectral_estimate_real_time) is shown in Figure 6. The hardware consists of a mix of analog and digital components, and for simplicity it is shown as sampling the Figure 5 NSE spectra using the real-time software algorithm for larger change in analysis window location of 1000 milliseconds. Changes from split peak to single peak can be observed over time.
data at exactly 1.0 kHz (1 millisecond intervals) using a 1 kHz clock. To update parameters from w = 81-325 corresponding to a 12-3 Hz frequency range, which consists of 245 values, a 245 kHz clock was used. For ease of implementation, an off-the-shelf 250 kHz oscillator could just as readily be used. The data stream (signal) is tracked and then held with an edge of the 1 kHz clock (labeled clk 1). The analog signal is then valid (hold) until the next 1 kHz pulse. The edge of the 1 kHz clock waveform is also used to reset a counter (wcounter). This counter provides the current value of segment length w being analyzed, and it counts in increments of 1 from 81 to 325 on each pulse from the 245 kHz clock. The wcounter output (labeled w) is used as addressing information for the memory chips that are present on the electronics board, as well as for data input elsewhere in the circuitry. In the schematic diagram, addressing information is indicated by dotted lines, while data is noted by solid lines.
In Table 2, the characteristics of each memory chip on the board are shown. Stored in EMEM are the ensemble means for each segment length w from 81-325. Each ensemble mean is a vector of length w. Therefore the addressing information sent to the EMEM chip, when a newly calculated ensemble mean value is to be stored, must specify not only the segment length w (received from wcounter), but also the current index number i for that w that is being addressed (received as 'index' from the index memory chip-IMEM). Each time a new input sample point is received, the index number for each w is incremented by 1, with wraparound. As an example, for w = 100, the index begins at 1 and is increased to 2, 3, 4, …, 100 for each new input sample point (i.e., once every millisecond). When the index reaches 100 it must then wraparound to 1. This reset is done by use of the counter and comparator that are associated with IMEM. The index value contained in IMEM for segment length w is presented to the counter (icounter). Using a delayed pulse from the 245 kHz clock, the icounter is then incremented by 1. The icounter output after incrementation is compared to the current value of w (received from wcounter). If the icounter output is greater than w (i.e., its value is w + 1), it is reset to 1. The icounter output is then stored in IMEM at the address value w that is currently being accessed.
The digital ensemble mean value E(k) contained within EMEM for a particular w and index is output as E(k-1) to a digital-to-analog converter (DAC). The DAC output is multiplied by the constant c1. The new input sample point (analog signal in the diagram) is multiplied by the constant c2. The constants c1 and c2 are used to weight the moving average of the ensemble mean (line 16 of spectral_estimate_real_time), according to the following equation: where, as computed in lines 6-7 of spectral_estimate_real_time: with N being the window length, 8,192 sample points, 'int' is the integer function, and: Eq.'s 19-22 are analogous to Eq.'s 3 and 12-14 for the software implementation. The values of c1(w) and c2(w) are constant, and it is supposed that they can be stored in C1MEM and C2MEM prior to real-time spectral analysis, using a separate device to write to the chips, to maintain the standalone quality of the electronic board design. The two products on the right-hand-side in Eq. 19 are then summed. This is illustrated in Figure 6, with the result E(k) being input to an analog-to-digital converter (ADC). The resulting digital value of E(k) replaces E(k-1) in EMEM for the w and index currently being pointed to by the addressing information. The moving average filter output is also squared, which is then used as one of the elements inputted to a summing circuit to revise the value of the ensemble power P. The update of P at discrete time k is given by the equation: where Eq. 23 is analogous to Eq. 16 used for software calculation. The quantity P(k-1) is accessed from the ensemble power memory chip (PMEM). The quantity E(k-1) 2 is obtained by squaring the output of the EMEM DAC. The output of this summing circuit is the analog power P at time k. This value is digitized using an ADC and stored as the new value of P at w in PMEM. The square root of P(k) is also obtained, which is then scaled according to Eq. 8c or 8d (for simplicity shown using Eq. 8c, √N / w), where √N = 90.51, and a DAC is used to convert w to analog form (for simplicity, not shown in Figure 6). The result is digitized and stored in the spectral memory (SMEM). SMEM contains 245 digital values representing the current spectral content from w = 81 to 325 (12-3 Hz). The SMEM digital output is then sent to a display. For simplicity, the spectral display device is not shown. To prevent the need to tie the spectral analyzer board to a computer, the output of SMEM could be sent, for example, to a dot matrix display with appropriate circuitry.
In total according to Figure 6, approximately 26 integrated circuit chips, plus delay, display circuitry, and power supply (not shown) can be used to implement the realtime NSE algorithm on a prototyping electronics board. The circuitry must run at nearly 250 kHz; thus the settling time for the integrated circuits should be less than 4 microseconds, which necessitates the purchase of high quality components. The DAC and ADC integrated circuits should be parallel input and output, respectively, for faster throughput; thus single rather than dual or quad packages should be used for these devices. For reduced cost, an 8-bit digital resolution can be selected. For memory access, latching, and counting, sufficient delay will be required to allow inputs to each integrated circuit chip in Figure 6 to settle. This can be done for example using 555 timer circuits, and the associated resistors and capacitors, which for simplicity are not shown in the schematic.
Since the fast clock runs at approximately 250 kHz, the generation of each spectral point must be completed within four microseconds. As shown in Table 3, all components of Figure 6 except the display unit would be purchasable in high performance versions for approximately $500 at the current market pricing for individual components. The sample rate with real-time spectral update of the circuitry, outlined in Figure 6, could be improved to a maximum value of approximately 4 kHz, using a 1 MHz fast clock (this is the maximum run speed according to the longest settling times of 1 μs listed in Table 3). [μs is the microsecond symbol].

Summary
In this study the NSE method was implemented for real-time analysis as a software algorithm and as a block diagram for a prototypical hardware electronics board. To be implemented in real time, the spectral estimate is updated within the time interval needed for the data stream to shift by the analysis window moving size of 1 sample point (M = 1), which for the application presented was 1 millisecond. The data stream consisted of sequences of retrospectively analyzed fractionated atrial electrograms from AF patients. The NSE algorithm was found to be implementable in real time using a few lines of software code. The mean time for calculation of one power spectrum was 3.29 μs for NSE versus 504.5 μs for the conventional radix-2 FFT implementation. Thus for real-time spectral analysis, the NSE algorithm was found to be over 150 × faster than the conventional FFT. Based on these values, for a 1 millisecond sampling period, the NSE algorithm has the potential to spectrally analyze over 300 data channels while the DFT algorithm can only analyze a single channel. Thus for NSE, although 216 sequences were actually analyzed, the number of sequences could be increased to 300 while maintaining real-time calculation within the 1 millisecond window moving size. Whereas for the conventional radix-2 FFT implementation, only one channel can be analyzed when the moving size M = 1. The rapid speed of NSE for real-time analysis is due to the low computational overhead in calculating the update as compared with the FFT recalculation. The NSE was found to be implementable in hardware using approximately 26 integrated circuits at a cost of approximately $500. No computer controller or digital signal processor would be needed to run the hardware implementation-it was developed as a standalone spectral analyzer board. Based on the settling times of the mixed-signal circuitry in Figure 6, the sampling rate could possibly be increased from 1 kHz to 4 kHz while maintaining a real-time spectral output. To analyze more channels simultaneously, the circuitry of Figure 6 would be repeated.

Optimized implementations
The FFT implementation used for this study was not optimized, but was selected because the software code is in the public domain. To estimate the approximate speed advantage if an optimized FFT implementation were to be used, benchmark data can be compared for the Numerical Recipes implementation versus a 'Fastest Fourier Transform in the West' (FFTW) implementation termed 'fftw3' [21]. FFTW is a structured library of optimized C code blocks termed codelets, and is portable to any platform having a C compiler. It is used for computing the FFT in one or more dimensions, with any input size, for either real or complex data. Over a variety of platforms, FFTW's performance has been shown to be improved over most other publicly available FFT software, and is competitive with commercial software [21]. Benchmark data was available from two CPUs with performance that is comparable to the CPU of the laptop computer used in this study-the 2 GHz PowerPC (970 (G5), 32 bit mode, gcc-3.4) and the 3.0 GHz Intel Core Duo (32-bit mode). For the Power PC implementation, benchmark data was available for computing the 1-dimensional FFT, single precision, with complex numbers. The Numerical Recipes in Fortran version, which was used in our study, was benchmarked at 600 megaflops, while the fftw3 version ran at 6000 megaflops. For the Intel Core Duo implementation, benchmark data was available to compute the 1-dimensional FFT, single precision, with real numbers. Numerical Recipes in Fortran was benchmarked at 1000 megaflops, while the fftw3 version ran at 8400 megaflops. Thus using a highly optimized FFT set of library functions and compiler, it is possible to obtain an approximate 10 × speed advantage as compared with the Numerical Recipes FFT implementation. Supposing that a 2 × speed advantage can be gained by including shortcuts to update the FFT in real time, it may therefore be possible to obtain overall as much as a 20 × speed increase in running a real-time FFT using optimized hardware and software as compared with the standard FFT used for our study. However, even if the best optimization protocols were implemented to speed up the FFT by approximately 20 ×, and no steps were taken to optimize the NSE algorithm, NSE would still be expected to run approximately 150/20~7.5 × faster than the FFT.

Real-time spectral analysis
A comparison of four software algorithms for real-time spectral analysis of AF signals was recently studied [22]. The estimators used were the FFT, Blackman-Tukey (BT), Autoregressive (AR), and Multiple Signal Classification (MUSIC) methods. The window size was 4096 discrete sample points. Calculations were made for a number of segments and averaged. The investigators found the time for spectral calculation for 3000 signals. Scaling for 216 signals, time for calculation for FFT, BT, AR, and MUSIC was respectively, 1.6, 24.5, 14.9, and 781.3 seconds. Supposing the spectra of 10 segments were averaged, then the FFT computation time for one segment would be 0.16 seconds, which is comparable to our value for static implementation of 0.12 seconds. The study therefore suggests that the FFT is the fastest spectral estimator of the group, and that the runtime for the FFT is similar to that found in our study.
Hardware implementations of real-time spectral analyzers can be purchased either in the form of an electronics board which is inserted onto a computer backplane bus, or as a standalone mainframe unit. Electronic boards are less expensive to purchase but have limited portability because they require a dedicated computer for control of the board. Real-time standalone units can be quite expensive. Models such as the RSA5000 Spectrum Analyzer (Tektronix Inc., Beaverton, OR) can be used to view real-time spectra in the frequency range from DC to gigahertz, but the cost is $25,000 or more. The prototypical spectral analyzer described in our study, when sampling is at a projected maximum of approximately 4 kHz, would have a power spectrum range from DC to 2 kHz. It is implementable without the need for interfacing with a computer and would be a small, highly portable device. Although the bandwidth is limited, the range of frequencies is more than sufficient to analyze most electrophysiologic signals. Another use, for example, would be to display the frequency content of music or voice in real time.
If the real-time NSE is implemented in software it would be possible to improve the calculation time, for example, by using vectorized packages in Visual C++, such as with the OptiVec shareware package (OptiCode-Dr Martin Sander Software Dev) which was written for 64-bit platforms including Windows 7 and 8. The OptiVec is a set of assembler-written functions for both floating-point and integer data types that perform vectorized and matrix operations. As suggested by the documentation, runtime speeds can be improved by 2-3 ×. Use of such a program could therefore enable a substantial increase in the number of data channels that are spectrally analyzed by the NSE method in real-time at 1-4 kHz sampling, and/or to increase the sampling rate. Higher sampling rates will result in a greater temporal resolution, and can also be used to analyze higher frequency components, which could be useful for other applications.

Clinical correlates
Analysis of multichannel electrogram data via spectral estimation has been shown useful to distinguish patients with paroxysmal versus persistent atrial fibrillation based upon frequency gradients [7]. Hence acquisition and real-time analysis of multichannel data is potentially assistive to characterize and localize the substrate changes that occur during atrial fibrillation. It may also be helpful to detect optimal ablation sites to prevent recurrence of arrhythmia [23,24]. In a prospective setting, rather than 216 different patient sequences as was done in this study, the data would be obtained from a multichannel electrode, for example using a noncontact [25] or basket [26] electrode. As the number of multichannel electrode recordings available in these devices is increased, the possibility of analysis of all of the channels in real time becomes more remote, yet would be necessary for optimally targeting arrhythmogenic regions for ablation. Thus, the implementation of a real-time spectral estimator for multichannel data is potentially important to improve clinical outcome.
The dominant frequency or DF is a useful parameter for catheter ablation of arrhythmogenic sites in AF patients [7]. When real-time spectral analysis is used to determine the DF during electrophysiologic mapping in patients with ongoing AF, it is possible to predict the long-term maintenance of sinus rhythm [27]. Real-time spectral analysis has also been used to show a strong correlation of the DF calculated from both unipolar and bipolar surface electrograms, to optical mapping data in ovine hearts [28]. The findings of that study suggest that estimation of activation rate via real-time spectral analysis of surface electrograms and DF calculation, correlates well with optical mapping measurement of the actual electrical activation rate. Thus, real-time spectral analysis of surface electrograms, as are acquired during clinical electrophysiologic study, provides physiologically relevant parameters of electrical activation that can be used for targeting arrhythmogenic regions.
Although our study was limited to bipolar electrogram analysis, the morphology of unipolar AF electrograms is reflective of specific patterns of conduction including broad fibrillation waves, collision of distinct wavelets, slow conduction, and pivoting fibrillation waves [29]. Thus unipolar recordings may be particularly useful to differentiate between AF types and to identify regions with different structural properties.

Limitations
The NSE and DFT real-time algorithms were tested on retrospective data. Implementation on prospective multichannel data would be desirable to determine the speed of the spectral estimators in this setting. Still, in principle there should be no difference in speed by using 216 retrospective patient data sequences versus 216 prospective multichannel data sequences that are acquired simultaneously. A conventional FFT software algorithm rather than an optimized algorithm was used for testing. However, the conventional FFT is likely to be implemented for most analyses because of its simplicity and ready availability. The NSE and FFT algorithms were tested with one software compiler and one computer and operating system. Use of a different software compiler and computer may result in somewhat different computation times. For the hardware implementation, the design was illustrated only as a schematic block diagram. Actual implementation of the hardware design may necessitate changes to the components used.

Conclusions
The NSE algorithm is implementable with low computational cost and complexity in both hardware and software. Real-time 1 millisecond updates of the spectral estimate can be done even when many tens or hundreds of data sequences are being acquired and analyzed during the same discrete time interval, as would be the case when using a multichannel basket electrode. The algorithm can be implemented as a standalone spectral analyzer board at a price of approximately $500 plus the cost of the display unit, without the need to interface with a computer. Although it is possible to improve the efficiency of the DFT real-time update and the compiler by as much as 20 × [1,2], such implementations would not match the improvement gained by using the NSE real-time algorithm, which was found by comparison of software algorithms to be approximately 150 × faster. The sampling rate using the hardware implementation of the NSE could possibly be further increased, being limited only by the settling times of the on-board components as compared with the fast clock speed. The realization of a fast spectral analysis algorithm is potentially helpful increase the number of multichannels analyzed, to characterize spectral transients in biomedical data as they occur, and to update in real time the detailed trends and gradients in the frequency content of biomedical data that may be present over longer sequences. This may be particularly helpful when probing the tissue substrate for anomalous regions, as is the case during electrophysiologic study of AF patients.