Divide time series into regions

Hello, i have a graph with rapidly changes and i need to divide it into regions like are shown in the picture

I can divide it considered the amplitude of each regions but that's not so appropreate because i have also big dicreases like in region B that reaches the previous region, and that's something that troubles me. I thought maybe i could use standard deviation between local peaks and local minima. What are youur thoughts and do you know any other simple (or not so) methods ? Thanks in advance. Note that this graph sould be divided unsmoothed , and then only each region will be processed.

Answer 1

[Star Strider]

I don’t have your data so I can only describe an approach. I would use a low-pass filter with a very low frequency passband to isolate the d-c offset for each section. (Ideally, you should get a series of ‘steps’ as the output.) Then use those values to define the different regions. You will probably have to devise a way to deal with the transitions between the steps to define them correctly.

[Manolis Michailidis]

If that helps the y axis are frequencies so for instance for the region A the low-pass frequency 2500 Hz and bandwith 200 Hz are good values? I am not mush experienced in filters. Anyway thanks i will give it a try and will let you know , thanks.

[Star Strider]

I would filter the y-value as the signal, with respect to the x-value as time.

Begin by taking the fft of your data to determine the signal frequency you want to exclude. (Use the code between the first two plots as your guide.) This will be all your data, since you want to keep the d-c offset only. (My filter design strategy is here.)

If you upload a sample of your data, I can see if my idea works with it.

[Manolis Michailidis]

ok i tried to upload here but i get an error so here is my data http://www40.zippyshare.com/v/JBmDamN0/file.html and as far the sampling frequency well i don't know how to compute it, because this graph is the trajectory of the fundamental frequency of a signal which fs=44.1kHz

[Star Strider]

I need to be clear on the details of what the x and y data are, and their units. Do all your data have the same x sampling frequency?

[Manolis Michailidis]

Ok let me explain a bit. Using autocorrelation i am tracking the fundamental frequency of a signal, all my signals have Fs=44.1kHz. So i am ploting the fo and the y axis is frequency in Hz and the x time in seconds always. I am uploading another mat file with the time vector. http://www71.zippyshare.com/v/3vWYH3n6/file.html So

plot(ts,fo)

will show you the graph, i hope it helped.

[Manolis Michailidis]

Alternativly i tryied to compute the dc taking the mean for each pair of consecutive peaks and then subtructing the fo from this result. I plotted them and well another non working solution , as you can see that the regions can't be seppareted.

Anyway , thanks for your anwers and effort.

[Star Strider]

Was working on it. This is not a trivial problem!

My idea essentially worked:

D = load('Manolis Michailidis fo.mat');
fo = D.fo;
x = 0:length(fo)-1;

figure(1)
plot(x, fo)
grid

Ts = 1;                                         % Sampling Time Interval
Fs = 1/Ts;                                      % Sampling Frequency
Fn = Fs/2;                                      % Nyquist Frequency
Ffo = fft(fo)/length(fo);                       % Fourier Transform
Fv = linspace(0, 1, fix(length(Ffo)/2)+1)*Fn;   % Frequency Vector
Iv = 1:length(Fv);                              % Index Vector

figure(2)
plot(Fv, abs(Ffo(Iv)))
grid
axis([0  0.01    ylim])

Wp = 0.005/Fn;                                  % Normalised Passband Frequency
Ws = Wp/0.8;                                    % Normalised Stopband Frequency
Rp =  1;                                        % Passband Ripple
Rs = 10;                                        % Stopband Ripple
[n,Ws] = cheb2ord(Wp, Ws, Rp, Rs);              % Design Chebyshev Type II LPF
[b,a]  = cheby2(n,Rs,Ws,'low');
[sos,g] = tf2sos(b,a);                          % Use ‘SOS’ For Stability
fof = filtfilt(sos,g,fo);                       % Filter Signal

xidx = find((fof >= 2900) & (fof <= 3100));     % Interpolate To Find ‘x’ For ‘fo=3000’
xbrk = [find(diff([0 xidx]) > 100) length(xidx)];
for k1 = 1:length(xbrk)-1
    idxrng = xidx(xbrk(k1):xbrk(k1+1)-1);
    xt = interp1(fof(idxrng), x(idxrng), 3000, 'linear','extrap');
    segx(k1) = find(x <= xt, 1, 'first') + idxrng(1)-1;
end

figure(3)
plot(x, fo)
hold on
plot(segx, 3000*ones(length(segx)), '+r', 'MarkerSize',15)
hold off
grid

The ‘segx’ vector are the x-index ‘breakpoints’ in your data, based on using an ‘fo’ value of 3000 to distinguish the segments. Using the low-pass-filtered signal made defining the segments much easier.

It might be possible to make this more efficient, but I opted to leave the analysis section of the code in it, since not all of your data may have the same characteristics.

I finished this before I saw your uploaded time vector, so replace the ‘x’ vector in my code with your time vector. The ‘Ts’ assignment then becomes:

Ts = mean(diff(t));

where ‘t’ is your time vector. You will have to define the numerator in the ‘Wp’ assignment to correspond to the correct frequency derived from your time vector, but the code should otherwise work.

[Manolis Michailidis]

Thank you very much, please edit it to your previous answer so i can accept it.

[Star Strider]

My pleasure.

What do you want me to edit?

[Manolis Michailidis]

it's all fine , thanks again

[Star Strider]

O.K.

My pleasure.

Answer 2

[Ashim]

You can also use alternatively the findchangepts() command from Matlab to find the different regions> here again, there are multiple options based on std dev or mean, or rms, or threshold. Try it out

[Star Strider]

Note that the findchangepts (link) function was introduced in R2016a.

[Manolis Michailidis]

I finished my project but anyway i will test also this function just for curiocity, thank you again