Separating data into different regions

Hi!

Would anyone have any ideas on how to obtain the final processed signal from the raw data in the attached image?I tried setting thresholds to separate the signals the raw signal into the three regions,but the signal is too noisy. I tried filtering it but that introduces overshoots at the transition and adds to the difficulty to separate the signal. Any pointers on how to approach this problem would be greatly appreciated! Thanks!

Answer 1

[Star Strider]

Without your signal (and necessary parameters like sampling frequency), it’s impossible to define a specific approach, much less write code.

My approach would be to bandpass (or lowpass) filter it, then perhaps use the Signal Processing Toolbox envelope function (in R2015b and later) or equivalently, the absolute value of the Hilbert transform of your bandpass (or lowpass) filtered signal, or the Savitzky-Golay filter to define your filtered signal. You can then just threshold it.

[George Mathai]

Thanks, Star! I have attached the data which was sampled every second. I managed to use the Savitsky Golay filter to get some pretty good results except in a couple of regions where it overshoots the raw data. Thanks so much for you help! Will try envelope and the Hilbert transform.

[Image Analyst]

I can't really tell the difference between brown and red with these single pixel wide lines on my monitor. Maybe use magenta or yellow or green or something. I still think it would be best to throw away elements that are obviously noise before smoothing with sgolayfilt(). Are there any elements like that, or ALL of those spiky things valid signal?

[Star Strider]

This is my approach, using only a lowpass filter and simple thresholding:

D = load('George Mathai data_mat.mat');
sig = D.data_mat;
Lsig = length(sig);
Fs = 1;                                                         % Sampling Frequency
Fn = Fs/2;                                                      % Nyquist Frequency
Ts = 1/Fs;                                                      % Sampling Time
Tv = linspace(0, Lsig, Lsig)*Ts;                                % Time Vector                                                     

ftsig = fft(sig)/Lsig;                                          % FFT
Fv = linspace(0,1, fix(Lsig/2)+1)*Fn;                           % Frequency Vector
Iv = 1:length(Fv);                                              % Index Vector

figure(1)
plot(Fv, abs(ftsig(Iv)))
grid

Wp = [0.02]/Fn;                                                 % Design Lowpass Filter
Ws = [1.5].*Wp;
Rs = 25;
Rp =  1;
[n,Wn] = buttord(Wp, Ws, Rp, Rs);
[b,a] = butter(n, Wn);
[sos,g] = tf2sos(b,a);                                          % Comvert To ‘SOS’ For Stability

figure(2)
freqz(sos, 2048, Fs)                                            % Filter Bode Plot

fsig = filtfilt(sos,g,sig);                                     % Filter Signal
gt_lv = fsig > 5;                                               % Threshold Signal (Empirical Threshold)                                              

figure(3)
plot(Tv, sig, '-b', 'LineWidth',0.5)
hold on
plot(Tv, fsig, '-r', 'LineWidth',1)
plot(Tv, gt_lv*90, '-g', 'LineWidth',1)
hold off
grid
axis([0  1000    -5  100])
xlabel('Time (sec)')
ylabel('Amplitude (units)')
legend('Original Signal', 'Filtered Signal', 'Thresholded Signal', 'Location', 'NW')

You can delete the fft call through the code for figure(1). I kept them in so you can see my filter design strategy.

The code produces this plot for figure(3):

I am not certain how best to get the small isolated threshold regions incorporated into the larger ones. You will need to find out how best to include them based on the separation criteria you define. One possibility would be to have a minimum value for the width of a particular threshold region (in my code, the ‘gt_lv’ assignment defines the thresholds), and then combine it with the others if its width is less than a particular value by setting the gap between the gap between it and the one following it to a logical value of 1. Defining them based on the width of the gap between the threshold areas is not itself a sufficient criterion, since it merges all adjacent ones. I experimented with a number of ways to define them.

[Image Analyst]

If you have the Image Processing Toolbox, you can combine the smaller regions into a neighboring region (which is hopefully larger) by using bwareaopen(). Something like

minSeparation = 10; % Or whatever distance from neighbor you want
gt_lv = ~bwareaopen(~gt_lv, minSeparation);

However if two big regions are closer than that, it will merge them. So you might have to first check the length of the regions before merging. You can do that with regionprops(), but then it's starting to get a bit complicated and the "separation criteria" that Star mentioned is likely going to have to be a full blown function with lots of tests and stuff, not a simple one-liner like I did.

Answer 2

[Image Analyst]

It's a little hard to tell what the signal is and what the noise is. But I guess the huge spikes are noise. If you just use a Savitzky-Golay filter with those included, then they will influence your result, which you do not want. It looks like the majority of the noise can be eliminated by just thresholding and setting to null:

threshold = 100; % whatever works.
% Find noise indexes:
indexesOfNoise = y > threshold;
% Remove x and y at those indexes by setting equal to null
y(indexesOfNoise) = [];
x(indexesOfNoise) = [];
% Now do sgolayfilt()
filteredSignal = sgolayfilt(y,k,f);

You might have to also remove spikes that go too low, like down to the x axis if those are not valid signal.

Perhaps you'd like to look at Dr. Abouali's File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/38841-matlab-implementation-of-harmonic-analysis-of-time-series--hants-