离散信号的周期性判定，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 <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 frequencyint 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;
}

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

离散信号的周期性判定，C++实现相关推荐

【信号与系统】傅里叶变换的离散型与周期性
傅里叶变换的离散型与周期性文章目录傅里叶变换的离散型与周期性前言 1.连续时间与连续频率 2.连续时间与离散频率 3.离散时间与连续频率 4.离散时间与离散频率前言通过傅里叶级数,即周期函数 ...
离散信号（一） | 信号的采样和恢复+时域、频域采样定理
离线信号是指在时间上是离散的,即只在某些不连续的规定时刻给出信号的瞬时值,而在其它时刻无意义的信号.连续时间信号的采样是离散信号产生的方法之一,而计算机技术的发展以及数字技术的广泛应用是离散信号分析. ...
信号系统｜信号的分类｜确定信号与随机信号连续信号与离散信号周期信号与非周期信号能量信号与功率信号奇异信号
目录信号所有分类: 一.确定信号与随机信号二.连续信号与离散信号三.周期信号与非周期信号四.能量信号与功率信号五.奇异信号信号所有分类: 1.确定信号与随机信号 2.连续信号与离散信号 3 ...
学习笔记之——基于matlab的数字通信系统（2）之离散信号的傅里叶分析
关于连续信号的傅里叶分析,可以参考博文<学习笔记之--基于matlab的数字通信系统(1)&连续信号的傅里叶分析> 目录离散时间信号的傅里叶变换(DTFT) 连续时间信号的抽样- ...
模拟信号与离散信号之间的频率关系（由模拟信号采样得到的离散信号）
先看看模拟信号与离散信号之间的关系,这里的离散信号是指由模拟信号采样得到的离散信号(这是得到离散信号的方式之一),我想看看它们频率之间的关系? 以正弦信号为例: 从上面的手稿中可以看出,由模拟信号经过 ...
离散信号与系统分析(上)
离散信号与系统分析一.利用MATLAB产生离散信号 1.前言部分 stem(X,Y):在X的指定点处画出数据序列Y: stem(X,Y,'filled'):以实心的方式画出茎秆: axis([xmi ...
离散信号内插matlab,离散信号和系统实验报告.doc
离散信号和系统实验报告三.实验效果分析(包括仪器设备等使用效果) 实验中注意中英文的切换是应注意,特别是标点符号,括号等需要注意.需要仔细键入程序. 教师评语指导老师年月日江西师范大 ...
数字信号处理学习笔记[1] 离散信号奇异信号抽样定理
文章目录 2 离散信号和抽样定理 2.1 离散信号奇异信号 2.2 连续信号的离散化,正弦波的抽样问题 2.3 带限信号与奈奎斯特频率用卷积考察抽样定理 2.4 离散信号的频谱和抽样定理 2 离散 ...
【数字信号处理】Python离散信号卷积的代码实现/时域直接法/列表法/信号与系统
Python离散信号卷积的代码实现(时域直接法) 1.卷积卷积是一种积分变换的数学方法,在许多方面得到了广泛应用.卷积是两个变量在某范围内相乘后求和的结果.如果卷积的变量是离散的数组/数列,则卷积的 ...

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

离散信号的周期性判定，C++实现相关推荐

最新文章

热门文章