离散信号的周期性判定，C++实现

问题：给定一列离散信号，判定是否为周期信号
解决方案：
理论上，计算信号的自相关，如果是周期性信号，则其自相关序列依然为周期性信号切几乎不会衰减；否则，则会出现逐渐衰减至0。
实际情况下由于噪声的存在，偶尔自相关的最大值不出现在$\tau$=0处，而且如果$\tau$值取的比较大的时候即便是周期信号也会出现峰值衰减很大的情况，所以一般$\tau$值取N/2 N=len(vec)。具体理论请参见《数字信号处理》第四版P90的2.6.3小节。
算法的核心代码是自相关的计算，在实际编码的时候对书上的公式做了一个小小的处理，将$r_{yy}(l)=\frac{1}{M}{\sum }_{n=0}^{M-1}[x(n)x(n-l)]$中的常量M换成变量m，m={N, N-1, …, N/2}， 即按照Matlab的xcorr的’unbiased’模式计算。
在判断的时候，网上有资料给出的方案是取自相关序列峰值序列的中间值，然后判断下一峰值和这一值的比率，在0.9范围内即可。但是，实际中这做还是有问题，估计作者没有用实际采集到的数据做测试，实际数据中仅仅这样判断误判率很高，存在问题的自相关序列图没保存，就没法上图了，反正肯定有问题！我用了取自相关序列第一个值的[0.9,1.1]这个区间做依据，如果峰值序列的90%都在这个范围内即认为ok，这样也还是存在误判，但是比网上的资料效果要好点了。当然，还得结合实际情况，同时利用先验经验比如频率、峰值等一并判断。

CycleCheck.h

#pragma once

#include <vector>
#include <iostream>
using namespace std;

typedef struct WavePeakFeature
{
	vector<int> peaks_index;
	vector<float> peaks_value;
	vector<float> freq_list;	
}WavePeakFeature;

class CCycleCheck
{
public:
	CCycleCheck();
	CCycleCheck(int);
	~CCycleCheck();

	bool CycleCheck(vector<int>input_data);

	// sampling frequency
	int sample_freq;

public:
	static const int default = 0;
	static const int unbiased = 1;

private:
	// computer vec's autocorrelation	
	vector<float> xcorr(vector<int>vec, int model);

	// computer the difference for the vec[i]-vec[i-1]
	vector<int> diff(vector<int>vec);
	vector<float> diff(vector<float>vec);

	// return the sign(x)
	vector<int> sign(vector<int>vec);
	vector<int> sign(vector<float>vec);

	vector<int> median_filter(vector<int> vec, int window);

	WavePeakFeature find_peaks(vector<int> vec);
	WavePeakFeature find_peaks(vector<float> vec);

	bool is_cycle(WavePeakFeature &feature_parame);

	// filter wave peaks based on pulse wave characteristics
	// peaks_index: peak's index in original data vector
	// peaks_value: the value of the peaks	 
	WavePeakFeature* peaks_filter(WavePeakFeature &feature_parame);
};

CycleCheck.cpp

#include "stdafx.h"
#include "CycleCheck.h"
#include <math.h>
#include <numeric>
#include <algorithm>
#include <exception>

CCycleCheck::CCycleCheck()
{
}

CCycleCheck::CCycleCheck(int sample_freq)
{
	this->sample_freq = sample_freq;
}

CCycleCheck::~CCycleCheck()
{
}

bool CCycleCheck::CycleCheck(vector<int> input_data)
{	
	bool res = false;
	try
	{
		vector<int> data = median_filter(input_data, 3);
		vector<float> xc = xcorr(data, CCycleCheck::unbiased);		
		WavePeakFeature peaks_feature = find_peaks(xc);
		WavePeakFeature *features = peaks_filter(peaks_feature);
		res = is_cycle(*features);
	}
	catch (exception &e)
	{
		res = false;
	}
	return res;
}

bool CCycleCheck::is_cycle(WavePeakFeature &feature_parame)
{
	bool res = false;

	vector<float> peaks_value = feature_parame.peaks_value;	
	
	if (peaks_value.size() < 3)
	{
		return false;
	}

	// 理论上自相关的第一项值最大
	float base_value = peaks_value[0];
	float min = base_value * 0.9;
	float max = base_value * 1.1;
	float num = 0.0;

	// 峰值相关周期特性
	std::for_each(std::begin(peaks_value), std::end(peaks_value), [&](const double ele) {
		num += ele > min && ele < max ? 1 : 0;
	});
	
	if (num > peaks_value.size()*0.9 )
	{
		res = true;
	}
	return res;
}

WavePeakFeature CCycleCheck::find_peaks(vector<int> vec)
{
	WavePeakFeature peak_feature;
	vector<int> vec_sign = sign(diff(vec));
	vector<int> vec_diff = diff(vec_sign);
	for (int i = 0; i< vec_diff.size(); i++)
	{
		if (vec_diff[i] == -2)
		{
			peak_feature.peaks_index.push_back(i + 1);
			peak_feature.peaks_value.push_back(vec[i + 1]);
		}
	}
	return peak_feature;
}

WavePeakFeature CCycleCheck::find_peaks(vector<float> vec)
{
	WavePeakFeature peak_feature;
	vector<int> vec_sign = sign(diff(vec));
	vector<int> vec_diff = diff(vec_sign);
	for (int i = 0; i < vec_diff.size(); i++)
	{
		if (vec_diff[i] == -2)
		{
			peak_feature.peaks_index.push_back(i + 1);
			peak_feature.peaks_value.push_back(vec[i + 1]);
		}
	}
	return peak_feature;
}

vector<int> CCycleCheck::median_filter(vector<int> vec, int kernel_size)
{
	vector<int> median;
	vector<int>::iterator iter = vec.begin();
	for (int i = 0; i < vec.size() - kernel_size;i++)
	{
		vector<int> window(iter + i, iter + kernel_size);
		for (int m = 0; m < kernel_size; ++m)
		{
			int min = m;
			for (int n = m + 1; n < kernel_size; ++n)
				if (window[n] < window[min])
					min = n;			
			int temp = window[m];
			window[m] = window[min];
			window[min] = temp;
		}
		median.push_back(window[kernel_size / 2 + 1]);
	}
	return median;
}

vector<float> CCycleCheck::xcorr(vector<int> vec, int model)
{
	int vec_len = vec.size();
	vector<float> res;
	vector<int> zero_vec(vec_len / 2, 0);
	vector<int> padding_vec;
	padding_vec.insert(padding_vec.end(), vec.begin(), vec.end());
	padding_vec.insert(padding_vec.end(), zero_vec.begin(), zero_vec.end());
	
	if (model == CCycleCheck::unbiased)
	{
		for (int i = 0; i < vec_len / 2; i++)
		{
			float sum = 0;
			for (int j = 0; j < vec_len; j++)
			{
				sum += vec[j] * padding_vec[i + j];
			}
			res.push_back(sum / (vec_len - i));
		}
	}
	else
	{		
		for (int i = 0; i < vec_len / 2; i++)
		{
			float sum = 0;
			for (int j = 0; j < vec_len; j++)
			{
				sum += vec[j] * padding_vec[i + j];
			}
			res.push_back(sum);
		}		
	}
	return res;
}

// computer the difference for the vec[i]-vec[i-1]
vector<int> CCycleCheck::diff(vector<int>vec)
{
	vector<int> res;
	for (int i = 1; i < vec.size(); i++)
	{
		res.push_back(vec[i] - vec[i - 1]);
	}
	return res;
}

vector<float> CCycleCheck::diff(vector<float>vec)
{
	vector<float> res;
	for (int i = 1; i < vec.size(); i++)
	{
		res.push_back(vec[i] - vec[i - 1]);
	}
	return res;
}

// return the sign(x)
vector<int> CCycleCheck::sign(vector<int> vec)
{
	vector<int> res;
	for (vector<int>::iterator data = vec.begin(); data != vec.end(); data++)
	{
		if (*data > 0)
		{
			res.push_back(1);
		}
		else if (*data < 0)
		{
			res.push_back(-1);
		}
		else
		{
			res.push_back(0);
		}
	}
	return res;
}

vector<int> CCycleCheck::sign(vector<float> vec)
{
	vector<int> res;
	for (vector<float>::iterator data = vec.begin(); data != vec.end(); data++)
	{
		if (*data > 0)
		{
			res.push_back(1);
		}
		else if (*data < 0)
		{
			res.push_back(-1);
		}
		else
		{
			res.push_back(0);
		}
	}
	return res;
}

// filter wave peaks based on pulse wave characteristics
// peaks_index: peak's index in original data vector
// peaks_data: the value of the peaks	 
WavePeakFeature* CCycleCheck::peaks_filter(WavePeakFeature &feature_parame)
{
	WavePeakFeature * peak_feature = new WavePeakFeature;
	//做一些过滤比较好，公司项目保密，省略不影响使用
	return peak_feature;
}

第一行两个周期序列，第二行两个杂波