3
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • 3
    $\begingroup$ what in your plot is signal, what is artifact? $\endgroup$
    – Marcus Müller
    Oct 27, 2016 at 20:42
  • $\begingroup$ I've updated the question. The artifacts are those high-amplitude biphasic pulses. $\endgroup$
    – bluevoxel
    Oct 28, 2016 at 12:48
  • 1
    $\begingroup$ 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.) $\endgroup$
    – robert bristow-johnson
    Oct 29, 2016 at 2:59

2 Answers 2

Reset to default
1
$\begingroup$

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

Improve this answer
$\endgroup$
12
  • 1
    $\begingroup$ uhm, so you're subtracting p after adding it. i don't think that speaks to the original problem. $\endgroup$
    – robert bristow-johnson
    Oct 29, 2016 at 2:56
  • 1
    $\begingroup$ what if you only have x and you don't have p? $\endgroup$
    – robert bristow-johnson
    Oct 29, 2016 at 3:06
  • 2
    $\begingroup$ 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? $\endgroup$
    – robert bristow-johnson
    Oct 29, 2016 at 3:50
  • 1
    $\begingroup$ "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... $\endgroup$
    – msm
    Oct 29, 2016 at 4:15
  • 2
    $\begingroup$ 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. $\endgroup$
    – robert bristow-johnson
    Oct 29, 2016 at 4:23
1
$\begingroup$

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

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – a concerned citizen
    May 21, 2017 at 6:18
  • $\begingroup$ @aconcernedcitizen Note that bandwidth is tunable with the alpha factor provided in the implementation (see dsp.stackexchange.com/questions/31028/…). 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). $\endgroup$
    – Dan Boschen
    May 21, 2017 at 10:45
  • $\begingroup$ 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). $\endgroup$
    – Dan Boschen
    May 21, 2017 at 10:49
  • $\begingroup$ 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. $\endgroup$
    – a concerned citizen
    May 21, 2017 at 16:03
  • 1
    $\begingroup$ Nice! And 15 more characters so I can say nice! $\endgroup$
    – Dan Boschen
    May 21, 2017 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.