小波變換檢測信號突變點的MATLAB實現

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

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

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

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

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幅度譜來做對比。從肉眼來看，我們可以發現原始信號和添加突變信號的頻域差別不大，只是突變信號的頻譜在高頻部分多了些抖動。所以要從傅里葉變換的頻域來檢測突變信號是不合適的。具體原因在第二部分會有總結，主要是兩個變換選取“基”的不同而導致的。

1.4 對信號做小波變換

重頭戲小波變換來了，這里我們用兩種小波變換的方法檢測突變點，一是連續小波變換；二是離散小波變換，這里只會簡略說明一下圖像，可以結合第二部分原理一起查看。

1.4.1 連續小波變換

我們對突變信號進行連續小波變換，實現代碼和圖像如下：

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較小時，可以看到一條較亮的白條，可以暫且理解成小波在這個位移和尺度上與信號相關性較大。在某個位置出現小波與信號相關性激增的原因就是信號在這個位置出現了突變，於是我們就愉快的找到了兩個突變點的位置。

1.4.2 離散小波變換

由於連續小波變換的位移參數(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，它類似於原始信號濾去了高頻成分的樣子，它是怎么得來的你看了第二部分就知道了！

1.5 總結

在這一部分中我們直抓要害，知道了怎么通過小波變換來檢測信號的突變點，MATLAB的函數用起來就是爽有木有。但是能應用是一回事，我們還是盡量多了解一下小波變換的原理為好。小波是怎么構造的，它的性質有什么？連續小波變換的圖像是怎么計算出來的呢?離散小波變換的每一層又是怎么算出來的呢？只有學習了它們背后的支撐運算的數學公式，我們才能算真正理解了它。

---

二：小波變換的基礎知識

關於基礎知識這一部分，主要都是參考的[2]里面的內容，我只是提煉了一些我認為重要的內容寫出來，然后有自己小小的補充。

2.1 連續小波變換(CWT)

首先直接給出連續小波變換公式：

\[C(a,b) = \frac{1}{\sqrt{a}} \int_R f(t) \psi^*(\frac{t-b}{a})dt \]

連續小波變換的結果就是得到很多個小波系數\(C(a,b)\)，它有兩個自變量參數，分別是尺度參數a和位移參數b。注意在CWT中a和b都是連續變化的。

2.1.1 尺度參數與位移參數

尺度變換：
如何對小波進行尺度變換呢？其實就是簡單的拉伸或者壓縮，下圖表現為一個小波在不同的尺度參數下的變換，可以看到尺度參數a越小，小波越壓縮，頻率越高；尺度參數a越大，小波越拉伸，頻率越低。

位移變換：
位移變換意味着將小波延遲或提前，如下圖將小波\(\psi(t)\)往右移位長度k得到\(\psi(t-k)\)。

2.1.2 連續小波變換具體過程

連續小波變換其實就是把小波從小尺度拉伸到大尺度，然后把不同尺度的小波依次從0位移到信號的完整長度，並不斷地計算它們的積分的過程。詳細來說就是下面幾個步驟：

1. 將小波\(\psi(t)\)放到原始信號\(f(t)\)的開頭處進行比較。
2. 計算小波系數C，C其實也表示了小波與信號這一部分的相關程度。C越大，說明相似度越高。

3. 將小波向右平移，距離為b，小波函數變為\(\psi(t-b)\)。並且重復步驟1和2，直到小波位移完整個信號\(f(t)\)。

4. 擴展小波的尺度，比如擴展一倍，小波函數變為\(\psi(\frac{t}{2})\)。然后重復步驟1~3。

5. 重復步驟1~4直到小波已經拓展到規定的最大尺度。

進行完這五個步驟之后，我們就能夠得到在不同的小波尺度和不同的信號位置上的小波系數了。如何更直觀的理解這個系數呢？於是我們可以把\(C(a,b)\)畫出來，x軸代表信號的時間(或說小波位移的長度b)；y軸表示尺度a；圖中每個點的顏色淺深表示C(a,b)的大小。下面是一個系數圖，其實在第一部分中也有我自己信號的CWT系數圖，忘了的可以回去看一下。

當然也可以看系數圖的側面視角，圖像如下：

回到文章主題，檢測信號突變點上。如果信號存在突變點的話，總有一些小尺度的小波在經過突變點這個位置的時候，此時計算出的\(C(a,b)\)會比周圍的點都大。從系數圖中我們可以想到這些位置會更亮一些，於是就找到了突變點。

2.2 離散小波變換(DWT)

\(f(t)\)的二進制連續小波變換為：

\[C(2^j,b) = \frac{1}{\sqrt{2^j}} \int_R f(t) \psi^*(\frac{t-b}{2^j})dt \tag{1} \]

當\(b=k2^j\)，式(1)為離散小波變換；當\(b=k\)，式(1)為平穩小波變換(二進小波變換)。
平穩小波變換只對尺度參數進行了離散化，而對時間域上的平移參數保持整數連續變化，平穩小波變換未破壞信號在時間域上的平移不變量。

DWT有多種實現方式，除了按照上述公式做的離散小波變換，其實還有一種更為常用的離散小波變換的方法，叫做Mallat算法，其實我們在第一部分將信號分成3層近似系數和細節系數就是用的這種方法。

2.2.1 半子帶濾波

對於大多數信號來說，信號的低頻部分是對信號整體特征的刻畫；而高頻部分則表現了信號的細枝末節。比如我們的聲音，如果把聲音的高頻部分去去掉，雖然聽起來有點不一樣，但我們仍然能夠聽得出說的什么內容，但如果把足夠多的低頻分量去掉，那就完全不知道在說什么了。

在小波分析中，我們把信號的低頻部分稱之為近似系數A(Approximations)，把信號的高頻部分稱作細節系數D(Details)，於是我們可以通過一個半子帶濾波器來得到A和D。

原始信號S，通過兩個互補的濾波器，生成兩個信號A和D。

值得注意的一點是，假設原始實數字信號S有1000個采樣點的長度，那么經過濾波后，A和D分別都會有1000個點，它們加在一起有2000個點，就比S還長了。由於之后重構S的需要，現在我們使用下采樣的方法將A和D各自降成500個點，方法是每兩個點取一個點就好了。下采樣之后的信號分別是cA1和cD1。

2.2.2 離散小波分解

在得到近似系數cA1之后，對cA1也進行半子帶濾波加下采樣，得到cA2和cD2，重復這個操作，可以最多進行\(log_2{N}\)層小波分解。

三：完整代碼

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
% x = x + 0.1 .* randn(1,1000); % 添加噪聲
x(233) = x(233) + 0.5;    % 添加突變點
x(666) = x(666) + 0.1;

figure(1); 
plot(x)
xlabel('采樣點數(1ms/點)');ylabel('幅值'); 
title('突變信號');       

%% 傅里葉變換觀察突變信號頻譜
Y = fft(x,1024); % 對信號進行傅里葉變換
f = Fs * (0 : (L / 2)) / L;
P2 = abs(Y / L);
P1 = P2(1 : L / 2 + 1);
figure(2)
plot(f,P1) 
title('突變信號的單邊幅度頻譜')
xlabel('f(Hz)')
ylabel('|P1(f)|')
axis([0,100,0,0.5])

%% 連續小波變換(CWT)
figure(3)
cw1 = cwt(x,1:32,'sym2','plot'); % 對信號做連續小波變換，並作出系數圖像
title('連續小波變換系數圖');


%% 離散小波變換(DWT) Wallat算法
%法一：直接用wavedec()進行3層分解，再重構生成近似系數和細節系數
[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層細節系數  

figure(4); 
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('時間'); 

%% 離散小波變換(DWT) Wallat算法
% 法二：用dwt()一層一層分解生成近似系數和細節系數
[ca11,cd1] = dwt(x,'db4');      % 第1層分解
[ca22,cd2] = dwt(ca11,'db4');   % 第2層分解
[ca3,cd3] = dwt(ca22,'db4');    % 第3層分解

figure(5)
subplot(511);plot(x);    % 畫出各層小波系數，注意由於函數自帶下采樣，所以每層系數長度會減半
title('原始信號x');            
subplot(512);plot(ca3);         
title('近似系數ca3');
subplot(513);plot(cd3);
title('細節系數cd3');
subplot(514);plot(cd2);
title('細節系數cd2');
subplot(515);plot(cd1);
title('細節系數cd1');
xlabel('時間'); 

% 為了與法一做對比，把系數又上采樣回1000點，和figure(4)圖像比較
cd_1 = dyadup(cd1);         % 上采樣到1000個點
cd_2 = dyadup(dyadup(cd2)); % 上采樣到1000個點
cd_3 = dyadup(dyadup(dyadup(cd3))); % 上采樣到1000個點
ca_3 = dyadup(dyadup(dyadup(ca3))); % 上采樣到1000個點

figure(6)
subplot(411);plot(ca_3);ylabel('近似信號a3');   %畫出各層小波系數
title('小波分解示意圖(方法二)');
axis([0 1000 -2 2]);
subplot(412);plot(cd_3);ylabel('細節信號d3');
axis([0 1000 -0.1 0.1]);
subplot(413);plot(cd_2);ylabel('細節信號d2');
axis([0 1000 -0.2 0.2]);
subplot(414);plot(cd_1);ylabel('細節信號d1');
axis([0 1000 -0.5 0.5]);
xlabel('時間'); 

參考文獻

[1] 張德豐.基於小波的信號突變點檢測算法研究[J].計算機工程與科學,2007(12):98-100.

[2] Misiti, Michel, et al. "Wavelet Toolbox 4–User’s Guide The MathWorks." Inc., Massachusetts, USA (2007).

[3] 吳茂. MATLAB R2016a 通信系統建模與仿真 28 個案例分析. 清華大學出版社, 2018.