之前在不經意間也有接觸過求突變點的問題。在我看來，與其說是求突變點，不如說是我們常常玩的"找不同"。給你兩幅影像，讓你找出兩個影像中不同的地方，我認為這其實也是找突變點在生活中的應用之一吧。回到找突變點位置上，以前自己有過一個傻傻的方法：就是直接求前後兩個取樣的的差分值，最後設定一個閾值，如果差分值大於這個閾值則該點是突變點。但這個方法問題很大，實際中突變點幅值有大有小，你怎麼能確定閾值到底是多少呢？還有可能訊號本來的差分值就比你那突變點的差分值還要大。所以這種方法在訊號或噪聲稍微複雜一點就行不通了。
這幾天看到了一種船新的訊號突變點檢測的方法-基於小波變換的訊號突變點檢測。於是乎去學習了一下小波變換的相關知識和應用，學習的不是很深入，但也模模糊糊感覺到了小波變換確實是檢測突變點的一大利器，下面分為兩個大部分總結一下我所學習到的小波變換求突變點的實現過程和相關知識理論。
小波變換求訊號突變點實現
我喜歡直接從應用入手，或者應用結合理論。一步一步分析程式碼，看資料和影像的變化比一步一步推公式有趣的多(雖然可能是錯誤的呀)。於是在這裡我就先直接上程式碼和影像了，這樣先讓我們對整個過程有個感性的認識。
原始訊號的生成
首先生成原始訊號，這裡隨便什麼訊號都可以，那我就生成一個正弦訊號吧，具體訊號引數見程式碼註釋。
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的函式用起來就是爽有木有。但是能應用是一回事，我們還是儘量多瞭解一下小波變換的原理為好。小波是怎麼構造的，它的性質有什麼？連續小波變換的影像是怎麼計算出來的呢?離散小波變換的每一層又是怎麼算出來的呢？只有學習了它們背後的支撐運算的數學公式，我們才能算真正理解了它。