SYSTEMS AND METHODS FOR TREATING CARDIAC ARRHYTHMIA

Abstract

Apparatus for monitoring activation in a heart comprises a probe 101, a plurality of electrodes 102 supported on the probe and extending over a detection area of the probe, the detection area being arranged to contact a detection region of the heart. Each of the electrodes 102 is arranged to detect electrical potential at a respective position in the heart during movement of a series of activation wave fronts across the detection region. A processor is arranged to analyse the detected electrical potentials to identify a propagation direction of at least one of the wave fronts, and to generate an output indicative of that direction.

Claims

1. An apparatus for monitoring activation in a heart, the apparatus comprising a probe, a plurality of electrodes supported on the probe and extending over a detection area of the probe, the detection area being arranged to contact a detection region of the heart, wherein each of the electrodes is arranged to detect electrical potential at a respective position in the heart during movement of a series of activation wavefronts across the detection region, a display screen, and a processor configured to: process and analyse the detected electrical potentials to identify a plurality of wavefronts passing the detection region; for each of the plurality of wavefronts, identify a transition period during which the wavefront crosses the detection region and define a series of update intervals occurring during the transition period; for each of the update intervals to determine a propagation direction; and to control the display screen to generate a display indicative of the direction.

2. The apparatus according to claim 1 wherein processor is configured to process the detected electrical potentials by filtering them to remove noise.

3. The apparatus according to claim 2 wherein the processor is configured to record values of the electrical potential at each of the electrodes as a respective unipolar electrogram, and to perform filtering on the unipolar electrograms to filter out noise from the unipolar electrograms.

4. The apparatus according to claim 3 wherein, in filtering the unipolar electrograms, the processor is configured to filter out from the unipolar electrograms noise features which are common to a plurality of the unipolar electrograms.

5. The apparatus according to claim 3 wherein, in filtering the unipolar electrograms, the processor is configured to retain in the unipolar electrograms, features which are unique to each of the unipolar electrograms.

6. The apparatus according to claim 3 wherein, in filtering the unipolar electrograms, the processor is configured to define a multi-electrode model which defines the unipolar signal at each of the electrodes, and to identify a best fit of the model to the unipolar electrograms.

7. The apparatus according to claim 6 wherein the model defines the difference between the two unipolar signals in each of a plurality of pairs of the unipolar signals and the processor is configured to determine from the unipolar electrograms measured values of said differences, and to use the modelled and measured differences in identifying the best fit.

8. The apparatus according to claim 1 wherein the processor is configured to identify at least one of the wavefronts as belonging to one of a plurality of classifications depending on at least one of its shape, its position, and its direction of travel.

9. An apparatus for monitoring activation in a heart, the apparatus comprising a probe, a plurality of electrodes supported on the probe and extending over a detection area of the probe, the detection area being arranged to contact a detection region of the heart, wherein each of the electrodes is arranged to detect an electrical potential at a respective position in the heart during movement of a series of activation wavefronts across the detection region, a display screen, and a processor configured to: record the detected electrical potential at each of the electrodes as a unipolar electrogram; forming an electrogram combination which is a linear combination of at least two of the unipolar electrograms and which is affected more by a change in just one of the two electrical potentials than it is by the same change in both of the two electrical potentials; filter out noise in the unipolar electrograms based on the degree to which the noise is more present in the individual unipolar electrograms than it is in the electrogram combination.

10. The apparatus according to claim 9 wherein the processor is configured to filter out the noise by: defining a model of the electrical potentials at all of the electrodes; defining a voltage combination which is a linear combination of at least two of the electrical potentials in the model; finding a best fit of the model and the voltage combination to the unipolar electrograms and the electrogram combination.

11. The apparatus according to claim 10 wherein the processor is configured to define a measurement matrix recording the unipolar electrograms and the electrogram combination, and an observer matrix which relates the measurement matrix to the model.

12. The apparatus according to claim 11 wherein the observer matrix defines a plurality of voltage combinations each of which is the linear combination of a different pair of the electrical potentials.

13. The apparatus according to claim 10 wherein the voltage combination is a difference between two of the electrical potentials.

14. The apparatus according to claim 11 wherein the processor is configured to define a weighting for each of the unipolar electrograms and the electrogram combination.

15. The apparatus according to claim 14 wherein the weightings are frequency dependent.

16. The apparatus according to claim 9 wherein the processor is configured to analyse each of the unipolar electrograms to identify any tissue activation features that they contain which are indicative of one of the activation wavefronts.

17. The apparatus according to claim 16 wherein the processor is configured to analyse each of the unipolar electrograms using a convolution with an inverse function of time.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0045] FIG. 1 is a diagram of a system according to an embodiment of the invention;

[0046] FIG. 2 is a set of plots showing variation in wave front propagation direction during atrial fibrillation;

[0047] FIG. 3 is a plot of a windowing function used in the system of FIG. 1;

[0048] FIG. 4 is a plot illustrating an algorithm for determining relative activation times used in the system of FIG. 1;

[0049] FIG. 5 shows an image displayed in the system of FIG. 1 showing multiple propagation directions;

[0050] FIG. 6 shows a set of three plots generated by the system of FIG. 1 showing most recent wave front directions at three different times during passage of a further wave front;

[0051] FIG. 7 shows a polar histogram plot of wave front directions generated by the system of FIG. 1;

[0052] FIGS. 8A-8D show the main steps in a method according to an embodiment of the invention;

[0053] FIG. 9 shows examples of wave front patterns detected an analysed by the system of FIG. 1;

[0054] FIG. 10 shows the identification of a series of wavefront positions during transition of a wavefront across a region of the heart; and

[0055] FIG. 11 shows how different types of wavefront can be identified;

[0056] FIG. 12 is a plot of measured unipolar electrograms and corresponding filtered unipolar electrograms;

[0057] FIG. 13 is a plot of the time derivative of a unipolar electrogram of FIG. 12 and the time derivative of the corresponding filtered unipolar electrogram;

[0058] FIG. 14 is plot of three filtered unipolar electrograms similar to those of FIG. 12, showing only a section centred on an activation event; and

[0059] FIG. 15 shows convolutions of the plots of FIG. 12 with a plot of 1/t.

DETAILED DESCRIPTION

[0060] Referring to FIG. 1, a cardiac monitoring system according to an embodiment of the invention comprises a catheter 100 having a probe 101 at one end with a set of electrodes 102 located on it. Each of the electrodes 102 is connected independently through the probe catheter 100 to a computer 104 which is arranged to acquire, store and analyse the voltages detected by the electrodes 102. Specifically the computer 104 comprises a memory 106 and a processor 108. The processor is arranged to sample the voltages detected by the electrodes 102 at a regular sample rate and store the values of the sampled voltages, which form a time series of sample values, in the memory 106, and then to analyse the stored voltage values so as to analyse the activation of the heart in the area contacted by the probe 101. Specifically the data can be analysed to identify focal targets within areas of irregular activation. The processor 108 is arranged to generate from the sampled voltage data, an image data set which it then provides to the display screen 110 which displays an image showing the activation pattern in the heart so that a user can interpret it. The probe 101 can be moved from region to region within the heart to focus attention in the regions where the targets are suspected.

[0061] The catheter 100 may further comprise an ablation tip 114 which is connected to a radio frequency (RF) power source. The ablation tip 114 can therefore be used for ablation of regions of the heart which are found to be sources of atrial fibrillation, whether focal sources or rotors. The catheter may for example be a Smart-Touch catheter (Biosense-Webster) or a Tacticath catheter (Abbott). Alternatively separate catheters may be used, one such as the AFocus catheter for diagnosis or location of the source of fibrillation, and the other for example a Navistar catheter for ablation.

[0062] The variety of directions of activation wave fronts at a particular location is illustrated for example in FIG. 2. In the region shown, there are two focal sources 200, 202 shown both of which act as sources of activation wave fronts 204. Due to the irregularity of, and lack of synchronisation between, the two sources, the activation waves from the two sources interact in a manner which varies considerably over time. Therefore if a point such as that 206 shown in FIG. 2 is considered, the direction in which wave fronts will cross that point will in general be different for each wave front. While some, such as the first, second and fourth shown, may be generally in the direction away from one of the sources, others, such as the third shown, will be in different directions as a result of the interference between waves from the two sources.

[0063] The system is therefore arranged to analyse the signals from the probe electrodes 102 so as to detect each of the different directions of propagation across each point, and then to analysis those as will be described in more detail below so as to locate, and enable treatment of, the rotors or other problematic regions.

[0064] The data acquisition, data processing, and image display will now be described in more detail. The processor 108 is arranged to perform each of these steps. For any particular position of the catheter 100, a stream of raw signal data is acquired from each of the numerous electrodes 102 of the catheter. The position of each electrode 102 is known through one of a variety of methods well known to those skilled in the art, such as those marketed as CARTO™ or NavX™. The following steps are then carried out by the system under the control of the processor 108. They are described here in sequence but they can occur almost simultaneously so that the operator sees the activation pattern at any position of the catheter 100 after only a minimal delay, for example less than a second.

[0065] Firstly the electrical data is acquired. For this step the catheter 100 and computer 104 are arranged to acquire unipolar or bipolar electrogram data. A standard definition of unipolar electrogram data for a particular site is the potential difference between an intracardiac electrode at that site and a reference potential, for example at Wilson's central terminal, or any other combination of skin surface electrodes. Alternatively, a unipolar electrogram can be defined as the potential difference recorded between an intracardiac electrode and an electrode placed within the body at a site outside the heart, for example in the inferior vena cava, a large vein adjacent to the heart in which an electrode can very conveniently be located. Therefore for unipolar electrogram data a further electrode, not shown, is also provided and connected to the computer to provide the reference signal in known manner.

[0066] Alternatively bipolar electrogram data can be used, being defined as the potential difference between two of the intracardiac electrodes 102. In this case no further reference electrode is needed.

[0067] Whether unipolar or bipolar electrogram signals are used, the electrical signal (voltage) from each electrode (or electrode pair) is sampled at a regular sample frequency and the sampled values stored in memory for analysis.

[0068] Then, the electrogram data obtained is filtered to remove noise and baseline artefact. A variety of filtering algorithms are well known to those skilled in the art. It is possible to apply one or more in sequence, using software programs coded operating on the microcomputer system 104 as in this embodiment. In other embodiments the processing is performed by hardware circuitry specifically designed or customised for filtering, known as digital signal processing hardware. For example a simple band pass filter may be used, which may be at 10-250 Hz. An example of filtered electrogram signals is shown in FIG. 3, in which the electrogram data for electrodes numbered 12, 13, 14, 15 and 16 are shown, together with the results of analysis of pairs of the electrograms which are shown on lines 12.5, 13.5, 14.5 and 15.5, as will be described in more detail below.

[0069] Referring to FIG. 4 for the following discussion, the electrode locations at which the electrograms are obtained will be referred to as ‘nodes’. The nodes are triangulated, i.e. grouped into sets of three nodes 102a, 102b, 102c which form a triangle. Each edge 104 of the triangle extends between a pair of the three nodes, e.g. 102a and 102b, with a corresponding pair of electrograms. Each of these three pairs of electrograms is then processed to determine a time delay between the activation times of the two electrodes in the pair. From the three time delays, the direction and speed of the wave front passing the group of three nodes can be determined. For example each direction may be calculated and defined as components in two orthogonal directions x, and y, of a vector of unit length in the direction of propagation.

[0070] In order to determine each of the time delays an autocorrelation algorithm is used, performed by the processor 108 to compare the activation times of the two electrodes in each electrogram pair. Autocorrelation is used because closely spaced electrograms usually have a morphology that is similar. Autocorrelation provides a way of determining relative times of activation between the pair.

[0071] Let represent an electrogram time series with n samples:

e=[E.sub.1,E.sub.2,E.sub.3, . . . ,E.sub.N]

[0072] A windowing function is applied to the time series in order to select information that corresponds to a particular time and avoids multiple activations falling within the same window.

w=[w.sub.−k,w.sub.−k+1, . . . ,w.sub.−1,w.sub.0,w.sub.1, . . . ,w.sub.k−1,w.sub.k]

[0073] This is shown in FIG. 4, with k=25. The total width of the windowing function is 51 (2k+1) samples, corresponding to approximately 100 ms at a sample rate of 500 Hz. As can be seen from FIG. 5, the windowing function is centred on one sample, indicated as sample 0 in FIG. 5, and has its highest value for sample 0. It falls away to zero at the samples at either end of the window, i.e. −k and +k.

[0074] Next, the windowed electrogram e.sub.w(at time t, corresponding to one of the sample times, is obtained by multiplying the electrogram values, centred on the sample at time t, with the values of the windowing function:

ew(t)=[e.sub.t−k.Math.w.sub.−k,e.sub.t−k+1.Math.w.sub.−k+1, . . . ,e.sub.t−1.Math.w.sub.−1,e.sub.t.Math.w.sub.0,e.sub.t+1.Math.w.sub.1, . . . ,e.sub.t+k−1.Math.w.sub.k−1,e.sub.t+k.Math.w.sub.k]

[0075] Then windowed electrograms from two neighbouring nodes are selected—e.sub.1w and e.sub.2w. Next, cross correlation R(t,Δ) between the two windowed electrograms is determined for all sample times t and all possible time delays (i.e. differences in time of activation between the two electrodes in the pair) Δ within a range:

[00001]$R\left(t,\Delta \right)=e_{1w}\left(t-\frac{\Delta }{2}\right).Math.{\left(e_{2w}\left(t+\frac{\Delta }{2}\right)\right)}^{\tau }$

[0076] For example, for a catheter where electrograms are recorded 10 mm apart, allowing for slow conduction of 0.5 m/s, this corresponds to a maximum transit time of 20 ms. Therefore, values of A may be limited to ±20 ms.

[0077] Next, to find values for t and Δ local maxima in R(t, Δ) are determined. Only maxima above a sensitivity threshold are considered. At each local maximum, R(t.sub.max, Δ.sub.max), the relative timing of two electrodes e.sub.1 and e.sub.2 is given by Δ.sub.max. Approximate timing T.sub.e1 and T.sub.e2 of activation at each of the two electrodes in the pair can be determined as follows:

[00002]$T_{e1}=t_{\max }-\frac{{\Delta }_{\max }}{2}{T}_{e2}=t_{\max }+\frac{{\Delta }_{\max }}{2}$

[0078] Information from multiple electrograms can, optionally, be compared by performing the above analysis on, for example, all possible pairs from groups of three electrograms, with their positions, after they have been triangulated.

[0079] The node positions 102a, 102b, 102c on the catheter are then transformed onto a 2 dimensional surface and a grid defined, with a spacing that is substantially smaller than the inter-electrode distance. At each point g on the grid, the nearest edge 104 is determined and the relative distances from the two nodes at the ends of that edge are calculated. The activation timing T.sub.g at each grid point g is then calculated as:

[00003]$T_{ℊ}=t_{\max }-\frac{{\Delta }_{\max }}{2}+\alpha .{\Delta }_{\max }$ [0080] where α has a value of 0 if the grid location g is next to node 1 and has a value of 1 if the location is next to node 2. Each grid value is assigned to 1 at times where there exists t−T.sub.g<100 ms.

[0081] This produces a set of activation times, one of each grid point g. However because of the method of calculation, the times will not be representative of a smoothly propagating wave front. Standard image smoothing algorithms are then used to create a visual display of smooth wavefront propagation. This may be performed, for example, by applying a box filter and then thresholding. Using the smooth wavefront data, a Sobel edge detector, or other suitable edge detector, is used to identify the wavefront direction at each grid point.

[0082] From the smoothed timing information describing wavefront propagation times at each grid point, direction data at each grid point is calculated.

[0083] Let the direction of a wavefront passing grid point, g, at time, t, be d(g, t)=x(g, t). i+y(g, t). j, where x(g, t) and y(g,t) are the magnitudes in two orthogonal directions (i and j) and √{square root over (x(g, t).sup.2+y(g, t).sup.2)}=1. That direction is determined using triangulation from any group of three grid points. For example it may be calculated for each grid point using a triangular group of that grid point and two adjacent grid points g. Now the wavefront propagation across an area of interest can be described by taking the most recent activation at each point within a specified time range. The direction of the latest wave front to pass a grid point, which is a vector of unit length

d.sub.latest(g,T)=find d(g,t) with t<T

[0084] Now, for any defined surface (usually the surface covered by the electrodes that have been analysed) can be subjected to analysis.

[0085] Wavefront direction, at a given sample time, for a particular wavefront travelling over the analysed region G(T) is calculated as the integral of latest directions for all grid points over a fixed time period, for example 200 ms, which will typically be approximately long enough so that each transition period is covered completely by one integration period, though this will of course not always be the case:

[00004]$G\left(T\right)=\underset{s}{\oint }{d}_{\mathrm{latest}}\left(ℊ,T\right)\mathrm{ds}$

[0086] where, the wavefront direction is given by a tan(G.sub.y(T)/G.sub.x(T)) and the wavefront ‘coherence’ is given by ∥G(T)∥.

It will be appreciated that, for a given wave front, G(T) is a vector sum of directional vectors d.sub.latest over the transition period during which the wave front is passing through the analysed region. Therefore, assuming the wave front moves in one direction the length of the vector G(T) will increase over the transition period. If there is variation in the direction of travel of the wavefront, then the direction of the vector G(T) will also vary over that period.

[0087] This can be displayed as a series of dots on an x-y plot having an origin, with the dots each being located at a point which is displaced from the origin by a distance and direction corresponding to the length and direction of the integral G(T) at the time the dot is added to the display. A new dot may be displayed at regular update intervals during the transition period, and each dot displayed may be displayed for a display period, which is much longer than those regular update intervals, so that the dots are superimposed on the image during the transition period. At the end of its display period each dot may be removed from the display, or faded out. Alternatively the dots may be each be displayed continuously so that the number of dots displayed increases until the end of the measurement. This addition of a series of dots generates a line of dots for each wavefront that starts at the origin and is extended after each update interval in the direction of travel of the wavefront during that update interval. Therefore if the direction of travel is constant, the line will be straight, whereas if the direction of travel varies over the transition period, the line will be curved. This display therefore gives information about the direction and coherence of wavefronts. When a wavefront has completed its transition of the analysed region no further dots will be added to the line representing that wavefront. When a new wavefront is detected entering the analysed region, e.g. when the value of the integral G(T) changes from zero on one update interval to a non-zero value in the next update interval, a new dot is displayed at the origin and a new line of dots started for that wavefront. FIG. 6 shows the appearance of the plot as a wavefront passes the recording electrodes in a direction from bottom-right towards top-left. Panels are shown at 15 ms intervals, which is equal to the sample interval. The dots in the example shown are roughly circular, but it will be appreciated that they can be of a variety of shapes such as square or triangular.

[0088] In a modification to the process described above, the integral G(T) is calculated over much shorter time periods, for example once for each sample period. Then a new dot is displayed for each sample period, the location of the dot being determined not relative to the origin, but relative to the location of the previous dot in the line, with the new dot being offset from the previous dot in the direction of the wavefront direction just in the latest sample period.

[0089] Referring to FIG. 7, the information on wavefront direction can also be amalgamated into a polar histogram plot, weighted according to the number and coherence of wavefronts in each direction. In order to do this the wavefront direction is calculated at regular intervals, for example using the integral G(T) as described above either once for each sample interval, or once every two or three sample intervals. This generates a set of sample wavefront directions which can be displayed on a polar histogram as shown in FIG. 7.

[0090] The vector quantities d.sub.latest at each of the grid points any one time define a vector field over the detected region and can therefore be analysed to determine the vector operators curl and divergence at points within that region. The curl can be defined as:

C(g,T)=curl(d.sub.latest(g,T))

and the divergence can be defined as

D(g,T)=div(d.sub.latest(g,T))

[0091] These are calculated for each point g on the grid from the values of d.sub.latest. The divergence and curl are calculated at each grid point g using the vector that represents the last wavefront direction at neighbouring grid points. (Thus the speed of the wavefront is not used.) The divergence and curl are calculated using standard algorithms. In a neighbourhood, let x and y be the distance along two orthogonal vectors (i and j) on the heart surface. At each grid point g, the wavefront direction, d, can de represented as the sum of two vectors:

d(x,y)=ui+vj

[0092] Where d(x,y) has been normalised to have a magnitude of 1.

[0093] Divergence and Curl are then calculated for each grid point g as:

[00005]$\mathrm{divergence}-\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$$\mathrm{curi}=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$

[0094] The wavefront divergence and curvature may be used to highlight locations where there is rotor (or rotational) activity and also where there are focal sources or wavefront collision. Specifically maxima, or high values, of curl are associated with rotors, and maxima or values of high divergence are associated with focal sources. Therefore these maxima or high values can be located as described below and used to indicate the location of tissue that can be ablated.

[0095] Referring to FIGS. 8A-8D, the operation of the system can be described as follows. Firstly at step A, shown in FIG. 8A, the electrogram data for each of the electrodes 102 on the probe 101 is acquired for a first position of the probe 101 and a corresponding region of the heart. Then at step B, shown in FIG. 8B, the electrogram data for that position over a sample time period is analysed as described above to determine the direction of propagation, and optionally also the coherence, of each wave front passing through that region in the total sampling time is determined. In the example shown, the probe 101 has a spiral array of electrodes 102, with the electrode positions or nodes shown in FIG. 8B. At step C, shown in FIG. 8C, data relating to the directions of wavefront propagation is displayed as an image on the display screen 110, for example as vectors on a polar plot, or a moving dot pattern as described above with reference to FIG. 6, for the region in which the probe 101 is located. The probe 101 is then moved to a second position and then subsequent positions, and in each position further data acquired, analysed and displayed. The data for each position is stored, together with probe location data indicative of the location of the probe and hence the position on the heart of the region currently being examined. At step D, shown in FIG. 8D, from this location data, the individual wave front direction images are mapped onto an image of the heart. The individual wavefront direction images may be shown as vector plots or dot images, or indeed as histogram plots or simple directional indicators as shown in FIGS. 8C and 8D.

[0096] On the image of the heart, the location of sources of fibrillation may then be determined by the user from the directional information displayed. Alternatively the processor 108 may be arranged to determine the location of sources and to control the display 110 to indicate the location of the sources on the image. For example the processor may be arranged to calculate the curl of the wavefront direction vector at positions on the heart as described above, locate a maximum of the curl, and identify the position of that maximum as the position of a rotor. The processor may then be arranged to control the display 110 to highlight, for example using arrows or colour or an outline, the position on the heat of the rotor. Alternatively the processor may be arranged to calculate the divergence of the wavefront direction vector at positions on the heart as described above, locate a maximum of the divergence, and identify the position of that maximum as the position of a focal source. The processor may then be arranged to control the display 110 to highlight, for example using arrows or colour or an outline, the position on the heat of the focal source.

[0097] Once the location of the rotor has been determined, the ablation tip 114 of the catheter, or a separate ablation catheter, is used to ablate heart tissue at the location of the rotor.

[0098] Referring to FIG. 9, the system of FIG. 1 can be arranged to differentiate between further different wave front patterns using the divergence and curvature (curl) of the wave front directions. The first set of four images in FIG. 9 shows the detection of a relatively uniform wave front passing across the region of the heart that is being contacted by the probe 101. The second set of four images shows the detection of a wave front propagating out from a focal point in the imaged region of the heart. The third set of four images shows a wave front that is rotating about a point generally at the centre of the imaged region of the heart. The fourth set of four images shows the detection of a wave front that is turning about a point outside (in this case below as seen in FIG. 9) the imaged region of the heart. Finally the last set of four images in FIG. 9 shows the detection of a pair of wave fronts colliding at the centre of the imaged region of the heart.

[0099] There are various ways of identifying the time of activation at each node, i.e. the position of each electrode, and identifying, monitoring, and categorizing the activation wavefronts as they pass through the monitored region. Referring to FIG. 10, in one approach the processor is arranged to identify at each of a series of update intervals those electrodes at which an activation has been detected. A map of the monitored region of the heart can then be defined which identifies an area 300 in which electrodes have experienced an activation and another area 302 in which the electrodes have not experienced an activation. The top row of FIG. 10 shows a sequence of maps for a corresponding sequence of update intervals in which the activated area 300, which is shown in black, moves or expands across the monitored region from bottom left to top right. In each map, the boundary 304 between the activated area 300 and the non-activated area 302 defines the wavefront.

[0100] The bottom row of FIG. 10 shows the current wavefront after each update interval, together with the locations of the wavefront after each previous update interval. The wavefronts may be displayed as an image which updates after each update interval and adds the latest wavefront position to those already displayed. This would result in an image that changes from the left hand image to the right hand image of the bottom row of FIG. 10 as the wavefront crosses the monitored region, moving from position 304a to 304e through the intermediate positions of 304b, 304c, 304d in between. For emphasis the latest wavefront position after each update interval may be indicted more brightly than the others as shown in FIG. 10. Each wavefront may be further faded in each subsequent image update in which it is displayed.

[0101] As described above the system may be arranged to characterize the wavefronts. In one method of doing this, the concept of child and parent wavefronts is used. A wavefront is a child wavefront of a wavefront that precedes it, for example 304e is a child of 304d, and a wavefront is a parent wavefront of a wavefront that follows it, so 304a is a parent of 304b, and 304b is a parent of 304c. The point at which a wavefront has an end, either within the image or where the wavefront meets the edge of the image (i.e. the edge of the monitored region) is defined as an endpoint 306.

[0102] These features allow various wavefront shapes, or types of activation, to be identified, as shown in FIG. 11.

[0103] Plane wavefront (a) satisfies all of: [0104] A wavefront with both endpoints at the edge of the field of view (i.e. on the outer circle). [0105] Either 1 parent or no parents. [0106] Either 1 child or no children.

[0107] Focal wavefront (b) satisfies all of: [0108] Either 1 parent wavefront which is also focal, or no parents [0109] The endpoint of the wavefront meets the wavefront (i.e. the wavefront is a loop).

[0110] Rotating wavefront (c) satisfies all of: [0111] One endpoint is within the field of view (i.e. one endpoint is ‘floating’ within the outer circle).

[0112] Collision wavefront (d) satisfies all of: [0113] 2 or more parent wavefronts.

[0114] These categorizations may be used displayed to a user in various ways, for example using different colours on a display for different categorizations of wavefront, or using labelling, depending on the purpose of the investigation being carried out.

[0115] Referring to FIG. 12, the raw (un-filtered) electrogram data from the individual electrodes, plotted as voltage V as a function of time t, includes both the underlying unipolar voltages, which include sharp drops caused by the activations, and noise which is superimposed on the underlying voltages, and shown as a separate signal varying about each of the underlying voltage signals. The noise is typically of a relatively high frequency, but needs to be filtered out without removing the sharp drops in the underlying voltages. In order to do this, it can be assumed that the noise in the signals will be, in general, common to all, or several, of the individual signals. This means that features having (approximately) the same shape occur at (approximately) the same time in some or all of the individual signals. In contrast, the sharp variations in voltage which are used to identify the activations occur at different times in the different signals and are therefore not common, at least in timing, to all of the signals. A filtering method is therefore used which takes into account not just the individual unipolar electrograms, but measurements of the difference between voltages at a pair of electrodes. This enables the system to distinguish between noise features and activation features, based on the similarities and differences between the individual unipolar electrograms. This enables the system to filter the individual electrogram signals so as to filter out the noise features and retain the activation features. While a difference between two voltages is a simple and effective measure, it will be appreciated that other linear combinations of two or more signals, for example including different weightings applied to the different signals, could also be used.

[0116] Firstly to understand the advantages of including the voltage differences in the analysis, it is helpful to consider the situation in which the time varying aspect of the voltages is ignored. Consider the situation where there are two electrodes (A and B). The voltage at each electrode is measured, and also the voltage difference between the two electrodes is measured directly. At a particular time the following data is recorded: V.sub.a=0.8, V.sub.b=1.2, V.sub.b-a=0.2. It is known that the direct measurement V.sub.b-a is subject to much less noise than V.sub.a and V.sub.b. The best estimate of the underlying state is that V.sub.aEst=0.9, V.sub.bEst=1.1. This estimate preserves the very accurate measurement of the difference between the voltages (0.2) and has assumed an equal magnitude of error in V.sub.a and V.sub.b.

[0117] Then considering a time-varying system with just two electrodes, measurements can be made of the potential difference between the electrodes and the skin—two unipolar voltage signals, u1 and u2. Measurement can also be made of the bipolar voltage, b.sub.1=u.sub.1-u.sub.2. The aim is to estimate unipolar signals utilising the information from all the measurements that have been made: in clinical practice, the bipolar measurements (i.e. differences between unipolar electrograms) have much less noise because of common mode noise rejection.

[0118] The aim is to estimate the underlying state, i.e. the actual voltages at each of the electrodes:

X=(u.sub.1,u.sub.2)  [1]

[0119] The noisy measurements Y are used, which are made up of the noisy measurements of u.sub.1, u.sub.2 and b.sub.1:

{tilde over (Y)}=(ũ.sub.1,ũ.sub.2,{tilde over (b)}.sub.1)=XH+noise  [2]

[0120] Where H is an observer matrix:

[00006]$\begin{array}{cc}H=\left(\begin{array}{ccc}1& 0& 1\\ 0& 1& -1\end{array}\right)& \left[3\right]\end{array}$

[0121] It will be noted that the observer matrix defines the linear combination, here a simple difference between the two electrograms, which is used in addition to the individual unipolar electrograms. For the more general case with more than two electrodes, the observer matrix defines for which pairs of electrodes the voltage difference is measured and used.

[0122] In a corresponding manner, an observer matrix can be used in a system with multiple (more than two) electrodes and time varying signals to estimate the unipolar voltages. For example a method based on a Savitsky Golay filter may be used. This estimates each of the unipolar signals to be modelled by a polynomial having a number of terms with respective weights, and estimates the weights using a least squares fit to the measured data, which includes the unipolar electrograms and measurements of the difference between pairs of the electrodes.

[0123] Consider a timeframe z that consists of M samples at interval dt around t.sub.0.

[00007]$\begin{array}{cc}z=\frac{t-t_0}{\mathrm{dt}}=\left(\frac{-\left(M-1\right)}{2},\frac{-\left(M-3\right)}{2},.Math.0,.Math.,\frac{\left(M-3\right)}{2},\frac{\left(M-1\right)}{2}\right)& \left[4\right]\end{array}$

[0124] For each individual electrogram recording assume it is well modelled by a local polynomial with order n:

y.sub.k(z)=β.sub.0+β.sub.1z+β.sub.2z.sup.2+β.sub.3z.sup.3+β.sub.4z.sup.4+ . . . +β.sub.nz.sup.n  [5]

[0125] This can be written as:

y.sub.k=Vβ.sub.k  [6]

[0126] Where V is a Vandermonde matrix:

[00008]$\begin{array}{cc}V=t+\mathrm{dt}\times \left(\begin{array}{ccccc}1& \frac{-\left(M-1\right)}{2}& {\left(\frac{\left(M-1\right)}{2}\right)}^2& .Math.& {\left(\frac{-\left(M-1\right)}{2}\right)}^n\\ 1& \frac{-\left(M-3\right)}{2}& {\left(\frac{\left(M-3\right)}{2}\right)}^2& .Math.& {\left(\frac{-\left(M-3\right)}{2}\right)}^n\\ .Math.& .Math.& .Math.& & .Math.\\ 1& 0& 0& .Math.& 0\\ .Math.& .Math.& .Math.& & .Math.\\ 1& \frac{\left(M-3\right)}{2}& {\left(\frac{\left(M-3\right)}{2}\right)}^2& .Math.& {\left(\frac{\left(M-3\right)}{2}\right)}^n\\ 1& \frac{\left(M-1\right)}{2}& {\left(\frac{\left(M-1\right)}{2}\right)}^2& .Math.& {\left(\frac{\left(M-1\right)}{2}\right)}^n\end{array}\right)& \left[7\right]\end{array}$$\mathrm{and}$$\begin{array}{cc}{y}_k={\left(\frac{y_{-\left(M-1\right)}}{2},\frac{y_{-\left(M-3\right)}}{2},.Math.,\frac{y_{\left(M-3\right)}}{2},\frac{y_{\left(M-1\right)}}{2}\right)}^T,{\beta }_k={\left({\beta }_0,{\beta }_1,{\beta }_2,{\beta }_3,{\beta }_4.Math.{\beta }_n\right)}^T& \left[8\right]\end{array}$

[0127] Now, if all of the measurements are considered they can be represented as:

Y=VB,Y=(y.sub.1,y.sub.2,y.sub.3, . . . ,y.sub.K),B=(β.sub.1,β.sub.2β.sub.3, . . . ,β.sub.K)  [9]

[0128] The measurements y.sub.k (k=1,2, . . . ,K) are related to the underlying states x.sub.j (j=1,2 . . . J) via the observer matrix H:

Y=XH,X=(x.sub.1,x.sub.2, . . . ,x.sub.j),H=J×K matrix  [10]

[0129] Where j is the number of electrodes and k is the number of measured parameters which may be individual unipolar voltages, differences between pairs of voltages, or other linear combinations of unipolar voltages. It follows that

x.sub.j(z)=α.sub.0+α.sub.1zαα.sub.2z.sup.2+α.sub.3z.sup.3+α.sub.4z.sup.4+ . . . +α.sub.nz.sup.n  [11]

And

X=VA,X(x.sub.1,x.sub.2, . . . ,x.sub.j),A=(α.sub.1,α.sub.2,α.sub.3, . . . ,α.sub.j)  [12]

With

B=AH  [13]

[0130] Then the aim is to find the J×(N+1) polynomial coefficients of A such that R, the weighted least squares residual in the estimated signals ŷ is minimised:

[00009]$\begin{array}{cc}R=\underset{k=1}{\overset{K}{.Math.}}{R}_k,R_k=\underset{z=\frac{-\left(M-1\right)}{2}}{\overset{\frac{\left(M-1\right)}{2}}{.Math.}}{w_{k,z}\left({\hat{y}}_{k,z}-{\tilde{y}}_{k,z}\right)}^2& \left[14\right]\end{array}$

[0131] This can be re-written in matrix form with W.sub.k being a matrix of weights w.sub.k,z on its diagonal:

R.sub.k=(ŷ.sub.k−Vβ.sub.k).sup.TW.sub.k(ŷ.sub.k−Vβ.sub.k)  [15]

[0132] To minimize R, we need

[00010]$\frac{\mathrm{dR}}{d{\alpha }_j}$

to be zero for all j. First differentiate with respect to β.sub.k:

[00011]$\begin{array}{cc}\frac{\mathrm{dR}}{d{\beta }_k}=\underset{k=1}{\overset{K}{.Math.}}-2V^T{W}_k{y}_k+2V^T{W}_k{V}^T{\beta }_k& \left[16\right]\end{array}$

[0133] Then using the chain rule

[00012]$\begin{array}{cc}\frac{\mathrm{dR}}{d{\alpha }_j}=\frac{\mathrm{dR}}{d{\beta }_k}\frac{d{\beta }_k}{d{\alpha }_j}=\underset{k=1}{\overset{K}{.Math.}}\left(-2V^T{W}_k{y}_k+2V^T{W}_k{V}^T{\beta }_k\right)h_{\mathrm{jk}}& \left[17\right]\end{array}$

[0134] Where h.sub.jk is the element of H form row j and column k.

[0135] Substituting β.sub.k=AH.sub.k, H.sub.k=(h.sub.1k, h.sub.2k, h.sub.3k, . . . , h.sub.jk).sup.T and setting

[00013]$\frac{\mathrm{dR}}{d{\alpha }_j}=0,$

we obtain

[00014]$\begin{array}{cc}\underset{k=1}{\overset{K}{.Math.}}{h}_{\mathrm{jk}}{V}^T{W}_k{y}_k=\underset{k=1}{\overset{K}{.Math.}}\left(V^T{W}_k{V}^T{\mathrm{AH}}_k\right)h_{\mathrm{jk}}\mathrm{for}\mathrm{all}j& \left[18\right]\end{array}$

[0136] In order to solve equation [18], rearrange A and Y as vectors:

[00015]$\begin{array}{cc}\stackrel{.Math.}{Y}=\left(\begin{array}{c}{y}_1\\ y_2\\ y_3\\ .Math.\\ y_K\end{array}\right)\left(K\times M\mathrm{column}\mathrm{vector}\right),\stackrel{.Math.}{A}=\left(\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ {\alpha }_3\\ .Math.\\ {\alpha }_J\end{array}\right)\left(J\times \left[n+1\right]\mathrm{column}\mathrm{vector}\right)& \left[19\right]\end{array}$

[0137] [18] can then be re-written as

[00016]$\begin{array}{cc}\left(\begin{array}{cccc}{h}_{11}{V}^T{W}_1& h_{12}{V}^T{W}_2& .Math.& h_{1K}{V}^T{W}_K\\ h_{21}{V}^T{W}_1& h_{22}{V}^T{W}_2& .Math.& h_{2K}{V}^T{W}_K\\ .Math.& .Math.& \ddots & .Math.\\ h_{J1}{V}^T{W}_1& h_{J2}{V}^T{W}_2& .Math.& h_{\mathrm{JK}}{V}^T{W}_K\end{array}\right)\stackrel{.Math.}{Y}=\left(\begin{array}{c}\underset{k=1}{\overset{K}{.Math.}}{h}_{1k}\left(h_{1k}{V}^T{W}_{k}V,h_{2k}{V}^T{W}_{k}V,.Math.,h_{\mathrm{Jk}}{V}^T{W}_{k}V\right)\\ \underset{k=1}{\overset{K}{.Math.}}{h}_{2k}\left(h_{1k}{V}^T{W}_{k}V,h_{2k}{V}^T{W}_{k}V,.Math.,h_{\mathrm{Jk}}{V}^T{W}_{k}V\right)\\ .Math.\\ \underset{k=1}{\overset{K}{.Math.}}{h}_{\mathrm{Jk}}\left(h_{1k}{V}^T{W}_{k}V,h_{2k}{V}^T{W}_{k}V,.Math.,h_{\mathrm{Jk}}{V}^T{W}_{k}V\right)\end{array}\right)\stackrel{.Math.}{A}& \left[20\right]\end{array}$

[0138] Finally an expression can be obtained for in terms of

=Φ  [21]

[0139] Where Φ is a matrix with [n+1]×J rows and [n+1}×K columns

[00017]$\begin{array}{cc}\Phi ={\left(\begin{array}{c}\underset{k=1}{\overset{K}{.Math.}}{h}_{1k}\left(h_{1k}{V}^T{W}_{k}V,h_{2k}{V}^T{W}_{k}V,.Math.,h_{\mathrm{Jk}}{V}^T{W}_{k}V\right)\\ \underset{k=1}{\overset{K}{.Math.}}{h}_{2k}\left(h_{1k}{V}^T{W}_{k}V,h_{2k}{V}^T{W}_{k}V,.Math.,h_{\mathrm{Jk}}{V}^T{W}_{k}V\right)\\ .Math.\\ \underset{k=1}{\overset{K}{.Math.}}{h}_{\mathrm{Jk}}\left(h_{1k}{V}^T{W}_{k}V,h_{2k}{V}^T{W}_{k}V,.Math.,h_{\mathrm{Jk}}{V}^T{W}_{k}V\right)\end{array}\right)}^{-1}\left(\begin{array}{cccc}{h}_{11}{V}^T{W}_1& h_{12}{V}^T{W}_2& .Math.& h_{1K}{V}^T{W}_K\\ h_{21}{V}^T{W}_1& h_{22}{V}^T{W}_2& .Math.& h_{2K}{V}^T{W}_K\\ .Math.& .Math.& \ddots & .Math.\\ h_{J1}{V}^T{W}_1& h_{J2}{V}^T{W}_2& .Math.& h_{\mathrm{JK}}{V}^T{W}_K\end{array}\right)& \left[22\right]\end{array}$

[0140] The matrix Φ is dependent only upon the observation matrix H and the weights W.sub.k, that have been chosen for each signal measurement. For continuous signals the estimate of the underlying states, {circumflex over (X)} can be calculated efficiently as a convolution operation with Y.

[0141] The weights W.sub.k can also be considered as a series of windows for each signal measurement. Selecting a long timeframe (M—see equation [4]), a low polynomial order (n—see equation [5]) and a uniform window W.sub.k=[1 1 1 . . . 1] will result in a low-pass filter being applied to signal k. Conversely a window with weights that are high only in the central regions will preserve higher frequencies. For example, using W.sub.k=[0 0 1 1 1 0 0] results in a window length of 3 being applied to signal k and allows window lengths of up to 7 to be used for other signals.

[0142] As mentioned above, the observer matrix H can be selected to include useful linear combinations of the unipolar voltages. In one example, as well as the individual unipolar voltages, the difference between the voltage at each electrode and the adjacent electrode is included, i.e. the unipoles and ‘sequential bipoles’. For four electrodes 1 to 4, this would be V(1) V(2) V(3) V(4) V(2-1) V(3-2) V(4-3). In some cases it is advantageous to measure more of the bipole ‘options’ because the noise in any bipole signal may be dependent on how far apart the two electrodes are. If the electrode positions are known, then the weighting of the errors used in the minimisation can be arranged to relate to the expected noise level.

[0143] The weights W.sub.k, that are chosen for each signal measurement can be used to select the relative weight of the different signal measurements. In this method they can also provide filtering in that they can be used to select the relative weight of different frequency components for each signal measurement as described above. For example in the bipolar electrograms there is a reduction of noise but also a loss of the high frequency signal components. However by incorporating bipolar electrograms in the model it is possible to recover the high frequency information. This is done by using weights (e.g. W.sub.k in equation 22) that are different for the bipolar signals and preserve the high frequency information. The measurement model enables the combination of this high frequency information with the lower frequency information that is present in the filtered unipolar electrograms.

[0144] FIG. 13 shows the results of a filter as described above, in which the gradient (time differential) of one of the unipolar electrograms is plotted, as a function of time, for the raw unfiltered and the filtered electrograms. As can be seen the filtered data shows much less noise and the steep negative gradients which are indicative of activation events are much more clearly identified.

[0145] Once the electrogram data has been filtered, the time at which an activation event occurs at each electrode, i.e. the time at which a wavefront passes that electrode, needs to be identified. FIG. 14 shows examples of electrograms of three different activation events, shifted so that they occur at about the same time and are centred on V=0. It is noted that the shape of the electrogram plot around the time of an activation event is similar to a plot of 1/t as a function of t. Therefore in order to locate the activation events, a convolution with 1/t is applied to the filtered electrograms of FIG. 12 (for example a Hilbert filter) which results, after filtering out low frequency noise using a digital linear phase filter, in the plots shown in FIG. 15 with peaks clearly identifying the location in time of the activation events.