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

以前在不經意間也有接觸過求突變點的問題。在我看來，與其說是求突變點，不如說是咱們經常玩的"找不一樣"。給你兩幅圖像，讓你找出兩個圖像中不一樣的地方，我認爲這其實也是找突變點在生活中的應用之一吧。回到找突變點位置上，之前本身有過一個傻傻的方法：就是直接求先後兩個採樣的的差分值，最後設置一個閾值，若是差分值大於這個閾值則該點是突變點。但這個方法問題很大，實際中突變點幅值有大有小，你怎麼能肯定閾值究竟是多少呢？還有可能信號原本的差分值就比你那突變點的差分值還要大。因此這種方法在信號或噪聲稍微複雜一點就行不通了。

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

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

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

原始信號的生成

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

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'也是一個小波的名字。處理完畢後獲得一、二、3層的細節係數(details)和第3層的近似係數(Approximations)。畫出這些係數的圖像以下：

由上圖可明顯看出，除去開頭和結尾的一些比較大的點外，在第一、二、3層的細節信號中，最大值點偏偏是第233點和第666點，因而也能夠愉快的能夠肯定這兩個點便是突變信號的位置了。

這裏還能夠稍微注意一下近似信號a3，它相似於原始信號濾去了高頻成分的樣子，它是怎麼得來的你看了第二部分就知道了！

總結

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

---

小波變換的基礎知識