之前在不经意间也有接触过求突变点的问题。在我看来,与其说是求突变点,不如说是我们常常玩的"找不同"。给你两幅图像,让你找出两个图像中不同的地方,我认为这其实也是找突变点在生活中的应用之一吧。回到找突变点位置上,以前自己有过一个傻傻的方法:就是直接求前后两个采样的的差分值,最后设置一个阈值,如果差分值大于这个阈值则该点是突变点。但这个方法问题很大,实际中突变点幅值有大有小,你怎么能确定阈值到底是多少呢?还有可能信号本来的差分值就比你那突变点的差分值还要大。所以这种方法在信号或噪声稍微复杂一点就行不通了。

这几天看到了一种船新的信号突变点检测的方法-基于小波变换的信号突变点检测。于是乎去学习了一下小波变换的相关知识和应用,学习的不是很深入,但也模模糊糊感觉到了小波变换确实是检测突变点的一大利器,下面分为两个大部分总结一下我所学习到的小波变换求突变点的实现过程和相关知识理论。

小波变换求信号突变点实现

我喜欢直接从应用入手,或者应用结合理论。一步一步分析代码,看数据和图像的变化比一步一步推公式有趣的多(虽然可能是错误的呀)。于是在这里我就先直接上代码和图像了,这样先让我们对整个过程有个感性的认识。

原始信号的生成

首先生成原始信号,这里随便什么信号都可以,那我就生成一个正弦信号吧,具体信号参数见代码注释。

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

添加突变点

第二步我们要人为添加突变点了,为了看起来直观就暂时不添加噪声了。此处我们添加两个突变点,将第233个点的幅度在原本基础上增加0.5,将第666个点的幅度在原本基础上增加0.1,代码和添加后信号图像如下:

x(233) = x(233) + 0.5;
x(666) = x(666) + 0.1;

可以看到一个突变点很明显,而另一个却不是那么的明显,可能肉眼看的话都会忽略掉这个突变点。

对信号做傅里叶变换

可能有人要问,既然我们做的是小波变换,为什么又要对信号做傅里叶变换呢?其实我们确实可以不用做傅里叶变换的,但是为了与小波变换做对比,分析各自的优势和劣势,我们还是看一下该突变信号的傅里叶变换。

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])

补充

下面我们再给一个原始信号的fft幅度谱来做对比。从肉眼来看,我们可以发现原始信号和添加突变信号的频域差别不大,只是突变信号的频谱在高频部分多了些抖动。所以要从傅里叶变换的频域来检测突变信号是不合适的。具体原因在第二部分会有总结,主要是两个变换选取“基”的不同而导致的。

对信号做小波变换

重头戏小波变换来了,这里我们用两种小波变换的方法检测突变点,一是连续小波变换;二是离散小波变换,这里只会简略说明一下图像,可以结合第二部分原理一起查看。

连续小波变换

我们对突变信号进行连续小波变换,实现代码和图像如下:

cw1 = cwt(x,1:32,'sym2','plot'); % 对信号做连续小波变换
title('连续小波变换');

cwt(Continuous wavelet transform)函数表示进行连续小波变换,主要是处理一维的数据,比如我们这个数据。参数x是输入的信号;1:32表示尺度参数Scales的取值范围为(1:32);'sym2'表示我们用的小波是sym2小波;'plot'是画出连续小波变换系数的意思。运行图像如下:

不同于傅里叶变换只有w一个自变量,小波变换有两个自变量,分别是a(尺度参数)和b(位移参数)。从上图我们可以看出在小波位移到第233个点和第666个点,且a较小时,可以看到一条较亮的白条,可以暂且理解成小波在这个位移和尺度上与信号相关性较大。在某个位置出现小波与信号相关性激增的原因就是信号在这个位置出现了突变,于是我们就愉快的找到了两个突变点的位置。

离散小波变换

由于连续小波变换的位移参数(b)和尺度参数(a)都是连续变化的,特别是尺度参数的连续变化,会带来巨大的计算量,于是利用离散小波变换来处理信号,这里还是主要说代码和图像,具体实现原理在第二部分有粗浅介绍。

[d,a]=wavedec(x,3,'db4');           %对原始信号进行3层离散小波分解
a3=wrcoef('a',d,a,'db4',3); %获得第3层近似系数
d3=wrcoef('d',d,a,'db4',3); %获得第3层细节系数
d2=wrcoef('d',d,a,'db4',2); %获得第2层细节系数
d1=wrcoef('d',d,a,'db4',1); %获得第1层细节系数 subplot(411);plot(a3);ylabel('近似信号a3'); %画出各层小波系数
title('小波分解示意图');
subplot(412);plot(d3);ylabel('细节信号d3');
subplot(413);plot(d2);ylabel('细节信号d2');
subplot(414);plot(d1);ylabel('细节信号d1');
xlabel('时间');

wavedec(wavelet decomposition)函数表示进行离散多辨小波分解,x是待处理的输入信号;3表示进行3层分解,'db4'也是一个小波的名字。处理完毕后得到1、2、3层的细节系数(details)和第3层的近似系数(Approximations)。画出这些系数的图像如下:

由上图可明显看出,除去开头和结尾的一些比较大的点外,在第1、2、3层的细节信号中,最大值点恰恰是第233点和第666点,于是也可以愉快的可以确定这两个点即是突变信号的位置了。

这里还可以稍微注意一下近似信号a3,它类似于原始信号滤去了高频成分的样子,它是怎么得来的你看了第二部分就知道了!

总结

在这一部分中我们直抓要害,知道了怎么通过小波变换来检测信号的突变点,MATLAB的函数用起来就是爽有木有。但是能应用是一回事,我们还是尽量多了解一下小波变换的原理为好。小波是怎么构造的,它的性质有什么?连续小波变换的图像是怎么计算出来的呢?离散小波变换的每一层又是怎么算出来的呢?只有学习了它们背后的支撑运算的数学公式,我们才能算真正理解了它。


小波变换的基础知识

小波变换检测信号突变点的MATLAB实现的更多相关文章

  1. 边缘检测算子和小波变换提取图像边缘【matlab】

    Roberts边缘检测算子:根据一对互相垂直方向上的差分可用来计算梯度的原理,采用对角线方向相邻两像素之差. 小波变换的方法比较适用于展现夹带在正常信号中的瞬间反常现象,具有方向敏感性.所以可以边缘检 ...

  2. 基于小波变换的数字图像处理(MATLAB源代码)

    基于小波变换的数字图像处理(MATLAB源代码) clear all; close all; clc;M=256;%原图像长度N=64; %水印长度[filename1,pathname]=uiget ...

  3. CM 安装CDH 错误: 安装失败。 无法接收 Agent 发出的检测信号。

    在安装CDH的时候出现错误提示: 安装失败. 无法接收 Agent 发出的检测信号. 日志提示错误: start >> raise socket.error(msg) >>er ...

  4. 对AM信号FFT的matlab仿真

    普通调幅波AM的频谱,大信号包络检波频谱分析 u(t)=Ucm(1+macos t)cos ct ma称为调幅系数 它的频谱由载波,上下边频组成 , 包络检波中二极管截去负半周再用电容低通滤波,可 ...

  5. matlab工具箱之人眼检测+meanshift跟踪算法--人眼跟踪

    Viola-Jones 人眼检测算法+meanshift跟踪算法 这次的代码是对视频中的人眼部分进行检测加跟踪,检测用的是matlab自带的人眼检测工具箱 下面是matlab官网介绍这个算法的一些东西 ...

  6. paper 69:Haar-like矩形遍历检测窗口演示Matlab源代码[转载]

    Haar-like矩形遍历检测窗口演示Matlab源代码 clc; clear; close all; % Haar-like特征矩形计算 board = 24 % 检测窗口宽度 num = 24 % ...

  7. (转)Haar-like矩形遍历检测窗口演示Matlab源代码

    from:http://blog.sina.com.cn/s/blog_736aa0540101kzqb.html clc; clear; close all; % Haar-like特征矩形计算 b ...

  8. matlab 霍夫变换—检测圆

    function [hough_space,hough_circle,para] = hough_Circle(BW,step_r,step_angle,r_min,r_max,p) % %%%%%% ...

  9. Matlab小波工具箱的使用2

    Matlab小波工具箱的使用2 (2011-11-11 09:32:57) 转载▼ http://blog.sina.com.cn/s/blog_6163bdeb0102dw7a.html#cmt_5 ...

  10. 基于FPGA的音频信号的FIR滤波(Matlab+Modelsim验证)

    1 设计内容 本设计是基于FPGA的音频信号FIR低通滤波,根据要求,采用Matlab对WAV音频文件进行读取和添加噪声信号.FFT分析.FIR滤波处理,并分析滤波的效果.通过Matlab的分析验证滤 ...

随机推荐

  1. ROS系统C++代码测试之gtest

    1. 安装gtestsudo apt-get install libgtest-dev 2.修改CMakeLists.txtfind_package(GTest REQUIRED)uncommend ...

  2. G - 小希的迷宫

    Time Limit:1000MS     Memory Limit:32768KB     64bit IO Format:%I64d & %I64u Submit Status Pract ...

  3. 设置Ubuntu 14.04右键终端的方法

    设置Ubuntu 14.04右键终端的方法如下: 首先要安装一个包,即可在右键里面添加一个"打开终端"的菜单. sudo apt-get install nautilus-open ...

  4. Reveal-Plugin-for-Xcode 自动结合 Reveal 进行 UI 分析

    下载地址:https://github.com/shjborage/Reveal-Plugin-for-Xcode 还记得之前我们如何使用 Reveal UI 分析工具进行实时查看 UI 的结构吗?如 ...

  5. UIImagePickerController显示中文界面

    iOS开发中,我们经常遇到获取拍照.相册中图片的功能,就必然少不了UIImagePickerController,但是我们发现当我们使用它的时候,它的页面是英文的,看着很别扭,国人还是比较喜欢看中文界 ...

  6. hdu 4274 Spy's Work(水题)

    Spy's Work Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others) Total ...

  7. LeetCode OJ 117. Populating Next Right Pointers in Each Node II

    题目 Given a binary tree struct TreeLinkNode { TreeLinkNode *left; TreeLinkNode *right; TreeLinkNode * ...

  8. 安卓开发 1配置环境 (JDK+Andriod stiuio)

    配置JDK+Andriod stiuio 所需要用到的JDK和AS 链接:https://pan.baidu.com/s/1smHCD1z密码:rxh2 1 JAVA_SDK   1.1 下载 htt ...

  9. JavaScript简介及作用

    JavaScript是一门脚本语言,是可以插入HTML页面的编程代码,插入HTML以后可以由所有现代浏览器运行 一.写如html输出 <body> <script> docum ...

  10. Remote error: VAR and OUT arguments must match parameter type exactly&#39;

    在XE10中 downloadfile(filename: string; out FileStream: TStream; out FileSize: int64)是没有问题的,升级到delphi ...