小波分析的理解


  小波變換是克服其他信號處理技術缺陷的一種分析信號的方法。

  小波由一族小波基函數構成,它可以描述信號時間(空間)和頻率(尺度)域的局部特性。采用小波分析最大優點是可對信號進行實施局部分析,可在任意的時間或空間域中分析信號。小波分析具有發現其他信號分析方法所不能識別的、隱藏於數據之中的表現結構特性的信息,而這些特性對機械故障和材料的損傷等識別是尤為重要的。如何選擇小波基函數目前還沒有一個理論標准,常用的小波函數有 Haar、 Daubechies(dbN)、 Morlet、 MeryerSymletCoifletBiorthogonal 小波等15種。

  但是小波變換的小波系數為如何選擇小波基函數提供了依據。小波變換后的系數比較大,就表明了小波和信號的波形相似程度較大;反之則比較小。  另外還要根據信號處理的目的來決定尺度的大小。如果小波變換僅僅反映信號整體的近似特征,往往選用較大的尺度;反映信號細節的變換則選用尺度不大的小波。由於小波函數家族成員較多,進行小波變換目的各異,目前沒有一個通用的標准。
    根據實際運用的經驗,Morlet小波應用領域較廣,可以用於信號表示和分類、圖像識別特征提取;墨西哥草帽小波用於系統識別;樣條小波用於材料探傷;Shannon正交基用於差
分方程求解。

 現在對小波分解層數與尺度的關系作如下解釋:


      是不是小波以一個尺度分解一次就是小波進行一層的分解?
      比如:[C,L]=wavedec(X,N,'wname')中,N為尺度,若為1,就是進行單尺度分解,也就是分解一層。    但是W=CWT(X,[2:2:128],'wname','plot')的分解尺度又是從21282為步進的,這里的“分解尺度”跟上面那個“尺度”的意思一樣嗎?
      [C,L]=wavedec(X,N,'wname')中的N為分解層數不是尺度,'wname'DB小波為例DB4, 4為消失矩,則一般濾波器長度為8, 階數為7.
     wavedec針對於離散,CWT是連續的。
       多尺度又是怎么理解的呢?
      多尺度的理解:  如將0-pi定義為空間V0,  經過一級分解之后V0被分成0-pi/2的低頻子空間V1pi/2-pi的高頻子空間W1,  然后一直分下去....得到 VJ+WJ+....W2+W1.    因為VJWJ是正交的空間且各W子空間也是相互正交的.  所以分解得到了是相互不包含的多個頻域區間,這就是多分辯率分析即多尺度分析.
      當然多分辨率分析是有嚴格數學定義的,但完全可以從數字濾波器角度理解它.當然,你的泛函學的不錯,也可以從函數空間角度理解.


       是不是說分解到W3W2W1V3就是三尺度分解?
       簡單的說尺度就是頻率,不過是反比的關系.確定尺度關鍵還要考慮你要分析信號的采樣頻率大小,因為根據采樣頻率大小才能確定你的分析頻率是多少.(采樣定理).然后再確定你到底分多少層.


       假如我這有一個10hz50hz的正弦混合信號,采樣頻率是500hz,是不是就可以推斷出10hz50hz各自對應的尺度了呢?我的意思是,是不是有一個頻率和尺度的換算公式?
       實際頻率=小波中心頻率×采樣頻率/尺度


         在小波分解中,若將信號中的最高頻率成分看作是1,則各層小波小波分解便是帶通或低通濾波器,且各層所占的具體頻帶為(三層分解)a1:0~0.5 d1: 0.5~1; a2:0~0.25 d2: 0.25~0.5; a3: 0~0.125; d3:0.125~0.25    

 

  可以這樣理解嗎?如果我要得到頻率為0.125~0.25的信號信息,是不是直接對d3的分解系數直接重構之后就是時域信息了?這樣感覺把多層分解純粹當作濾波器來用了,又怎么是多分辨分析??怎樣把時頻信息同時表達出來??
       這個問題非常好,我剛開始的時候也是被這個問題困惑住了,咱們確實是把它當成了濾波器來用了,也就是說我們只看重了小波分析的頻域局部化的特性。但是很多人都忽略其時域局部化特性,因為小波是變時頻分析的方法,根據測不准原理如果帶寬大,則時窗寬度就要小。那么也就意味着如果我們要利用其時域局部化特性就得在時寬小的分解層數下研究,也就是低尺度下。這樣我們就可以更容易看出信號在該段時間內的細微變化,但是就產生一個問題,這一段的頻率帶很寬,頻率局部化就體現不出來了。
       d3進行單支重構就可以得到0.1250.25的信號了,當然頻域信息可能保存的比較好,
但如果小波基不是對稱的話,其相位信息會失真。


      小波變換主要也是用在高頻特征提取上。

      層數不是尺度,小波包分解中,N應該是層數,個人理解對應尺度應該是2^N

     小波分解的尺度為a,分解層次為j   如果是連續小波分解尺度即為a。離散小波分解尺度嚴格意義上來說為a2^j,在很多書上就直接將j稱為尺度,因為一個j就對應者一個尺度a。其實兩者是統一的。


小波基:一般從線性相位,消失矩,相似性,緊支撐等來選擇。

Daubechies小波基的構造
%  此程序實現構造小波基
%  periodic_wavelet.m
function ss=periodic_wavelet;
clear;clc;
% global MOMENT;  %  消失矩階數
% global LEFT_SCALET;  %  尺度函數左支撐區間
% global RIGHT_SCALET;  %  尺度函數右支撐區間
% global LEFT_BASIS;  %  小波基函數左支撐區間
% global RIGHT_BASIS;  %  小波基函數右支撐區間
% global MIN_STEP;  %  最小離散步長
% global LEVEL;  %  計算需要的層數(離散精度)
% global MAX_LEVEL;  %  周期小波最大計算層數

[s2,h]=scale_integer;
[test,h]=scalet_stretch(s2,h);
wave_base=wavelet(test,h);
ss=periodic_waveletbasis(wave_base);
 

function [s2,h]=scale_integer;
%  本函數實現求解小波尺度函數離散整數點的值
%  sacle_integer.m
MOMENT=10;  %  消失矩階數
LEFT_SCALET=0;  %  尺度函數左支撐區間
RIGHT_SCALET=2*MOMENT-1;  %  尺度函數右支撐區間
LEFT_BASIS=1-MOMENT;    %  小波基函數左支撐區間
RIGHT_BASIS=MOMENT;     %  小波基函數右支撐區間
MIN_STEP=1/512;          %  最小離散步長
LEVEL=-log2(MIN_STEP);  %  計算需要的層數(離散精度)
MAX_LEVEL=8;  %  周期小波最大計算層數

h=wfilters('db10','r');  %  濾波器系數
h=h*sqrt(2); % FI(T)=SQRT(2)*SUM(H(N)*FI(2T-N)) N=0:2*MOMENT-1;
for i=LEFT_SCALET+1:RIGHT_SCALET-1
    for j=LEFT_SCALET+1:RIGHT_SCALET-1
       k=2*i-j+1;
       if (k>=1&k<=RIGHT_SCALET+1)
       a(i,j)=h(k);  %  矩陣系數矩陣
       else
       a(i,j)=0;
       end
    end
end
[s,w]=eig(a);  %  求特征向量,解的基
s1=s(:,1);
s2=[0;s1/sum(s1);0]; %  根據條件SUM(FI(T))=1,求解;
 
 
%  本函數實現尺度函數經伸縮后的離散值
%  scalet_stretch.m
function [s2,h]=scalet_stretch(s2,h);
MOMENT=10;  %  消失矩階數
LEFT_SCALET=0;  %  尺度函數左支撐區間
RIGHT_SCALET=2*MOMENT-1;  %  尺度函數右支撐區間
LEFT_BASIS=1-MOMENT;  %  小波基函數左支撐區間
RIGHT_BASIS=MOMENT;  %  小波基函數右支撐區間
MIN_STEP=1/512;  %  最小離散步長
LEVEL=-log2(MIN_STEP);  %  計算需要的層數(離散精度)
MAX_LEVEL=8;  %  周期小波最大計算層數

for j=1:LEVEL  %  需要計算到尺度函數的層數
   t=0;
   for i=1:2:2*length(s2)-3  %  需要計算的離散點取值(0123 -> 1/2, 3/2, 5/2)
      t=t+1;
      fi(t)=0;
      for n=LEFT_SCALET:RIGHT_SCALET;  % 低通濾波器沖擊響應緊支撐判斷
          if ((i/2^(j-1)-n)>=LEFT_SCALET&(i/2^(j-1)-n)<=RIGHT_SCALET) %  小波尺度函數緊支撐判斷
            fi(t)=fi(t)+h(n+1)*s2(i-n*2^(j-1)+1);  %  反復應用雙尺度方程求解
          end
      end
   end
   clear s
   n1=length(s2);
   n2=length(fi);
   for i=1:length(s2)+length(fi)  %  變換后的矩陣長度
      if (mod(i,2)==1)
      s(i)=s2((i+1)/2);  %  矩陣奇數下標為小波上一層(0,1,2,3)離散值
      else
      s(i)=fi(i/2);  %  矩陣偶數下標為小波下一層(1/2,3/2,5/2)(經過伸縮變換后)的離散值
      end
   end
   s2=s;
end
 
 

%  采用雙尺度方程求解小波基函數 PSI(T)
%  wavelet.m
function wave_base=wavelet(test,h);
MOMENT=10;  %  消失矩階數
LEFT_SCALET=0;  %  尺度函數左支撐區間
RIGHT_SCALET=2*MOMENT-1;  %  尺度函數右支撐區間
LEFT_BASIS=1-MOMENT;  %  小波基函數左支撐區間
RIGHT_BASIS=MOMENT;  %  小波基函數右支撐區間
MIN_STEP=1/512;  %  最小離散步長
LEVEL=-log2(MIN_STEP);  %  計算需要的層數(離散精度)
MAX_LEVEL=8;  %  周期小波最大計算層數

i=0;
for t=LEFT_BASIS:MIN_STEP:RIGHT_BASIS;  %  小波基支撐長度
    s=0;
    for n=1-RIGHT_SCALET:1-LEFT_SCALET  %  g(n)取值范圍
        if((2*t-n)>=LEFT_SCALET&(2*t-n)<=RIGHT_SCALET)  %  尺度函數判斷
          s=s+h(1-n+1)*(-1)^(n)*test((2*t-n)/MIN_STEP+1);  %  計算任意精度的小波基函數值 
        end
    end
    i=i+1;
    wave_base(i)=s;
end

 

 

 

 

 

 

 

一維數字濾波器filter():
        Y=filter(B, A, X)   由傳遞函數模型向量BA描述的濾波器對向量X中的元素進行濾波,並將結果數據存放在向量Y中。
      [Y, Zf]=filter(B, A, X, Zi)   給出了濾波器延時的初始和終止條件ZfZi
例子:
        人體心電信號在測量過程中往往受到工業高頻干擾,所以必須經過低通濾波處理后,才能判斷心臟功能的有用信息。下面給出一實際心電圖信號采樣序列樣本x(n),其中存在高頻干擾。在試驗中以x(n)作為輸入序列,濾除其中的干擾成分
{ x(n) } = {-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,
                 -2,-4,8,12,12,10,6,6,4,0,0,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0}
Matlab程序設計如下:
            X=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,4,0,0,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
               figure;
           plot(X);
           xlabel('時間');
           ylabel('幅值');
           wp=40;   ws=50;  rp=0.5;  rs=40;   Fs=200;
           [N, Wn] = buttord(wp/(Fs/2), ws/(Fs/2), rp, rs);
           [b, a]=butter(N, Wn);
           figure;
           [H, W]=freqz(b, a);
           plot(W*Fs/(2*pi), abs(H));   grid;
           xlabel('頻率/Hz');
           ylabel('幅值');
           Y=filter(b,a,X);
           figure
           plot(Y)
           xlabel('時間');
           ylabel('幅值');
           figure
           psd(X, [ ],200);
           figure
           psd(Y, [ ],200);

      end;

分析這段程序可知包括以下幾部分:
       1)首先繪制原始數據的圖形;
       2)設計一個Butterworth低通濾波器並繪制出它的幅頻響應曲線;
       3)用設計的濾波器對原數據進行濾波;
       4)繪制濾波以后的數據圖形;
       5)繪制原數據功率譜圖形;
       6)繪制濾波以后數據功率譜圖形。

       濾波器的主要目的是按照設計者的目的,突出或抑制一些頻段。在本程序中,設計了一個低通濾波器,主要是抑制高頻段突出低頻段;在心電圖信號分析中,要濾除工業高頻干擾,突出低頻部分.

 

  有時某些信號容易受到噪聲污染,導致無法直接辨別信號的發展趨勢。  由於信號的發展趨勢往往代表信號的低頻部分,因此通過信號的多尺度分解,在分解的低頻系數中可以觀察到信號的發展趨勢。

  由於噪聲的污染,從原始信號x中無法觀察信號的發展趨勢。通過進行五尺度的小波分解,在小波分解的低頻系數重構中可以明顯地看到原始信號的發展趨勢。這是因為信號的發展趨勢往往是信號的低頻成分,在小波變換中對應着最大尺度小波變換的低頻系數。此外還可以在低頻中理解它,在進行低頻成分的尺度分解時,隨着分解層數的增加,它所含的高頻成分會隨之減少,因此隨着尺度的增加,更多高頻的信號被濾掉,可以看到信號的發展趨勢。

 

1.監測信號的自相似性
       直觀上講,小波分解系數表示了信號與小波之間的“相似指數”,如果相似程度越高,則相似指數越大。因此如果一個信號的不同的尺度之間相似,則小波系數在不同的尺度上也應該相似。因此可以通過小波分解檢測信號的自相似性,即檢測信號的分形特征。實踐表明,通過小波分解可以很好地研究信號或圖像的分形特征。
       下面通過一個簡單的例子來說明小波分析在檢測信號自相似性中的應用,待檢測的信號是經過反復迭代生成的信號,因此具有自相似性。

       程序代碼如下:
       load vonkoch;
       x=vonkoch;
       subplot(211);
       plot(x);
       title('原始信號');
       subplot(212);
       %進行一維連續小波變換
      f=cwt(x,[2:2:128],'coif3','plot');
     從圖中可以看出分解后的小波系數在許多尺度上看上去都非常相似。


2.信號的奇異性檢測
       信號的突變點和奇異點等不規則部分通常包含重要信息。
       一般信號的奇異性分為兩種情況:(1)信號在某一時刻其幅值發生突變,引起信號的非連續,這種類型的突變稱為第一類型的間斷點;   2)信號在外觀上很光滑,幅值沒有發生突變,但是信號的一階微分有突變發生且一階微分不連續,這種類型的突變稱為第二類型的間斷點。
       應用小波分析可以檢測出信號中的突變點的位置、類型以及變化的幅度。下面介紹小波分析在信號奇異性檢測中的應用。

      1)第一類型間斷點的檢測
       下面舉例說明小波分析用於檢測第一類型的間斷點。

 
       在本例中,信號的不連續是由於低頻特征的正弦信號在后半部分突然有高頻特征的正弦信號加入,首先利用傅里葉變換分析對信號在頻域進行分析,發現無檢測突變點,接着利用小波分析進行分析,結果證明它能夠准確地檢測出了信號幅值突變的位置,即高頻信號加入的時間點。
    
       程序代碼如下:
       load  freqbrk;
       x=freqbrk;
       
      %對信號進行傅里葉變換
       f=fft(x,1024);
       f=abs(f);
      
      figure;
      subplot(211);
      plot(x);
      subplot(212);
      plot(f);

      %使用db6小波進行6層分解
       [c,l]=wavedec(x,6,'db6');
       figure(2);
      subplot(811);
      plot(x);
      ylabel('x');
   %對分解的第六層低頻系數進行重構
     a=wrcoef('a',c,l,'db6',6);
     subplot(812);
     plot(a);
     ylabel('a6');
     for i=1:6
                %對分解的第6層到第1層的高頻系數分別進行重構
                d=wrcoef('d',c,l,'db6',7-i);
                subplot(8,1,i+2);
                plot(d);
                ylabel(['d',num2str(7-i)]);
   end

 


        由圖中可以看出,由於傅里葉變換不具有時間分辨力,因此無法檢測信號的間斷點。而在小波分析的圖中,在信號的小波分解的第一層高頻系數d1和第二層高頻系數d2中,可以非常清楚地觀察到信號的不連續點,用db1小波比用db6小波要好。
       這個例子也表明小波分析在檢測信號的奇異點時具有傅里葉變換無法比擬的優越性,利用小波分析可以精確地檢測出信號的突變點。


      在信號處理中,信號中通常都包含噪聲,而噪聲的存在增加了辨別信號不連續點的難度。一般來說,如果信號小波分解的第一層能夠估計出噪聲的大體位置,則信號的間斷點就能夠在小波分解的更深層次上表現出來。

 

下面通過例子說明如何應用小波分析識別某一頻率區間上的信號:

 

 

        在本例中,使用小波分析一個由三個不同頻率的正弦信號疊加的信號,看是否能將這三個正弦信號區分開來,結果證明小波分析可以很好地識別某一頻率區間的信號。

        程序代碼如下:
        load sumsin;
        x=sumsin;
       figure;
       subplot(611);
       plot(x);
       ylabel('x');
       title('原始信號以及各層近似信號');
      
      %使用db3小波進行5層分解
       [c,l]=wavedec(x,5,'db3');
       for i=1:5
                %對分解的第5層到第1層的低頻系數分別進行重構
                a=wrcoef('a',c,l,'db3',6-i);
                subplot(6,1,i+1);
                plot(a);
                ylabel(['a',num2str(6-i)]);
       end
 

       figure;
       subplot(611)
       plot(x);
       ylabel('x')
       for i=1:5
              %對分解的第5層到第1層的高頻系數進行重構
               d=wrcoef('d',c,l,'db3',6-i);
               subplot(6,1,i+1);
               plot(d);
               ylabel(['d',num2str(6-i)]);
       end

 

 

分析:
       在本例中,該信號是由周期分別為200202的信號組成的,它們的采樣周期均為1,為方便起見,在此分別稱為低頻、中頻和高頻的正弦信號。從圖中可以看出,低頻、中頻和高頻信號分別對應於分解的近似信號a4、細節信號d4以及細節信號d1


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM