How do I count cycles of a periodic signal that contains noise?

I have a periodic signals that contains noise as shown below in graph 1 and 2. Graph 1 has 3 periods. Graph 2 has 5 periods. How do setup an algorithm that counts these periods and then is able to extract the 1st period and then the 2nd and so on?

graph 1

graph 2

Answer 1

[Star Strider]

There is no really elegant and ‘non-kludgy’ way, but having had to deal with similar situations, I suggest two options:

For Graph 1, the easiest is to use ‘diff’:

difsig1 = diff([0; sig1]):    % This assumes a column vector, the initial zero makes the indices of ‘difsig1’ match those of ‘sig1’

Characterize your noise limits as:

statsnoise = mean(difsig1)+1.96*[-std(difsig1)  std(difsig1)];

then use ‘find’ to determine the maxima and minima of ‘difsig1’ as being those outside the limits of ‘statsnoise’.

You can do essentially the same with Graph 2, although instead of using ‘diff’ and ‘find’, use ‘findpeaks’ to find the leading edge, then flip it

sig2flip = sig2-max(sig2);

and use ‘findpeaks’ again on ‘sig2flip’ to find the trailing edge.

These functions all use indices, and produce the indices of the requested data as their output. You will have to map these indices back to their respective times to calculate the periods, pulse widths, etc.

Answer 2

[Image Analyst]

Just threshold it. It looks like there is some value that separates the tops fo the signals and the bottoms.

highSignalLocations = signal > thresholdValue;
lowSignalLocations = signal <= thresholdValue;

You could use bwthresh() if you want, to find the threshold value. To count, if you have the Image Processing Toolbox you can use bwlabel

[highRegions numberOfHighRegions] = bwlabel(highSignalLocations);
[lowRegions numberOfLowRegions] = bwlabel(lowSignalLocations);
labeledRegions = lowRegions + highRegions;

numberOfRegions is the number of times you have a run above the threshold. Do you know that numberOfHighRegions will always be the same as numberOfLowRegions?

You can then extract out any of the segments easily, or all of them. Something like

region1 = signal(labeledRegions == 1);
region2 = signal(labeledRegions == 2);
region3 = signal(labeledRegions == 3);

and so on, up to numberOfRegions. You can put it in a for loop to get all the different length regions into cells of a cell array. If you need code for that, say so. But include your signal so I can test it on your actual data.

Answer 3

[Wayne King]

In R2012a, the Signal Processing Toolbox has functions specifically for the type of waveform you show: see pulseperiod, pulsesep, and pulsewidth