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.
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)
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.
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.
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)
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)
(a)
(b)
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)
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)
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 11 depicts the frequency response provided by (17), taking into account and , for (dark trace), (green trace), and (gray trace).
(a)
(b)
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.
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.
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.
4. Further Comparative Analyses and Discussion
Figure 15 shows the frequency response of (3), for some values of and .
(a)
(b)
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 .
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)
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 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)
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.