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

问题：给定一列离散信号，判定是否为周期信号

解决方案：

理论上，计算信号的自相关，如果是周期性信号，则其自相关序列依然为周期性信号切几乎不会衰减；否则，则会出现逐渐衰减至0。

实际情况下由于噪声的存在，偶尔自相关的最大值不出现在τtauτ=0处，而且如果τtauτ值取的比较大的时候即便是周期信号也会出现峰值衰减很大的情况，所以一般τtauτ值取N/2 N=len(vec)。具体理论请参见《数字信号处理》第四版P90的2.6.3小节。

算法的核心代码是自相关的计算，在实际编码的时候对书上的公式做了一个小小的处理，将ryy(l)=1M∑n=0M?1[x(n)x(n?l)]{r}_{yy}(l)=frac{1}{M}sum_{n=0}^{M-1}[x(n)x(n-l)]ryy?(l)=M1?∑n=0M?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

#include

using namespace std;

typedef struct WavePeakFeature

{

vector peaks_index;

vector peaks_value;

vector freq_list;

}WavePeakFeature;

class CCycleCheck

{

public:

CCycleCheck();

CCycleCheck(int);

~CCycleCheck();

bool CycleCheck(vectorinput_data);

// sampling frequency

int sample_freq;

public:

static const int default = 0;

static const int unbiased = 1;

private:

// computer vec's autocorrelation

vector xcorr(vectorvec,int model);

// computer the difference for the vec[i]-vec[i-1]

vector diff(vectorvec);

vector diff(vectorvec);

// return the sign(x)

vector sign(vectorvec);

vector sign(vectorvec);

vector median_filter(vector vec,int window);

WavePeakFeature find_peaks(vector vec);

WavePeakFeature find_peaks(vector 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

#include

#include

#include

CCycleCheck::CCycleCheck()

{

}

CCycleCheck::CCycleCheck(int sample_freq)

{

this->sample_freq = sample_freq;

}

CCycleCheck::~CCycleCheck()

{

}

bool CCycleCheck::CycleCheck(vector input_data)

{

bool res = false;

try

{

vector data = median_filter(input_data,3);

vector 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 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 vec)

{

WavePeakFeature peak_feature;

vector vec_sign = sign(diff(vec));

vector 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 vec)

{

WavePeakFeature peak_feature;

vector vec_sign = sign(diff(vec));

vector 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 CCycleCheck::median_filter(vector vec,int kernel_size)

{

vector median;

vector::iterator iter = vec.begin();

for (int i = 0; i < vec.size() - kernel_size;i++)

{

vector 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 CCycleCheck::xcorr(vector vec,int model)

{

int vec_len = vec.size();

vector res;

vector zero_vec(vec_len / 2,0);

vector 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 CCycleCheck::diff(vectorvec)

{

vector res;

for (int i = 1; i < vec.size(); i++)

{

res.push_back(vec[i] - vec[i - 1]);

}

return res;

}

vector CCycleCheck::diff(vectorvec)

{

vector res;

for (int i = 1; i < vec.size(); i++)

{

res.push_back(vec[i] - vec[i - 1]);

}

return res;

}

// return the sign(x)

vector CCycleCheck::sign(vector vec)

{

vector res;

for (vector::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 CCycleCheck::sign(vector vec)

{

vector res;

for (vector::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;

}

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

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!