峰值检测概述和思路
最近集训任务点是做去年的国赛信号题，即信号失真度测量。也就是输入一个周期信号，经过AD转换之后做FFT变换，因为变换之后也还是复数所以还要求一下幅值，得出幅值数组之后再找出各次谐波幅值，归一化，最后计算失真度。
乍一听很简单对吧？可是从冗长的幅值数组里面找出来峰值的确还要花一点心思的。
在网上找到的多是一些算法程序，不太实用，于是就自己写了。
今天听到学长的思路是先找直流分量（最高）和基波（第二高）的频率，这俩相减之后再去找恰巧在倍数频率上面的幅值，相反的是，我手搓了一个函数专门用来先找谐波峰值+归一化+计算失真度的（还真是信号调理不够算法来补），下面来一点点看具体程序捏。（完整程序会贴在最后面，先分段看捏）
差分！一阶段！
我们都知道，要分析信号峰值在哪里，就要看信号的变化规律。在连续函数里面，我们通常把一阶导数的零点叫做“极值点”，那么要分析离散信号的变化规律，就要用到离散数据的求导——差分了。
差分，一般在大数据里用在以时间为统计维度的分析中，其实就是下一个数值 ，减去上一个数值。
当间距相等时，用下一个数值，减去上一个数值 ，就叫“一阶差分”，做两次相同的动作，即再在一阶差分的基础上用后一个数值再减上一个数值一次，就叫“二阶差分"。
理论成立，实践开始！
首先要说明整个程序其实是一个函数，参数是float *mag_input（经过fft变换和复数幅度计算之后的幅度数组，大小是NUM也就是FFT点数）。
定义一个数组sign用来存差分并赋值后的内容，在循环里对于整个输入数组差分，并根据差值给sign赋值。
为什么要赋值呢？因为在实际信号中，总会有小的波动，比如这样一组数字[1,2,1,2,4,3,108,62,7,2,3,1]，根据我们的感性认识，可以知道108才是我们需要的那个峰值。但你能说2，4，3这三个数不是小小的峰值吗？程序需要判断什么是我们要的峰值，什么是不要的杂乱波动。
如果在上述数组差分过程中，采用0作为判断上升与下降的标准，那么必然会出现大大小小的峰值干扰计算结果，那如果我们采用10呢？也就是说，对于一个离散的数据来说，它周围正负10以内的波动都算作没有波动，那么很轻易就能够得出108才是那个唯一有用的峰值。
下面程序段里面采用5000是因为我今天测试所用的电压偏移量比较大，幅值差自然也很大，是根据真实数据改出来的。如果寻求一种更加自动化的方法找这个边界值的话，这边建议求一下频谱幅值的最大值，取0.05-0.1就差不多了（我猜测啊没试过）。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
| int *sign = (int *)malloc(NUM / 2);
for (int i = 1; i < NUM / 2; i++)
{
float diff = mag_input[i] - mag_input[i - 1];
if (diff > 5000)
sign[i - 1] = 1;
else if (diff < -5000)
sign[i - 1] = -1;
else
sign[i - 1] = 0;
}
|
虽然差分！但是二阶段！
如果用10作为边界值的话，上面那个数组求得的sign序列是（放一起方便看）
1
2
| 1 2 1 2 4 3 108 62 72 3 1
0 0 0 0 0 1 -1 -1 -1 0
|
再差分一下可以发现变成了：
1
2
3
| 1 2 1 2 4 3 108 62 72 3 1
0 0 0 0 0 1 -1 -1 -1 0
0 0 0 0 1 -2 0 0 1
|
第二次差分之后是负数的就是我们的极大值点了！
下面程序还有一点要注意的是，假如我的序列是[1,2,300,301,2,1]
那么我差分就会：
1
2
3
| 1 2 300 301 2 1
0 1 0 -1 0
1 -1 -1 1
|
300和301对应的第二次差分值都是负数啊！实际很可能会出现这种情况（已经遇到过了…），所以要在两个峰值相邻的时候决定选出哪一个！
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
|
int indMax[6];
int index = 0;
for (int j = 1; j < NUM / 2 && index < 6; j++)
{
int diff = sign[j] - sign[j - 1];
if (diff < 0)
{
if (index > 0 && indMax[index - 1] == j - 1)
if (mag_input[j] - mag_input[indMax[index - 1]] > 0)
indMax[index - 1] = j;
else
j++;
else
{
indMax[index] = j;
index++;
}
}
}
free(sign);
|
然后我们就得到了一个存着峰值下标的数组indMax！
重新审视下标数组
1
2
3
4
5
6
7
8
9
10
11
12
| if (indMax[0] < 3)
{
for (int i = 0; i < index; i++)
indMax[i] = indMax[i + 1];
indMax[5] = 0;
index = index - 1;
}
|
峰值好汉们入座！都可以归一！
这里我的思路是：既然已经存好了哪些是峰值，那就一个个除以基波的序号（因为基波的频率横坐标=数组下标序号*频率分辨率），得到每个峰值到底对应的是几次谐波，高于五次的就不要了。顺便归一化一下，也不是什么难事儿。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
| float *mag_result = (float *)malloc(6);
for (int i = 0; i < 6; i++)
mag_result[i] = 0;
mag_result[0] = 1;
for (int n = 1, order; n < index; n++)
{
order = roundf((float)indMax[n] / (float)indMax[0]);
if (order <= 5)
mag_result[order - 1] = mag_input[indMax[n]] / mag_input[indMax[0]];
else
break;
}
|
一起来算算失真度（THD）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
| float squart_sum = 0;
float THD = 0;
printf("归一化幅值：\r\n");
for(int i = 0;i<5;i++)
{
if(i>0)
squart_sum += mag_result[i] * mag_result[i];
printf("%.2f ",mag_result[i]);
}
THD = sqrt(squart_sum)/mag_result[0];
printf("\r\nTHD = %.2f%%\r\n",THD*100);
|
完整程序
有很多不足，发现问题的话欢迎指导我捏。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
| #include "main.h"
#include <stdlib.h>
#define NUM = 512;
void findPeaks(float *mag_input)
{
int *sign = (int *)malloc(NUM / 2);
for (int i = 1; i < NUM / 2; i++)
{
float diff = mag_input[i] - mag_input[i - 1];
if (diff > 5000)
{
sign[i - 1] = 1;
}
else if (diff < -5000)
{
sign[i - 1] = -1;
}
else
{
sign[i - 1] = 0;
}
}
int indMax[6];
int index = 0;
for (int j = 1; j < NUM / 2 && index < 6; j++)
{
int diff = sign[j] - sign[j - 1];
if (diff < 0)
{
if (index > 0 && indMax[index - 1] == j - 1)
if (mag_input[j] - mag_input[indMax[index - 1]] > 0)
indMax[index - 1] = j;
else
j++;
else
{
indMax[index] = j;
index++;
}
}
}
free(sign);
float *mag_result = (float *)malloc(6);
for (int i = 0; i < 6; i++)
mag_result[i] = 0;
if (indMax[0] < 3)
{
for (int i = 0; i < index; i++)
indMax[i] = indMax[i + 1];
indMax[5] = 0;
index = index - 1;
}
mag_result[0] = 1;
for (int n = 1, order; n < index; n++)
{
order = roundf((float)indMax[n] / (float)indMax[0]);
if (order <= 5)
mag_result[order - 1] = mag_input[indMax[n]] / mag_input[indMax[0]];
else
break;
}
printf("归一化幅值：\r\n");
for(int i = 0;i<5;i++)
{
if(i>0)
squart_sum += mag_result[i] * mag_result[i];
printf("%.2f ",mag_result[i]);
}
THD = sqrt(squart_sum)/mag_result[0];
printf("\r\nTHD = %.2f%%\r\n",THD*100);
free(mag_result);
}
|
