點擊下方卡片,關注「新機器視覺」公眾號
重磅乾貨,第一時間送達
轉自:阿柴的算法學習日記
開門見山,變換的本質就是將原有的信號分解為最簡單的表示。這個最簡單的表示抽象出來就是基。基在向量空間中的表現就是一組兩兩無關的向量,向量空間中任一一個向量都可以用它們表示;在函數空間中的變現就是一組兩兩無關的函數,函數空間的任一一個函數都可以用它們表示。舉一個很形象的例子,學計算機的人都知道計算機的任意文件在電腦中都是以0和1兩個數字組成的,那麼我們就可以把0、1稱為文件空間的基。進一步,我們還得定義何為兩兩無關。在此,我們只看與我們息息相關的向量空間與函數空間。在向量空間與函數空間中的兩兩無關被定義為內積為0,也就是下圖所示(上面為向量內積的定義,下面為函數內積的定義):
有了這些基礎,我們便可以開始去理解什麼是傅里葉變換。
我們首先來看傅里葉級數,它是傅里葉變換的基礎。傅里葉級數是法國數學家傅里葉發現的,他認為任何連續周期函數都可以由一組適當的正弦曲線組合而成。事實上從現在來看,這句話是錯誤的。它必須加上前提條件,也就是所謂的狄利克雷條件,但即便如此,這個發現在當時也是一個巨大的突破,它提供給人們一個全新視角來看待信號。由此,我們可以說:滿足狄利克雷條件的周期信號,可以分解成一組適當的正餘弦信號集。狄利克雷條件包括:
其實實際上大部分周期信號都滿足狄利克雷條件,也就是說我們平常接觸到的周期信號大多數都是可以進行傅里葉展開的。傅里葉級數實際上就是把這個空間中的一個向量通過基的線性組合的方式寫出來。
我們已經知道了,兩個正交的函數內積為 0 ,因此如果一個函數是由正交基的線性組合表示,我們可以很容易地求得它的係數。
設:
這就是傅里葉級數,實際上就是基的線性組合,1,.....,,......,
就是這個函數空間的基。它們都是正交的也就是說它們兩兩內積為0。
是 1 的係數,
和
分別是
和
的係數。現在我們想計算每一個基的係數。為此,我們試着計算一下它和其中一項的內積,比如

不過注意到不同的函數都是正交的,所以非零項只有.也就是說,
。這麼說似乎還是不明顯,我這麼寫:

類比一下歐幾里得空間。比如說二維歐式空間中,我取一組正交基
,那麼向量
寫成這組基的線性組合時係數分別是多少呢?答案就是


實際上從空間的角度,二者確實沒有什麼本質區別。
那麼現在我們已經理解了,傅里葉基其實就是這個函數構成的空間中的一組正交基,而傅里葉級數就是把空間裡的元素寫成基的線性組合。但我們注意到這裡的函數是的,如果換一個區間,結果會如何呢?
2.2、傅里葉變換
我們把區間換成,這樣相當於把每個函數「拉伸」了
,這樣每個基也拉伸了
倍,變成了
或
,
.
而傅里葉級數就變成了


後面一種形式根據歐拉公式由前一種形式簡化而來
下面我們考慮指數函數形式的傅里葉級數在時的情形。
1)在區間變為後,傅里葉係數變為:
2)在的傅里葉展開中,把傅里葉係數帶入,得到
3)上式的積分式是關於變量的,而外邊的項
與
無關。因此可以寫成
4)我們現在希望上式是黎曼和的形式,這樣就可以寫成定積分。為此取,這樣上式變成
5)當時,
,上式作為黎曼和的形式可以寫成定積分:
6)上式略作變形,得到
這樣,我們從傅里葉級數出發,在的情形下得到了一個式子。這個式子的本質仍然是一個求和式,只是因為極限所以寫成了積分式。它把
寫成了
的線性組合。等等!線性組合?!
確實如此啊。只是這裡的變量是取遍整個數軸的。相當於對實數軸上的每個點都對應了這個函數的一個「基」,而這個積分式就是這些基的「線性組合」。
括號里的那堆東西,也就是

就叫做這個函數的傅里葉變換。
傅里葉變換是這樣一個函數,它在處的函數值
表示函數
在
對應的基上的係數。至此我們就完成了傅里葉變換從空間角度的介紹。
2.3、傅立葉變換的可視化理解
傅里葉變換的形式為:(此式與上一小節推導出的式子本質上相同,只是去掉了前面的係數,對後面的討論沒有影響),我們知道,根據歐拉公式,
。也就是說,傅里葉變換的本質就是:將原始信號乘上一組三角函數(正餘弦),之後在整個時間域上積分。就這麼簡單!
你可能也聽過,傅里葉變換是將函數從從時域轉為頻域。那麼,為什麼這麼做就能將信號從時域轉為頻域呢?
我們來看一個信號:y = sin(3t),如下圖:

很好的周期性質,且每個周期的積分值都是0,是吧?如果對這個函數在積分,那就是基本是0,因為
包含了無數個周期。PS:雖然這個積分在高數上不可積,但是你應該明白這裡我要表達的意思:因為良好的周期性,且每個周期積分值是0,那麼最後在很長的一段時間區間上積分,得到的還是一個很小的數,近似為0。
我們來用一段較長的時間區間計算一下,,符合我們的預計。
現在,我們來將這個信號乘上一個sin(4t) ,則信號變為y1 = sin(3t)*sin(4t),如下圖:

具有一個較短的小幅震動的周期和一個較長的主體周期,對吧?且每個主體周期的積分值都是0。同以上討論,如果對這個函數在積分,基本還接近於0,因為
包含了無數個主體周期。
我們來用一段較長的時間區間計算一下,,也比較符合我們的預計。
之後呢,我們來將這個信號乘上一個sin(3.1t) ,則信號變為y2 = sin(3t)*sin(3.1t),如下圖:

同樣是有一個較短的小幅震動的周期和一個較長的主體周期,對吧?可以判斷,每個主體周期的積分值都是0(雖然(0,50)這個區間沒有完整地展示主題周期)。那麼,依然同以上討論,如果對這個函數在積分,基本還接近於0,因為
包含了無數個主體周期。
我們來用之前的時間區間計算一下,。
咦?這一次怎麼距離0這麼遠了呢?
原因就是:對於sin(3t)*sin(4t),它的主體周期較短,(0,50)是包含了好幾個主體周期的,也就是說(0,50)在某種程度上是類似於的。但是,對於sin(3t)*sin(3.1t),它的主體周期很長,(0,50)連它的一個完整的主體周期都沒有包含,那麼(0,50)是不能類似於
的,積分值自然就比較大。
我們此時可以這樣小小總結一下,對於信號y = sin(3t),它的頻率是3rad/s,(如果你喜歡用HZ,那就除以,就是
HZ,這裡使用rad/s,是為了與前面的傅里葉變換的公式中的w一致),而sin(4t)的頻率是4rad/s,sin(3.1t)的頻率是3.1rad/s。如果在
積分,那麼y1 = sin(3t)*sin(4t),y2 = sin(3t)*sin(3.1t)的積分值都是0,也就是說,sin(4t)和sin(3.1t)在這裡是沒差別的。但是!!!如果在一個有限區間內積分,由於sin(3.1t)的頻率3.1rad/s,距離原信號y = sin(3t)的頻率3rad/s更近,那麼sin(3.1t)和sin(3t)的乘積,也就是y2 = sin(3t)*sin(3.1t)的積分的絕對值會更大,也就是會離0更遠。這裡已經顯示出一定的頻率選擇性了。
最後,讓我們請出我們今天的主角,將這個信號乘上一個自己同頻率的sin(3t) ,則信號變為y3 = sin(3t)*sin(3t),如圖:

Amazing!!!發現了什麼?良好的周期性?還有呢?由於乘上了自己,任何時間的幅值都大於等於0了!不再滿足周期內積分值為0這個條件了!那麼,此時,我們對這個信號在積分,就會得到一個非常非常大的數字。這個很大很大的數字就告訴你,這個信號和你乘的信號是同頻率的!這就是可以知道信號中具有哪些頻率部分了,不是嗎?
我們還是來用之前的時間區間計算一下,。是不是比其他的積分值都大了好多?
好了,我們已經知道,將一個信號乘上一個特定頻率的sin函數,在上積分,可以將信號中與sin函數同頻率的部分篩選出來。那麼,原則上講,只要乘上所有頻率的sin函數,並積分,不就知道原始信號中的所有頻率部分了嗎?
但是這樣做需要把所有頻率乘進去,做無數次計算哈!算不出來的。所以,我們將所乘的sin函數的頻率作為符號變量w,來進行積分,即:
注意:這裡的w只是一個符號,一個變量,這樣的話,就只需要做一次積分,可以計算了。
計算出來X(w)之後,想知道特定的頻率w0對應的積分值,直接將w0帶入X(w)就立馬得到積分值。也就知道原信號中是否含有這一頻率的部分了。
好了,我們推導的這個式子,是不是與傅里葉變換的式子:
很像了呢?
這就是傅里葉變換的原理!乘上帶有符號變量的sin、cos函數,並積分,就知道原始信號中的所有頻率部分啦!
下面這幅圖可以幫助我們更好的理解。
三、二維傅里葉變換-在圖像上的應用
那首先來看一個例子~
1)它到中點的距離描述的是頻率
2)中點到它的方向,是平面波的方向
3)那一點的灰度值描述的是它的幅值

註:上述講解來源見水印
我們現在來看兩個信號,如下圖:


這兩個信號都是由sin(t)和sin(5t)組成的,y1是先出現了sin(5t),再出現了sin(t),y2是先出現了sin(t),再出現了sin(5t)。
我們對它們進行FT,看看他們包含怎麼樣的頻率,如下圖:
y1,FT
y2,FT
Amazing!發現了什麼?變換後的結果是一模一樣的,都在w=1rad/s和w=5rad/s出現了峰值!這就可以說明FT的缺點了——FT只能提供頻域信息,而完全丟失了時域信息!!!
不管某一頻率的信號出現的時間是早還是晚,FT都是將它一視同仁地乘上sin和cos(FT的變換基函數),然後在整個時間區間加和。因此,它不能提供某一頻率信號出現的時間。
比如,對於上面兩個信號,FT只能告訴我們,它們都有1rad/s和5rad/s的頻率,而不能告訴我們1rad/s和5rad/s分別出現在哪個時間段。
所以,怎麼辦呢???
先別看下面,你自己想三想,你肯定可以想到的!
想一想——
想二想——
想三想——
好了,你想的是對的!
那就是把信號分成左右兩半啊!左邊進行一次FT,右邊進行一次FT,很簡單吧!好了,這就是短時傅里葉變換(STFT)的基本原理。
所以,接下來我們要正式開始步入—短時傅里葉變換(STFT),看看它是如何解決這個問題的。
如上所述,我們將信號從中間截斷,左邊進行一次FT,右邊進行一次FT,分別來看看。
y1左
y1右
y2左
y2右
可以看出,y1的左半部分是5rad/s,右半部分是1rad/s,y2恰好相反。這就說明,在y1中,(0, 25)的信號是5rad/s的頻率,(25, 50)的信號是5rad/s的頻率,y2恰好相反。這就是短時傅里葉變換的基本原理。
但是數學嘛,能用一個公式表達的,就別用一段話表達,截斷、切開這些語句太不專業了。截斷、切開的操作,更專業的講叫作分窗,接下來我們來看一看。
首先,你可以想象一下,有一個窗子在這個信號上從左向右滑動,每次你都只能看到這個信號的一部分,所以我們把這個長度叫作窗長width。
現在我們來定義一個方窗函數:
,如下圖,即是width = 10 的一個方窗函數:

定義了方窗函數之後,我們只需要對方窗函數進行平移,再與原信號作乘,就相當於原來的截斷、切開的操作,因此這種操作更專業地叫作分窗。
那麼,將方窗函數向右平移了(s可能是sliding的意思吧),再與原信號相乘,由於方窗函數除了中心的width部分是1外,其他部分都是0,這就相當於提取出了原信號在
處,寬度為width的部分,這個信號分窗這個操作就可以寫成:
。
如下兩圖所示,將與
相乘,就相當於取出來了
中的(20,30)中的一段。


那麼,我們對原信號中被提取出來的這一部分進行FT,就可以寫成:

PS:這裡之所以要變成
,是為了保證做FT的時候相乘的基函數具有統一性。
如此,變換後的代表原信號在
處、寬度為width的部分的傅里葉變換,也就可以提取出來原信號在
處、寬度為width的部分,包含各個頻率部分的多少!帶入不同的
,也就是隨着窗子的滑動,就可以知道不同的時間段內頻率的成分。
我們採用width為10的方窗函數對

進行STFT,如下:
首先,方窗函數位於處,與原始信號相乘,選擇出(0,10)的信號。

對選擇出來的信號進行FT。可以看到,當時,選擇的時間區間為(0,10),這一部分只包含了
的頻率成分。

之後,方窗函數向右移動,與原始信號相乘,選擇出不同時間區間的信號,進行FT。這裡選擇進行展示。


可以看到,當時,選擇的時間區間為(20,30),這一部分即包含了
的頻率成分,也包含了
的頻率成分。
重複以上過程,我們可以將方窗函數選擇的不同時間區間的信號的FT的結果拼合起來,形成一張三維圖。由此,我們即可知道,在的時間區間內,信號具有怎麼樣的頻率成分。

通過width = 10的方窗的STFT結果,我們可以知道,對於信號:

在(0,10)、(10,20)時間區間內,具有20rad/s的頻率成分;在(20,30)時間區間內,具有1rad/s和20rad/s的頻率成分;在(30,40)、(40,50)時間區間內,具有1rad/s的頻率成分。
最後,進行三點重要的討論。
第一點,變換之後的是一個三維函數,它有兩個自變量,
和w。
指的是原信號在
處,w上一篇文章我們已經討論過了,就是頻率。所以,STFT提取出來的信息就是:原信號在
處、寬度為width的部分,包含的頻率信息。
原則上講,可以得到任一對應的頻率成分,如下圖。

但是是連續的,並不意味着你知道了每個時刻的頻率成分,你知道的還只是
這一段區間內的頻率信息。所以一般不需要計算所有的
,每隔width計算一次即可。
你或許會想,我把width縮小一些,不就可以知道更精確的時間範圍內的頻率了嗎?是的,你的猜想很對!但是,如此做也會帶來一些頻域分辨率的問題。這一點涉及到一些時域分辨率和頻域分辨率的知識,我們下一篇文章會着重講。
第二點,方窗函數是可以包含入變換基函數內部的,這組成了新的基函數,同時反映了STFT的本質。
我們來看, 如果定義

那麼

那麼,STFT的公式:

就可以寫成:

我們在上一篇文章里說過,變換就是將原信號乘上一個基函數,再積分的過程,那麼,SDFT的基函數就是

Amazing!所以,STFT的本質是什麼呢?
STFT的本質就是將FT的基函數乘上一個方窗函數,形成了一個新的基函數
!前面說的分窗、截斷之類的都是表象,STFT的本質是基函數的改變!
那麼,為什麼STFT的基函數可以用於分窗,而FT的基函數不行呢?我們來看,我用正弦函數sin(5t)表示原來的基函數,那麼FT基函數和STFT基函數如下:


原因就是:FT的基函數是在時域無限延伸的,因此,無論怎麼平移,都是任分布在整個時域的,起不到分窗的作用。而STFT的基函數只在時域一段不為0,在剩下的時域都是0,因此,STFT的基函數的平移,就相當於自動加了窗子啦!
這種只在時域一段不為0,在剩下的時域都是0的性質被稱為「緊支撐性」(compactly supported),具有這種性質的函數,平移之後與一個信號相乘,就相當於分窗操作。這一點很重要,我們之後講小波變換的基函數的時候還會講。
第三點,我們前面對於分窗操作使用的函數一直稱為「方窗函數」,這是一種最理想的窗函數。還有一些其他的窗函數,比如,漢寧窗、海明窗、高斯窗等。窗函數本質都是一個窗子而已,原理是一模一樣的,上面所有的討論也都成立,只是這些窗子會讓信號稍稍變形一丟丟而已。你就想像方窗函數就是一面平面鏡,其他的窗函數就是哈哈鏡就行了。
六、短時傅里葉變換(STFT)的缺點(為什麼要提出小波變換)
3.1、分辨率問題
首先,我們需要了解一下海森堡測不準原理:,
為信號的時間不確定度,
為信號的頻率不確定度。即,我們永遠無法同時確定一個信號的確切時間和確切頻率。
原因比較簡單,頻率其實就是時域周期性。如果我只給你一個數據點,問你這個數據點的頻率是多少,這肯定是做不到的。要確定頻率,就需要一個時域區間(包含幾個時域周期)的信號。
時域區間越寬,信號的時間定位越不准,時間不確定度越大,但是得到的頻率越准,頻率不確定度
越小;我們稱之為:低的時域分辨率,高的頻域分辨率。
時域區間越窄,信號的時間定位越准,時間不確定度越小,但是得到的頻率越不准,頻率不確定度
越大;我們稱之為:高的時域分辨率,低的頻域分辨率。


如上兩圖,對於第一個圖中的信號,要確定頻率,即使把
都拿來,還是不太準,因為只包含了一個周期;對於第二個圖中
的信號,要確定頻率,取個
就差不多了,因為已經包含了好幾個周期了。
我們來總結一下。
對於低頻信號,為了更好地確定頻率,我們希望,時域區間寬一些,即時間不確定度大一些,根據海森堡測不準原理,頻率不確定度
自然小一些;即低頻信號,我們希望:寬窗子,低的時域分辨率,高的頻域分辨率。
對於高頻信號,為了更好地在時域定位,我們希望,時域區間窄一些,即時間不確定度小一些,根據海森堡測不準原理,頻率不確定度
自然大一些;即高頻信號,我們希望:窄窗子,高的時域分辨率,低的頻域分辨率。

上圖所示是我們希望的動態分辨率。圖中每個小矩形的軸方向的寬度是時間區間大小,
軸方向的寬度是頻率區間大小。注意,每個小矩形的面積是相等的,這保證了時域分辨率乘上頻域分辨率是定值,最大程度滿足海德堡測不準原理。通過圖可以看出,我們希望,對於低頻信號:低的時域分辨率,高的頻域分辨率;對於高頻信號:高的時域分辨率,低的頻域分辨率。
對於整體低頻、局部高頻的信號,這種動態調整分辨率的規則特別有用。在實際信號中,頻率非常高的高頻信號往往是一種噪聲,只在局部出現,基本都滿足整體低頻、局部高頻這一條件。
最後,我們再來看兩張分辨率圖來強化一下對於分辨率的理解。

上圖是一張採集信號的分辨率圖。每個小矩形的軸方向的寬度是很小,
軸方向的寬度很大。也就是說,其時域分辨率很好,可以確切地確定每個信號採樣點的時間,但是其頻域分辨率很差,或者說完全不具有頻域分辨率。

上圖是一張傅里葉變換(FT)的分辨率圖。每個小矩形的軸方向的寬度是很大,
軸方向的寬度很小。也就是說,其頻域分辨率很好,可以比較精確地確定信號中的頻率成分,但是其時域分辨率很差,或者說完全丟失了時域分辨率。
傅里葉變換的這一特性,這一點我在上一篇文章里講過,這也是我們轉而使用短時傅里葉變換(STFT)的原因。
3.2、短時傅里葉變換(STFT)的缺點
我們來回憶一下STFT,STFT的窗長是固定的,即時域分辨率是固定的,根據海德堡測不準原理,其頻域分辨率也是固定的。其分辨率圖如下:

每個小矩形的軸方向的寬度和
軸方向的寬度是恆定的!也就是說,不論高頻低頻,其時域和頻域分辨率都不可調,這與我們之前討論的「對於低頻信號:低的時域分辨率,高的頻域分辨率;對於高頻信號:高的時域分辨率,低的頻域分辨率」這一原則不符合。
這種不符合會帶來什麼後果呢?
如圖所示正弦信號,0~250ms:300HZ,250~500ms:200HZ ,500~750ms:100HZ , 750~1000ms:50HZ。

選擇一個較窄的窗子做STFT,結果如下:

當窗子較窄的時候,STFT的時域分辨率還行,但是頻域分辨率不佳。
我們選擇一個寬一些的窗子做STFT,結果如下:

當窗子較寬的時候,STFT的頻域分辨率很好,基本可以確定頻率,但是時間軸上開始出現交疊了,也就是時域分辨率下降了。
我們選擇一個更寬的窗子做STFT,結果如下:

當窗子更寬的時候,STFT的頻域分辨率非常好了,但是時域分辨率已經很差了,時間軸上出現了大規模的交疊現象。
我們來總結一下,對於STFT,如果窗子的寬度選擇合適,是可以得到時域和頻域分辨率都「還可以」的結果的(由於STFT的分辨率固定,只能說「還可以」,不能說「滿意」,因為我們最想要的是動態分辨率)。但是,在變換之前,我們也不知道選擇多寬的窗子是合適的。
這就是STFT的缺點:
1、時間和頻率分辨率都固定,不能隨着頻率的高低實現動態可調;
2、選擇一個合適的窗寬十分困難。
這也正是小波變換提出的動機。
相關資料:
[1] https://zhuanlan.zhihu.com/p/66246381
[2] https://zhuanlan.zhihu.com/p/68323379
本文僅做學術分享,如有侵權,請聯繫刪文。
