3

How can I remove a biphasic, square-wave-like artifacts from the signal, so I'll keep as much as possible information "under" the artifact? Below is the signal with artifacts:

enter image description here

The artifacts are those two big spike-like pulses. They are caused by square-wave pulses applied to the skin in neighborhood of recording electrodes.

Peter K.
  • 25,714
  • 9
  • 46
  • 91
bluevoxel
  • 201
  • 1
  • 2
  • 6
  • 3
    what in your plot is signal, what is artifact? – Marcus Müller Oct 27 '16 at 20:42
  • I've updated the question. The artifacts are those high-amplitude biphasic pulses. – bluevoxel Oct 28 '16 at 12:48
  • 1
    are the interfering pulses added at regular intervals of time? is the interfering signal a "buzz" with a pitch or fundamental frequency? if so, a pitch-detector to determine the pitch of it, and a subtractive comb-filter that is tracking that frequency is the only automated way i can think of to do it. (perhaps an LMS adaptive filter can be cooked up to do it.) – robert bristow-johnson Oct 29 '16 at 02:59

2 Answers2

1

Maybe you find this useful:

N=500;                          %signal length  
t=0:N-1;
s=20*randn(size(t));
x=s;
p=200*[ones(1,5) -ones(1,5)];   %square artifact template
x(:,101:110)=x(:,101:110)+p;    %add artifact to signal at two random points
x(:,401:410)=x(:,401:410)+p;
figure;plot(x)
threshold = 1e5;                %peak detection threshold
c=xcorr(p,x);                   %cross-correlation
[pk l]=findpeaks(c,'MinPeakHeight',threshold);
l=N-l+1;                        %find the actual lag
y=x;
for j=1:numel(l)                %remove artifact
y(:,l(j):(l(j)+length(p)-1))=y(:,l(j):(l(j)+length(p)-1))-p;
end
hold on;
plot(y,'r--')
legend('corrupted signal','recovered signal')

enter image description here

msm
  • 4,285
  • 1
  • 13
  • 25
  • 1
    uhm, so you're subtracting p after adding it. i don't think that speaks to the original problem. – robert bristow-johnson Oct 29 '16 at 02:56
  • 1
    what if you only have x and you don't have p? – robert bristow-johnson Oct 29 '16 at 03:06
  • We never have p. We only have x. It is our task to find (estimate) p first. From the curve you can see it is not so difficult to do so. Here, it is implicit that p is perfectly estimated... BTW, "artifact" is something that we usually know about. – msm Oct 29 '16 at 03:13
  • 2
    your code does not reflect that at all. you construct p. you add it to x, you do correlation to find where you added it to x. then you subtract it. although you do something to find where the interfering pulse is, you do nothing to determine its width nor its amplitude (nor shape). you just subtract it. how would you know how big of a thing to subtract? – robert bristow-johnson Oct 29 '16 at 03:50
  • 1
    "The artifacts are those two big spike-like pulses. They are caused by square-wave pulses applied to the skin in neighborhood of recording electrodes." I am sure the OP knows what square pulse s/he is applying on the skin... – msm Oct 29 '16 at 04:15
  • 2
    well, we'll let the OP decide that. in "cleaning up" a signal, i generally understand that to mean we know something about the crap dirtying up the signal, but we don't know exactly what the crap is, otherwise we can just subtract it. – robert bristow-johnson Oct 29 '16 at 04:23
  • 1
    I think that there might be some confusion stemming from the fact that the code actually simulates a signal, then corrupts it with the pulse sequence and then recovers it with a threshold (?) Maybe the basis of the technique needs to be described better with words, i.e use a threshold to recover the disturbance (because it is orders of magnitude larger) and then subtract it from the main signal. Provided that this additive model is valid in this setting of course (?) – A_A Oct 29 '16 at 09:13
  • @msm The problem is, due to electrode-specific impedance, artifacts are different across the electrodes. Thats why it's unknown to me how p can be estimated. I can for example create an average of the artifact based on its all occurences in the signal and than use it as p. But it will be only an approximation. – bluevoxel Oct 30 '16 at 20:08
  • @bluevoxel Approximation is your only choice! Even if you use more advanced estimation of 'p', it will be only an approximation. But I think averaging as you suggested works here very well. – msm Oct 30 '16 at 23:56
  • @msm Not always. Substraction of averaged representation of the artifact from single, individual artifacts sometimes leave ugly spike-like post-artifacts. I tried to increase sampling frequency of the signal in order to reduce this naughty effect, but it is not a perfect solution. – bluevoxel Oct 31 '16 at 14:15
  • Ignoring the filtering impact of skin on your pulse is a mistake. I don't know how you performed it and what your approach in 'averaging' those pulses is. But the subtraction should be done within the same durations corresponding to the up- and downward peaks (with appropriate scaling). Of course if the duration is not appropriate, then there are some spikes (difference of the two duration). If you can't benefit from the known duration of the pulse in any sense, I would suggest assuming a parametric pulse whose four parameters are the two durations and the two amplitudes. @bluevoxel – msm Oct 31 '16 at 22:44
  • @bluevoxel You can use the code I gave you to locate the pulses. Then use any technique to find the four parameters (e.g. look for the zero-crossing in the given duration...), adjust the pulse using the acquired parameters, and subtract it. – msm Oct 31 '16 at 22:50
1

Another answer that may be useful due to its simplicity: if the spikes are occurring at a consistent rate, a moving average filter of length T where 1/T is the fundamental repetition rate of the spike, will have a null in frequency conveniently at the spike frequency and every harmonic of the spike (which is convenient given the spike will have many integer harmonics). Further and for the same reason, choose a sampling rate that is an integer multiple of the spike frequency to minimize folding artifacts (if you have that luxury). The moving average filter will also attenuate signal content (as a Sinc filter), so this may not be a desired solution. Optionally an interpolated exponential filter could be used with the same strategy except providing a narrow null just around the spike frequency and harmonic locations while minimizing attenuation of other frequencies. If there is interest I can provide specific details on that. Both approaches assume a consistent frequency of spikes occurring.

UPDATE: I detailed the "interpolated exponential filter" in this post: What kind of filter I should use to remove the oscillations in this autocorrelation function? showing the following example harmonic nulling filter that I believe would work quite well to remove the noise spikes if they do occur at a consistent repetition rate. This is an example response but the nulls would be placed at the repetition rate of the noise spikes in actual implementation, for that reason the signal would need to have a commensurate sampling rate to the spikes which can be done with resampling prior to implementing this filter.

harmonic nulling filter response

Dan Boschen
  • 50,942
  • 2
  • 57
  • 135
  • The zeroes at specific frequencies will have a very narrow bandwidth, so any jitter could compromise it, particularly for the magnitude graph you posted. At the cost of double the length, one can double the zeroes and place them symmetrically across the point of interest, while being very close together, such that the lobe between has sufficient attenuation, while providing some small enough bandwidth. I played with this for an active power filter, where I needed 50Hz+/-2Hz, it worked pretty well for a quick and dirty solution. – a concerned citizen May 21 '17 at 06:18
  • @aconcernedcitizen Note that bandwidth is tunable with the alpha factor provided in the implementation (see https://dsp.stackexchange.com/questions/31028/transfer-function-of-second-order-notch-filter/31030#31030). But cascading two (or more) with slightly offset zeros as you suggest is an interesting idea, thanks! Cascading if you did it that way is also robust in isolating the smaller order IIR filter sections similar to what is done in biquad implementations. (to minimize the possibility of issues with delay and precision in higher order IIR implementations). – Dan Boschen May 21 '17 at 10:45
  • But the harmonic approach above isn't limited to IIR filters-- the basic idea is to take any template low pass filter and then insert zeroes in between the coefficients to repeat the template at each frequency. (In the example I had inserted 10 zeros). – Dan Boschen May 21 '17 at 10:49
  • Don't forget the "quick and dirty" part, that's very important. :-) Also, I used that with a basic moving average (FIR). The resultant polynomial was very triangularish, with a slightly smoothed out top. I chose to do it because there wasn't much room left in the uC, otherwise I would have gone for an adaptive moving average, much lower delay. – a concerned citizen May 21 '17 at 16:03
  • 1
    Nice! And 15 more characters so I can say nice! – Dan Boschen May 21 '17 at 19:46