Substituting noise with zeros

I have a segment of data that has a significant amount of noise. I know the start and end point of the noise when it is graphed and I want to substitute these aberrant values with zeroes OR split the segment into two areas, deleting the noise in the middle.

[Star Strider]

Before you throw out parts of your data, first determine if the noise if frequency-limited or broadband. If it’s frequency-limited, you can likely filter it, retaining your underlying data.

If it’s broadband impulse noise a frequency-selective filter won’t work, but a Savitzky-Golay filter could work.

You can find out easily enough by taking the fft. (Note specifically the code between the top two plot figures.)

[atek]

It is not filter-able, the noise is coming from overloading the amplifier. I think I need to index it out but am not too sure how to go about doing so

[Star Strider]

If you can segregate the noise segments, I would replace them with NaN instead of zero. Since the noise came from railing your amplifier, those are likely not good data anyway.

We can’t help without actually having your data to work with. To upload it here, use the ‘paperclip’ icon and complete both the ‘Choose file’ and ‘Attach file’ steps.

[atek]

I was thinking of substituting with NaN's but that would interfere with the spectral analysis that I have been given code for. Substituting with NaN's makes more sense to me, but the data set is so large that I do not believe substituting with zeros will make much of a difference. I will upload the data shortly, thanks!

[atek]

Star Strider please see attached for the data I am having trouble with. Thank you for your help so far

[dpb]

The use of 0 will introduce discontinuities at the points where the signal breaks are; how significant that might be would depend largely on what the actual patterns are and to a lesser extent how long the overload sections are. I'd be concerned as well from a practical standpoint of how long after return from saturation it is before the sensor, whatever it is, is actually responding without bias or phase issues. Some are fairly tolerant to such transients others "not so much".

[atek]

I see your point - I had not considered whether the sensor has to come back to some functional baseline after saturation. Visualizing the data though, it seems to be OK.

We could take the data segment and split it into three components 1) before noise 2) noise substituted for NaNs 3) after noise - does this make more sense?

[Image Analyst]

Andy, you're making us work blind. Please attach a screenshot of your graph so we know what we're dealing with. Code to plot your .mat file would also be good.

[dpb]

Not enough info to answer the question of what the effect of simply ignoring the overload section is. Everything depends on what is the actual measurement and what you intend to get out of it, which as IA says, you're mum about. Trying to answer these kinds of questions independent of context (and very detailed context of both the DUT and the measurement system plus the analysis intent) is simply not possible with any certainty. The answer is, "it all depends"...some systems/analyses will make essentially no difference, others may be totally meaningless to continue. The real result is likely somewhere in the middle; the question is where in that continuum might it lay.

[atek]

Understood - here is the screenshot of plot(PleaseClean), I had previously uploaded the data itself in response to Strider's comment, please let me know if that is not what you are looking for

Sorry for the delay, Matlab keeps crashing on my laptop because of insufficient disk space - I sincerely appreciate the help!

[atek]

The goal of this is to generate a power spectrum using welch, the data is coming from an EEG electrode. I'm not too concerned about which method is "best" - I appreciate the points you are making but my level is very novice and I am looking for some starting code to help index out the signal from the overloaded amp

[Star Strider]

You have three columns in your .mat file.

1. What are they?
2. What am I looking at?
3. How best to plot them?
4. What’s your sampling frequency?

[atek]

Thank you for spoon-feeding me these questions Strider

1. The three columns represent three EEG electrodes that are in very close proximity to one another

2. They are in units of milli volts

3. Power spectrum

4. 1024

Answer 1

[Star Strider]

The good news is that I got a filter that filtered your signal up to the spike, because that’s the easiest way to deal with the discontinuities.

The bad news is that there’s nothing in your data. I cannot identify any EKG patterns such as QRS complexes in your (extremely noisy) data. Yours is one of those ‘I wish we could have discussed this earlier’ situations. I’ve gone into detail about the vectorcardiogram loop and its projections on the various EKG leads in many earlier posts. You didn’t describe your experimental set-up or design, so I don’t know if you used a ‘ground’ (actually common reference) electrode to subtract-out the noise in your differential amplifier set-up. The problem with putting the EKG leads too close together is that there is nothing for the vectorcardiogram loop to project onto. I’m not even sure where you put them.

I used a low-pass version of my usual band-pass filter design for EKGs on your signal, and could recover nothing in the filtered signal that I could identify as anything resembling a normal EKG. In figure(5), I even tried to see if I could recover information from a spatial representation of the signal, but there really is nothing in your data. I gave it my best effort!

My code:

D = load('Andy Tek PleaseClean.mat');

s  = D.PleaseClean;
Fs = 1024;                                                              % Sampling Frequency (Hz)
Fn = Fs/2;                                                              % Nyquist Frequency
L  = size(s,1);
t = linspace(0, 1, L)/Fs;                                               % Sampling Time (s)

Smean = mean(s);
Sstd  = std(s);
CE95  = bsxfun(@plus, Smean, bsxfun(@times,Sstd*1.96,[-1  1]'));

[maxs,idx] = max(s);
sLen = min(idx-10);
s = s(1:sLen,:);
t = t(1:sLen);

figure(1)
subplot(3,1,1)
plot(t, s(:,1))
axis([xlim    0.001*CE95(:,1)'])
grid
subplot(3,1,2)
plot(t, s(:,2))
axis([xlim    0.001*CE95(:,2)'])
grid
subplot(3,1,3)
plot(t, s(:,3))
axis([xlim    CE95(:,3)'])
grid

Q1 = s(1:5,:);                                                          % Look

fts = fft(s)/sLen;
Fv = linspace(0, 1, fix(sLen/2)+1)*Fn;
Iv = 1:length(Fv);

figure(2)
subplot(3,1,1)
plot(Fv, abs(fts(Iv,1))*2)
axis([0  100    0  1E-1])
grid
subplot(3,1,2)
plot(Fv, abs(fts(Iv,2))*2)
axis([0  100    0  1E-1])
grid
subplot(3,1,3)
plot(Fv, abs(fts(Iv,3))*2)
axis([0  100    0  1E-4])
grid

Wp =  50/Fn;
Ws =  75/Fn;
Rp = 10;
Rs = 50;
[n,Wn] = buttord(Wp, Ws, Rp, Rs);
[b,a] = butter(n,Wn);
[sos,g] = tf2sos(b,a);

figure(3)
freqz(sos, 2048, Fs)

sf = filtfilt(sos,g,s);
Q2 = sf(1:5,:);                                                         % Look

figure(4)
subplot(3,1,1)
plot(t, sf(:,1))
axis([xlim    0.001*CE95(:,1)'])
grid
subplot(3,1,2)
plot(t, sf(:,2))
axis([xlim    0.001*CE95(:,2)'])
grid
subplot(3,1,3)
plot(t, sf(:,3))
axis([xlim    CE95(:,3)'])
grid

figure(5)
plot3(sf(:,1), sf(:,2), sf(:,3))
grid on

[atek]

Thank you all! Strider - I was a bit disengenuous to say it was EEG, it is local field potential data from a neuronal probe. So that should explain the lack of signal relative to what you expected? I was not expecting such detailed response! Thank you very much

[Star Strider]

My pleasure!

That would definitely explain it. To do any such recording, you have to have a specific reference electrode, and as many other spatial electrodes as your instrumentation can accommodate. You also have to have them arranged in a particular pattern so that you can make sense of the data.

As with all of us here, I do my best to help. For me, that is particularly with respect to biomedical instrumentation and signal processing problems.

[atek]

This is a really great community - hopefully I'll be able to contribute solutions in a couple years. Thanks again Strider!

[Star Strider]

As always, my pleasure!

[atek]

Star Strider, I have a request. Could you show me how I would replace all sample points that are outside of the CI with the mean value for the segment?

[Image Analyst]

Same as what I showed below, you just use conv() instead of medfilt1().

meanSignal = conv(signal, ones(1, windowWidth)/windowWidth, 'same');

[atek]

Thanks for the input Image Analysis.

windowWith here would represent what? The dimensions of the matrix? My data set is [825634,3]

For my own curiosity I was trying to use bsxfun to accomplish swapping out the values above 3SD for the mean, would you be able to advise on that?

Thank you for your help!

[Image Analyst]

windowWidth is the moving window size over which you want to take the mean inside. For example if you had a 1-D array and a window width of 5, the element #13 of the output would be the average of elements 11, 12,13,14, & 15 of the input because those are the 5 elements enclosed by the window when the window is centered at element 13.

[atek]

So I want to have a 30 second window width and fs is 1024. My signal is 'PleaseClean' which has dimensions [825634,3]. I ran the following:

windowWidth = 1024*30 meanSignal = conv(PleaseClean, ones(1, windowWidth)/windowWidth, 'same');

But I received the following error Error using conv (line 28) A and B must be vectors.

Answer 2

[Image Analyst]

Can you just do a salt and pepper noise cleaning? Just run the signal through a median filter, and where there are huge spikes, replace them with the median filtered version

badElements = abs(signal) > 5; % or whatever constitutes noise.
% Get rid of spikes.  However, this alters the "good" elements also.
medianFilteredSignal = medfilt1(signal, 3);
% Replace only the bad elements, not the good ones.
signal(badElements) = medianFilteredSignal(badElements);