之前在不经意间也有接触过求突变点的问题。在我看来,与其说是求突变点,不如说是我们常常玩的"找不同"。给你两幅图像,让你找出两个图像中不同的地方,我认为这其实也是找突变点在生活中的应用之一吧。回到找突变点位置上,以前自己有过一个傻傻的方法:就是直接求前后两个采样的的差分值,最后设置一个阈值,如果差分值大于这个阈值则该点是突变点。但这个方法问题很大,实际中突变点幅值有大有小,你怎么能确定阈值到底是多少呢?还有可能信号本来的差分值就比你那突变点的差分值还要大。所以这种方法在信号或噪声稍微复杂一点就行不通了。
这几天看到了一种船新的信号突变点检测的方法-基于小波变换的信号突变点检测。于是乎去学习了一下小波变换的相关知识和应用,学习的不是很深入,但也模模糊糊感觉到了小波变换确实是检测突变点的一大利器,下面分为两个大部分总结一下我所学习到的小波变换求突变点的实现过程和相关知识理论。
一:小波变换求信号突变点实现
我喜欢直接从应用入手,或者应用结合理论。一步一步分析代码,看数据和图像的变化比一步一步推公式有趣的多(虽然可能是错误的呀)。于是在这里我就先直接上代码和图像了,这样先让我们对整个过程有个感性的认识。使用小波变换进行的信号的分解和重构
1.1 原始信号的生成
首先生成原始信号,这里随便什么信号都可以,那我就生成一个正弦信号吧,具体信号参数见代码注释。
clear all; close all; clc;
Fs = 1000; % 采样频率1000Hz
Ts = 1 / Fs; % 采样时间间隔1ms
L = 1000; % 采样点数1000
t = (0 : L - 1) * Ts; % 采样时间。1000个点,每个点1ms,相当于采集了1s
x = sin(2 * pi * 10 * t); % 原始正弦信号,频率为10Hz,振幅为1
1.2 添加突变点
第二步我们要人为添加突变点了,为了看起来直观就暂时不添加噪声了。此处我们添加两个突变点,将第233个点的幅度在原本基础上增加0.5,将第666个点的幅度在原本基础上增加0.1,代码和添加后信号图像如下:
x(233) = x(233) + 0.5;
x(666) = x(666) + 0.1;
可以看到一个突变点很明显,而另一个却不是那么的明显,可能肉眼看的话都会忽略掉这个突变点。
1.3 对信号做傅里叶变换
可能有人要问,既然我们做的是小波变换,为什么又要对信号做傅里叶变换呢?其实我们确实可以不用做傅里叶变换的,但是为了与小波变换做对比,分析各自的优势和劣势,我们还是看一下该突变信号的傅里叶变换。
Y = fft(x,1024);
f = Fs * (0 : (L / 2)) / L;
P2 = abs(Y / L);
P1 = P2(1 : L / 2 + 1);
plot(f,P1)
title('突变信号的单边幅度频谱')
xlabel('f(Hz)')
ylabel('|P1(f)|')
axis([0,100,0,0.5])
1.3.1补充
下面我们再给一个原始信号的fft幅度谱来做对比。从肉眼来看,我们可以发现原始信号和添加突变信号的频域差别不大&#