之前在不經意間也有接觸過求突變點的問題。在我看來，與其說是求突變點，不如說是我們常常玩的"找不同"。給你兩幅圖像，讓你找出兩個圖像中不同的地方，我認為這其實也是找突變點在生活中的應用之一吧。回到找突變點位置上，以前自己有過一個傻傻的方法：就是直接求前后兩個采樣的的差分值，最后設置一個閾值，如果差分值大于這個閾值則該點是突變點。但這個方法問題很大，實際中突變點幅值有大有小，你怎么能確定閾值到底是多少呢？還有可能信號本來的差分值就比你那突變點的差分值還要大。所以這種方法在信號或噪聲稍微復雜一點就行不通了。

這幾天看到了一種船新的信號突變點檢測的方法-基于小波變換的信號突變點檢測。于是乎去學習了一下小波變換的相關知識和應用，學習的不是很深入，但也模模糊糊感覺到了小波變換確實是檢測突變點的一大利器，下面分為兩個大部分總結一下我所學習到的小波變換求突變點的實現過程和相關知識理論。

小波變換求信號突變點實現

我喜歡直接從應用入手，或者應用結合理論。一步一步分析代碼，看數據和圖像的變化比一步一步推公式有趣的多(雖然可能是錯誤的呀)。于是在這里我就先直接上代碼和圖像了，這樣先讓我們對整個過程有個感性的認識。

原始信號的生成

首先生成原始信號，這里隨便什么信號都可以，那我就生成一個正弦信號吧，具體信號參數見代碼注釋。

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】

Roberts邊緣檢測算子:根據一對互相垂直方向上的差分可用來計算梯度的原理,采用對角線方向相鄰兩像素之差. 小波變換的方法比較適用于展現夾帶在正常信號中的瞬間反常現象,具有方向敏感性.所以可以邊緣檢 ...

基于小波變換的數字圖像處理(MATLAB源代碼)

基于小波變換的數字圖像處理(MATLAB源代碼) clear all; close all; clc;M=256;%原圖像長度N=64; %水印長度[filename1,pathname]=uiget ...

CM 安裝CDH 錯誤： 安裝失敗。 無法接收 Agent 發出的檢測信號。

在安裝CDH的時候出現錯誤提示: 安裝失敗. 無法接收 Agent 發出的檢測信號. 日志提示錯誤: start >> raise socket.error(msg) >>er ...

對AM信號FFT的matlab仿真

普通調幅波AM的頻譜,大信號包絡檢波頻譜分析 u(t)=Ucm(1+macos ?t)cos ?ct ma稱為調幅系數 它的頻譜由載波,上下邊頻組成 , 包絡檢波中二極管截去負半周再用電容低通濾波,可 ...

小波變換在matlab中的使用

對信號進行一層分解 clc; clear; % 獲取噪聲信號 load('matlab.mat'); sig = M(1,1:1400); SignalLength = length(sig); %使 ...

matlab工具箱之人眼檢測+meanshift跟蹤算法--人眼跟蹤

Viola-Jones 人眼檢測算法+meanshift跟蹤算法 這次的代碼是對視頻中的人眼部分進行檢測加跟蹤,檢測用的是matlab自帶的人眼檢測工具箱 下面是matlab官網介紹這個算法的一些東西 ...

paper 69：Haar-like矩形遍歷檢測窗口演示Matlab源代碼[轉載]

Haar-like矩形遍歷檢測窗口演示Matlab源代碼 clc; clear; close all; % Haar-like特征矩形計算 board = 24 % 檢測窗口寬度 num = 24 % ...

(轉)Haar-like矩形遍歷檢測窗口演示Matlab源代碼

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

matlab 霍夫變換—檢測圓

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

隨機推薦

[JS進階] 編寫可維護性代碼 (1)

今天的web應用大至成千上萬行的javascript代碼,執行各種復雜的過程,這種演化讓我們開發者必須得對可維護性有一定的認識!編寫可維護性代碼很重要,很多情況下我們是以他人的工作成果為基礎,確保代碼 ...

給jdk寫注釋系列之jdk1.6容器(5)-LinkedHashMap源碼解析

前面分析了HashMap的實現,我們知道其底層數據存儲是一個hash表(數組+單向鏈表).接下來我們看一下另一個LinkedHashMap,它是HashMap的一個子類,他在HashMap的基礎上維持 ...

使用BufferedReader的時候出現的問題

今天在使用BufferedReader的時候,出現了一個奇怪的問題 有時候換行的時候,行首會少一個字符 開始的代碼是這樣寫的 while( br.read()!=-1 ){ String str = ...

UVa 357 - Let Me Count The Ways

題目大意:也是硬幣兌換問題,與147.674用同樣的方法即可解決. #include #include #define MAXN 3000 ...

翻譯連載 | 附錄 C：函數式編程函數庫－《JavaScript輕量級函數式編程》 |《你不知道的JS》姊妹篇

總結

以上是生活随笔為你收集整理的matlab怎么分析突变点,小波变换检测信号突变点的MATLAB实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。