Journal of Healthcare Engineering
Volume 2018 (2018), Article ID 7901502, 14 pages
https://doi.org/10.1155/2018/7901502

Research Article

Enhancement of the Comb Filtering Selectivity Using Iterative Moving Average for Periodic Waveform and Harmonic Elimination

José L. Ferreira,¹ Yan Wu,¹ and Ronald M. Aarts1,2

¹Department of Electrical Engineering, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, Netherlands
²Philips Research Laboratories, Eindhoven, Prof. Holstlaan 4, 5656 AE Eindhoven, Netherlands

Correspondence should be addressed to José L. Ferreira; j.l.ferreira@tue.nl

Received 10 August 2017; Accepted 28 November 2017; Published 1 February 2018

Academic Editor: Saugat Bhattacharyya

Copyright © 2018 José L. Ferreira et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A recurring problem regarding the use of conventional comb filter approaches for elimination of periodic waveforms is the degree of selectivity achieved by the filtering process. Some applications, such as the gradient artefact correction in EEG recordings during coregistered EEG-fMRI, require a highly selective comb filtering that provides effective attenuation in the stopbands and gain close to unity in the pass-bands. In this paper, we present a novel comb filtering implementation whereby the iterative filtering application of FIR moving average-based approaches is exploited in order to enhance the comb filtering selectivity. Our results indicate that the proposed approach can be used to effectively approximate the FIR moving average filter characteristics to those of an ideal filter. A cascaded implementation using the proposed approach shows to further increase the attenuation in the filter stopbands. Moreover, broadening of the bandwidth of the comb filtering stopbands around −3 dB according to the fundamental frequency of the stopband can be achieved by the novel method, which constitutes an important characteristic to account for broadening of the harmonic gradient artefact spectral lines. In parallel, the proposed filtering implementation can also be used to design a novel notch filtering approach with enhanced selectivity as well.

1. Introduction

In biomedical signal processing and signal processing in general, comb filtering approaches represent an important class of filters that play a relevant role in different fields, such as extraction or elimination of periodic signal and harmonic components, speech and audio signal processing, decimation processes, prediction and estimation of geophysical signals, and power line rejection [1–5]. In its simplest form, a comb filter can be viewed as a combination of notch filters in which the null frequencies occur periodically across the filter bandwidth. Another very popular comb filtering approach is the conventional FIR moving average filter indicated in

whose representation in z-domain and discrete time realisation is shown, respectively, in (2) and Figure 1.

with , where is the sampling frequency and is the fundamental of the periodic null frequencies.

Figure 1: Discrete time realisation of the moving average filter described in (2).

The comb filter realisation indicated in (2) is widely employed because of its computational efficiency. As limitations, however, (2) provides a magnitude response with low attenuation in the filter stopbands as well as nonuniform gain and high attenuation in the pass-bands, as depicted in Figure 2(a). In addition, despite the piecewise linearity of the phase characteristic (Figure 2(b)), it can provoke increased phase delay for higher values of [1, 4].

(a)

(b)

Figure 2: Frequency response of for some values of : (a) magnitude response; (b) phase response.

Such characteristics are undesirable in some applications and are far from those of an ideal comb filter: zero gain at notch frequencies, uniform and unity gain in the pass-bands, and no effects on the signal phase. To make the comb filter realisation of (2) more selective or closer to the ideal behaviour, some strategies have been suggested in the literature. For instance, it can be achieved by the introduction of poles in the transfer function of (2), as indicated in [4, 6, 7]

where the value of the parameter r is contained in the interval 0 ≤ r < 1. As mentioned by Proakis and Manolakis [4], the insertion of poles in (2) has the effect of introducing a resonance in the vicinity of the null, thus provoking reduction of the bandwidth of the stopbands. Thereby, the zeros placed at the unit circle in the z-plane will have in their vicinity the poles . The comb filter described in (3) has been successfully used in several applications, such as in harmonic compensators and rectifiers in power systems [6, 7]. However, one of the limitations of this approach is the decreased attenuation in the stopbands, as the value of increases towards 1. Furthermore, there is a trade-off between the values of and , which is contingent to the performance requirement of the filter: on one hand, the use of higher values of makes this method computationally expensive regarding memory usage. On the other hand, the exponential decrease of the power makes the filter to be implemented by using lower resolution computer unit. Thus, there exists a compromise between the value of and the computer unit resolution [6].

Another proposed strategy to improve the selectivity of the comb filtering provided by (2) is the time-domain averaging approach. Time-domain averaging consists of a kind of comb filtering approach based on a coherent detection process whereby estimation and elimination of the periodic activity are carried out by averaging repetitive sequences of a periodic signal, , observed in the input signal, .

In (4), represents the nonperiodic component of , which could be a noise signal or some stochastic process. Under the assumption that and are uncorrelated, summing up subsequent segments corresponding to the periodic signal results in coherent summation of [8]. Thus, the estimate of obtained by time-domain averaging can be calculated by the following discrete representation:

Or in z-domain,

The frequency response associated with the nonperiodic component is derived from the subtraction between the discrete representation of and (5).

where is the magnitude response of the comb filtering for elimination of the periodic component.

Time-domain averaging is a well-established comb filtering approach which has been widely used to estimate and extract periodic signals encountered in phenomena involving some rotating machinery [2, 8, 9]. Time-domain averaging-based approaches have been also proposed to estimate and eliminate the gradient artefact from the EEG signal, such as the average artefact subtraction method [10, 11]. The gradient artefact consists of a periodic waveform voltage interference which is induced in the electrical potential recorded in the human scalp (scalp potential) by the rapidly varying magnetic field gradients and radiofrequency pulses used in MRI sequences during simultaneous acquisition of EEG and fMRI data [12, 13]. One limitation of time-domain averaging, however, is its high dependency on accurate sampling of the periodic waveform . The occurrence of jitter errors may result in imprecise sampling of the averaging waveforms, which can impair the effectiveness of the method. Thus, the period of the repetitive waveform must be an exact multiple of the sampling interval. In parallel, the period of must be precisely known, which requires an external trigger or reference signal provided by an additional hardware [8, 14]. In case of the suppression of the gradient artefact from the EEG signal, subject movements or small drifts may also compromise the performance of the algorithm, since they change the morphology and shape of the artefact, in such a way that it is not possible to obtain an accurate estimate of . Subject movements or small drifts also provoke broadening of the harmonic artefact spectral lines [15], whose attenuation may not be effectively accounted for by the time-domain averaging comb filter. As a consequence, residual artefacts are left behind in the corrected EEG after subtraction of the estimated periodic waveform .

In this paper, we present a novel comb filtering implementation to improve the selectivity of the comb filtering provided by the FIR moving average filter of (2). As described in Section 2, implementation of such a comb filter has been based upon an iterative filtering decomposition process [16], whereby an estimation of the filtered signal can be obtained by the iterative application of a FIR moving average filter-based approach named double average filter. Comparison between the novel comb filtering implementation and those existing methodologies to enhance the selectivity of (2) described above shows that the novel method could be used in scenarios in which those approaches are not effective, such as during broadening of the harmonic gradient artefact spectral lines. In addition, the iterative application of time-domain averaging is revealed to enable the use of a smaller number of averages during application of such a method, as shown in Sections 3 and 4.

2. Methods

In recent research [16, 17], iterative filtering decomposition has been proposed as an alternative implementation for empirical mode decomposition [18]. According to this methodology, a series of low-pass filters (or moving average filters) is used to decompose a signal in intrawave frequency modes or intrinsic mode functions (IMFs). Here, we have further exploited the estimation of the first IMF, , by application of the filter () in the input signal :

where corresponds to a FIR moving average-based filter. The convergence of the iterative filtering decomposition is ensured by the coefficients (masks) of the filter having a value between 0 and 1, which has been demonstrated by Lin et al. [16].

2.1. Design of a Novel Comb Filtering Approach for Elimination of Periodic Waveforms

As , initially, we investigated the forward-backward application of the moving average filter indicated in (1) in the input signal, [19, 20]. This procedure allows obtaining a filtered signal with zero-phase distortion, which is a characteristic of an ideal comb filter. The forward-backward application of (1) in can be expressed as

Equation (10) is also referred to as double average filter [16, 17], where the coefficients of correspond to a triangular window of length . By applying the z-transformation in (10), it results in the following transfer function:

whose discrete time realisation is depicted in Figure 3.

Figure 3: Discrete time realisation of (11).

The frequency response of is derived from (11) by setting . Hence,

Figure 4 depicts the magnitude response of , calculated according to (12), for some values of . It also shows the presence of spaced zeros at the frequency 2π/M. For a hypothetical value , becomes an all-pass band filter.

Figure 4: Magnitude response of for some values of .

The phase response of possesses a zero-phase characteristic, as a result of the forward-backward application of the moving average filter of (1).

Therefore, (11) describes a kind of FIR moving average-based filter that provides no distortion effects in the phase of the signal in the whole filter pass-band. Replacing (11) by in (9) and taking into account a number of iterations, it can be rewritten as [19]

Equation (14) corresponds to the transfer function that relates the extracted periodic waveform, , and the input signal, , as likewise indicated in (4). Therefore, after elimination of , the output of the proposed comb filter has been related to as [20]

In order to improve the attenuation in the stopbands, we investigated the application of within the cascade implementation indicated in

where is the number of cascades. Since (15) and (16) have been derived from (11), which has zero-phase distortion characteristic, they do not cause any distortion effects on the filtered signal phase either.

2.2. Iterative Application of Time-Domain Averaging

As an alternative for the filter in (9), we have also investigated the use of the time-domain averaging filter described in (8). To this end, we have taken into account a number of iterations of (9), which has been rewritten as

Thereby, by eliminating the estimated periodic signal, the output has been related to the input as

2.3. Using the Proposed Method to Design a Novel Notch Filtering Approach

As remarked by Braun [2], the frequency response provided by (7) and (8) corresponds to the convolution between the frequency response depicted in Figure 2 and a train of unit pulses separated by the period . Making use of this idea, we also investigated the convolution of a single pulse, , and the magnitude response indicated in Figure 2(a) to design a novel notch filtering approach approximated to the ideal case.

Figure 5 shows this convolution, where the unit pulse has been located at the frequency , and was calculated for .

(a)

(b)

Figure 5: Convolution between the unit pulse δ (located at the frequency ) and (). It results in a pass-band filter with a central frequency at .

As can be noticed in Figure 5(b), (19) corresponds to a pass-band filter with unit amplitude and central frequency at . Replacing in (9) by and taking into account a certain number of iterations, it results in the notch filter of

Like in (16), we investigated the application of (20) in a cascade implementation to enlarge the attenuation at the notch frequency, as indicated in

where represents the number of cascades as well.

3. Results

3.1. Frequency Characteristics of the Novel Comb Filtering Approach

All frequency responses depicted below have been calculated using a number of samples samples, so that was set as , with a frequency interval at 2π × (1/Ns). In these figures, only the frequencies ranging from 0 to π are shown. The simulations were performed in MATLAB environment (The MathWorks Inc., Natick, USA).

In Figure 6, the magnitude response of (15) is depicted, taking into account and and some values of . It can be observed that increasing of is followed by substantial increasing of the filter gain, which attained 0 dB as well as became more uniform in the different pass-bands. On the other hand, increasing of is also followed by a reduction in the attenuation in the filter stopbands. On the extreme case, when J → ∞, (15) tends to become an all-pass band filter, as shown in Figure 7(a). By fixing the value of and varying , (15) shows to provide more uniform gain in the pass-bands for smaller values of , as observed in Figure 7(b).

(a)

(b)

Figure 6: Magnitude response of , taking into account (a) and (b) and some values of .

(a)

(b)

Figure 7: (a) Magnitude response of for and some values of J (J → ∞); (b) magnitude response of for and some values of M: (blue trace), (red trace), and (green trace).

To demonstrate the enhancement of attenuation in the stopbands using (16), its magnitude response has been calculated taking into consideration different values of . Figure 8 depicts the magnitude response of , for and , and taking into account , and some values of . As can be noticed, (16) provides higher attenuation in the stopbands by increasing the value of .

(a)

(b)

Figure 8: Magnitude response of (16) for (a) and (b) , taking into account and some values of .

3.2. Iterative Application of Time-Domain Averaging

By setting in (6), the magnitude and phase response of the periodic component obtained by time-domain averaging can be, respectively, derived as

and

Figure 9 depicts the magnitude and phase response of (22) and (23), taking into account , for (dark trace), (green trace), and (gray trace).

(a)

(b)

Figure 9: Frequency response of for , taking into account (dark trace), (green trace), and (gray trace): (a) magnitude response; (b) phase response.

The main lobes of the magnitude response lie spaced at the frequency 2π/M and have unit amplitude. By increasing , the amplitude of the side lobes is reduced as well as the main lobes are narrowed. Regarding the phase response, it is not linear, but has approximated piecewise linearity in between the null frequencies of the lobes [4, 8, 9]. In turn, Figure 10 shows the frequency response associated with the nonperiodic component , as indicated in (8).

Figure 10: Magnitude response of ((8) for and some values of .

Figure 11 depicts the frequency response provided by (17), taking into account and , for (dark trace), (green trace), and (gray trace).

(a)

(b)

Figure 11: Frequency response of (17), taking into account and , for (dark trace), (green trace), and (gray trace): (a) magnitude response; (b) phase response.

It can be noticed that increasing in (17) has a similar effect to increasing in (22), by reducing the amplitude of the side lobes and narrowing the main lobes. Regarding the phase response, it remains piecewise linear after increasing . In Figure 12, the frequency response corresponding to (18) is depicted, for , , and some values of . It can be observed that as increases, the filter gain is approximated to unity along the pass-bands as well as the stopbands become narrower. Therefore, an increase of in (18) attests a similar effect to that provided by increasing in (8) alike.

Figure 12: Magnitude response of ((18), taking into account and and some values of .

3.3. Notch Filter Characteristics Obtained by Iterative Filtering

The magnitude response corresponding to is depicted in Figure 13, taking into account , some values of , and samples for as well. The notch filter stopband was located at . It can be observed that the bandwidth of the notch filter stopband is narrowed by increasing the value of , as well as (20) tends to become an all-pass band filter when J → ∞. This notch filter does not cause any effects on the phase of the signal.

Figure 13: Magnitude response of ((20), taking into account some values of and cutoff frequency at .

In Figure 14, we show the use of (20) within the cascade implementation indicated in (21). The variation of the bandwidth of the notch filter stopband around −3 dB provided by (21) is indicated, taking into account , some values of and , , and  kHz. As can be seen, combination of proper values of and allows obtaining stopband bandwidths ranging from 0 (when J → ∞) up to 4500 Hz (for ). Therefore, (21) can produce a notch filter with a large range of values for the stopband bandwidth around −3 dB.

Figure 14: Bandwidth of the notch filter stopband around −3 dB ((21), taking into account and some values of and .  kHz.

4. Further Comparative Analyses and Discussion

Figure 15 shows the frequency response of (3), for some values of and .

(a)

(b)

Figure 15: Frequency response of for some values of : (a) magnitude response; (b) phase response.

It can be noticed that as the value of increases towards 1, the gain of the filter approximates to unity and becomes more uniform along the different pass-bands, as well as the phase response approaches to a zero-phase distortion. On the extreme scenario of r → 1, the filter becomes an all-pass band filter. On the other extreme, when , (3) equals (2) [6].

As indicated in Figures 6 and 10, by increasing the parameters ((15) and ((8), the gain in the filter pass-bands provided, respectively, by and increases and attains unity (0 dB). Simultaneously, a reduction of attenuation in the stopbands is also observed. Nevertheless, rather than the response of shown in Figure 15, there is no effect on the phase response corresponding to and , which remains zero radians for any values of frequency. On the other hand, ((8) cannot be applied in a cascade implementation because of the ripple that occurs below and around 0 dB in the pass-bands (see stopband detail in Figure 10), which can compromise the gain uniformity along the pass-bands. As shows no ripple along the pass-bands and no phase distortion, these characteristics allow the application of in a cascade implementation ((15) to improve the attenuation in the stopbands. , in turn, cannot be applied within a cascade implementation without causing some distortion in the signal phase alike.

Figures 9–12 reveal that the iterative application of time-domain averaging, according to (18), represents an alternative to using higher values of , since it permits to obtain narrower comb filter stopbands as well as phase distortion approximately zero along the filter pass-bands. Thus, in addition to enhancing the comb filter selectivity, time-domain averaging applied with a number of iterations enables to use a smaller number of averaging periods. In Figure 16, we illustrate the variation of the number of averaging epochs according to the number of iterations, which was calculated by taking into consideration a certain bandwidth (18.01 Hz) in the stopbands around −3 dB and .

Figure 16: Variation of the number of averaging epochs () according to the number of iterations () in (18), taking into account a certain bandwidth (18.01 Hz) in the pass-bands around −3 dB.  kHz.

For instance, by using and averaging epochs, it was possible to produce the same bandwidth in the stopbands around −3 dB by setting and in (18). Therefore, when J ≠ 1 in (18), the attenuation of the noise component is carried out by a factor higher than [8].

Figure 17 shows the difference amongst , , , and regarding the bandwidth in the stopbands around −3 dB.

(a)

(b)

(c)

(d)

Figure 17: Bandwidth of the stopbands around −3 dB (), for (a) , (b) , (c) , and (d) (). kHz.

The analysis depicted in this figure has taken into account . In case of , above around , the bandwidth around −3 dB is equal for all stopbands. In turn, for (and ), such a bandwidth is approximately similar for all stopbands, irrespective of the value of (and ). On the other hand, for , the width of the stopbands around −3 dB depends on the frequency of the stopband as well as the value of . The attenuation in the stopbands (as well as the gain in the pass-bands) may lie below −3 dB for smaller values of (see Figure 6). As the gain in the pass-bands approximates to 0 dB when increases, the bandwidth of the stopbands around −3 dB becomes narrower. However, it is not the same for all stopbands, but it enlarges according to the fundamental frequency (fd) of the stopband (Figure 17(c)). Such a difference is even higher for smaller values of . Therefore, as shown in Figure 17(c), the higher the fundamental frequency of the stopband, the broader is its respective bandwidth. This characteristic has been demonstrated to be useful for obtaining a more effective attenuation in the frequency bins associated with the harmonic gradient artefact spectral lines, with similar preservation of the EEG signal than that provided by the AAS method which is based upon time-domain averaging [20]. Thereby, broadening of the gradient artefact spectral lines provoked by micromovements of the subject head with the fMRI scanner, which are mostly observed in higher frequencies, can be more effectively accounted for by because of the increasing stopband bandwidths around −3 dB according to fd.

In Figure 18, an exemplary setting of parameters r, N, J, and L are used to illustrate the bandwidth of the stopbands according to the fundamental frequency (fd). As predicted by (16), the stopband bandwidth is further enlarged when is applied within the cascade implementation (), especially in stopbands with higher fundamental frequency (Figures 17(d) and 18). As can also be observed in Figure 18, , , and do not provide a stopband around 0 Hz, whereas it is produced by and .

Figure 18: Bandwidth of the stopbands around −3 dB () for an exemplary setting of parameters r, N, J, and L.  kHz.

Figure 19 depicts the impulse response associated with , , , and , for some values of the parameters r, N, J, and L. After the impulse is applied, an overshoot occurs in each of these responses, followed by decaying peaks spaced at the window length . The peak overshoot as well as the duration of the decaying peaks depend on the value of the parameter: smaller r, N, and J produce an increase in the overshoot and decreased settling time. Rather, a higher value of provides an increase in the overshoot and increased settling time.

(a)

(b)

(c)

(d)

Figure 19: Impulse response (), for (a) , (b) , (c) , and (d) ().

Regarding the computational effort, we verified that and are less computationally demanding than , , and . This can be ascribed to the computational efficiency of (11), which precisely requires the double of the computations of (2), as indicated in Figures 1 and 3. The worst case of computational effort was for , whose application was seriously compromised because of expensive computational memory demand by higher values of [6, 7].

Since a comb filter can be also implemented as a combination of notch filters, the selectivity provided by (21) could be used to design a highly selective comb filtering approach with a variable bandwidth for the different filter stopbands. Regarding the use of values of M > 2 in (20), we observed that some ripple around and below 0 dB may appear in the frequency response of for smaller values of , which can cause nonuniform gain along the pass-bands within the cascade implementation indicated in (21). Case higher values of are used in this scenario, the filter gain becomes uniform, but the notch filter stopband would have narrower bandwidth than in (21).

In future work, (20) should have its performance compared with other notch filter approaches [21], as well as (21) should be used and evaluated in applications where notch or comb filtering is required, such as power line rejection in biomedical signals [1]. Additionally, other filters in (9) should be investigated and evaluated for the iterative comb filtering implementation proposed here. Equations (15) and (16) should also be used and have their performance evaluated in signal processing applications other than the suppression of the gradient artefact from the EEG recordings [19, 20], where broadening of harmonic spectral lines is observed. As a further suggestion for future work, the iterative application of time-domain averaging should be assessed by using other kinds of moving-averaging filters, such as the exponential averaging and running averaging [2]. Last, the novel comb and notch filtering implementation described in this work shall be used and evaluated in other biomedical signal processing applications, such as speech signal processing and during estimation of evoked potentials (EPs) and event-related potential (ERPs) responses [1, 22].

5. Conclusions

A number of biomedical and other signal processing applications require the use of comb filtering approaches that perform elimination or extraction of periodic waveforms with a high degree of selectivity. As regards the elimination of periodic waveforms, the comb filter should be able to suppress the harmonics associated with the periodic signal and, simultaneously, to preserve the stochastic component or noisy signal according to the level of quality required by the application. Thus, it is important to make sure that the performance of the comb filtering approach meets the selectivity requirements of the application. Since such requirements are not always met by the existing methods, investigation and proposal of novel approaches to improve the comb filtering selectivity have been often described in the literature.

In this work, we have demonstrated how iterative filtering can be used to improve the selectivity of comb filtering approaches as well as to design a novel notch filter, which are based on the conventional FIR moving average filter. The novel comb filtering approach implementation is revealed to provide unity gain in the pass-bands, no effects on the signal phase, and broadening of the stopband bandwidth around −3 dB according to the fundamental frequency of the stopband. This characteristic has been proven to be useful within a scenario of broadening of spectral lines, such as that observed during the occurrence of the gradient artefact in the EEG signal recorded simultaneously with fMRI data. Moreover, a cascade implementation of the proposed approach permits to further increase the attenuation provided in the stopbands. In parallel, the iterative application of time-domain averaging allows using a smaller number of averaging epochs in order to estimate the periodic signal component. Hence, when a number of iterations are taken into consideration, the noise component can be attenuated with a factor higher than the square roots of the number of averaging epochs. Last, the novel notch filter implementation by iterative decomposition shows to provide a selective filtering with a large range of values for the stopband bandwidth around −3 dB.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this article and regarding the funding that they have received.

Acknowledgments

The authors are grateful to René M. H. Besseling for kindly providing review support. This work has been made possible by a grant from the European Union and Erasmus Mundus – EBW II Project and by a grant from CNPq – Science without Borders Program.

References

R. M. Rangayyan, Biomedical Signal Analysis: A Case-Study Approach, Wiley, New York, 2002.

S. Braun, “The synchronous (time domain) average revisited,” Mechanical Systems and Signal Processing, vol. 25, no. 4, pp. 1087–1102, 2011.

J. Meher, P. Meher, and G. Dash, “Improved comb filter based approach for effective prediction of protein coding regions in DNA sequences,” Journal of Signal and Information Processing, vol. 2, no. 02, pp. 88–99, 2011.

J. G. Proakis and D. G. Manolakis, Digital Signal Processing: Principles, Algorithms, and Applications, Prentice-Hall, Upper Saddle River, 3rd edition, 1996.

Z. Zeng, Y. Xie, Y. Wang, Y. Guan, L. Li, and X. Zhang, “An improved harmonic current detection method based on parallel active power filter,” in Proceedings of the 2nd Asia Conference on Power and Electrical Engineering, ACPEE 2017, pp. 1–6, IOP Publishing, Shanghai, China, March 2017.

M. Tahir and S. K. Mazumder, “Improving dynamic response of active harmonic compensator using digital comb filter,” IEEE Journal of Emerging and Selected Topics in Power Electronics, vol. 2, no. 4, pp. 994–1002, 2014.

A. Prodić, J. Chen, D. Maksimović, and R. W. Erickson, “Self-tuning digitally controlled low-harmonic rectifier having fast dynamic response,” IEEE Transactions on Power Electronics, vol. 18, no. 1, pp. 420–428, 2003.

S. Braun, “The extraction of periodic waveforms by time domain averaging,” Acustica, vol. 32, no. 2, pp. 69–77, 1975.

P. D. McFadden, “A revised model for the extraction of periodic waveforms by time domain averaging,” Mechanical Systems and Signal Processing, vol. 1, no. 1, pp. 83–95, 1987.

P. J. Allen, O. Josephs, and R. Turner, “A method for removing imaging artifact from continuous EEG recorded during functional MRI,” NeuroImage, vol. 12, no. 2, pp. 230–239, 2000.

R. Becker, P. Ritter, M. Moosmann, and A. Villringer, “Visual evoked potentials recovered from fMRI scan periods,” Human Brain Mapping, vol. 26, no. 3, pp. 221–230, 2005.

K. Anami, T. Mori, F. Tanaka et al., “Stepping stone sampling for retrieving artifact-free electroencephalogram during functional magnetic resonance imaging,” NeuroImage, vol. 19, no. 2, pp. 281–295, 2003.

W. X. Yan, K. J. Mullinger, M. J. Brookes, and R. Bowtell, “Understanding gradient artefacts in simultaneous EEG/fMRI,” NeuroImage, vol. 46, no. 2, pp. 459–471, 2009.

H. Mandelkow, P. Halder, P. Boesiger, and D. Brandeis, “Synchronisation facilitates removal of MRI artefacts from concurrent EEG recordings and increases usable bandwidth,” NeuroImage, vol. 32, no. 3, pp. 1120–1126, 2006.

G. S. Spencer, EEG-fMRI: Novel Methods for Gradient Artefact Correction [Ph.D. Thesis], University of Nottingham, Nottingham, UK, 2015.

L. Lin, Y. Wang, and H. Zhou, “Iterative filtering as an alternative algorithm for empirical mode decomposition,” Advances in Adaptive Data Analysis, vol. 1, no. 4, pp. 543–560, 2009.

Y. Wang, G.-W. Wei, and S. Yang, “Iterative filtering decomposition based on local spectral evolution kernel,” Journal of Scientific Computing, vol. 50, no. 3, pp. 629–664, 2012.

N. E. Huang, Z. Shen, S. R. Long et al., “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 454, no. 1971, pp. 903–995, 1998.

J. L. Ferreira, R. M. Aarts, and P. J. M. Cluitmans, “Optimized moving-average filtering for gradient artefact correction during simultaneous EEG-fMRI,” in Proceedings of the 5th ISSNIP-IEEE Biosignals and Biorobotics Conference (2014): Biosignals and Robotics for Better and Safer Living (BRC), pp. 1–6, IEEE, Salvador, Brazil, May 2014.

J. L. Ferreira, Y. Wu, R. M. H. Besseling, R. Lamerichs, and R. M. Aarts, “Gradient artefact correction and evaluation of the EEG recorded simultaneously with fMRI data using optimised moving-average,” Journal of Medical Engineering, vol. 2016, Article ID 9614323, 17 pages, 2016.

M. Vlček and P. Zahradník, “Fast analytical design algorithms for FIR notch filters,” IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 51, no. 3, pp. 608–623, 2004.

S. Sanei and J. A. Chambers, EEG Signal Processing, Wiley, Hoboken, 2007.