頻譜分析演算法

⑴ 有什麼演算法處理生物信號

頻譜分析演算法處理生物信號。

處理生物信號產品特點：

1、外置式，採用USB2.0 全速傳輸方式，即插即用,方便於科研與教學。

2、4個通道性能指標完全一致的光電隔離放大器，硬體參數全程式控制可調。

3、交、直流具有相同的增益：量程±0.5V——±20μV。

處理生物信號系統原理：

由於生物信號種類繁多，信號的強弱不一（有些生物電信號非常微弱，比如兔減壓神經放電，其信號強度為微伏級，如果不進行信號的前置放大，根本無法觀察），頻率混疊（由於在生物信號中夾雜有眾多聲、光、電等干擾信號。

比如電網的50Hz信號，這些干擾信號的幅度往往比生物電信號本身的強度還要大，如果不將這些干擾信號濾除掉，那麼可能會因為過大的干擾信號致使有用的生物機能信號本身無法觀察），因此信號採集前往往需要放大和濾波處理。

⑵ 基準面旋迴的小波和頻譜分析

高解析度層序地層學理論的核心內容是「在基準面旋迴變化過程中，由於可容納空間與沉積物補給通量比值的變化，相同沉積體系域或相域中發生沉積物的體積分配作用和相分異作用所導致的沉積物保存程度、地層堆積樣式、相序、相類型及岩石結構和組合類型發生的規律性變化」(Cross，1994)。由於基準面的變化是海平面、構造沉降、沉積物補給、沉積負荷補償、沉積壓實與沉積地形等各種因素變化的綜合反映，因此，通過碎屑岩的厚度變化、粒度大小、有機質含量和沉積物類型及岩石的結構構造特徵表現出來，而這些地質現象又被高解析度的測井曲線記錄下來，這就為利用數學方法定量分析旋迴信息提供了依據。隨著計算機技術的發展，數字信號處理、地學信息系統分析使得地質問題的定量化成為地學研究的熱點之一。頻譜分析和小波分析是進行基準面旋迴定量分析的重要技術方法之一，二者均能從復雜的疊加信號提取相對獨立的天文周期信息，為基準面旋迴的分析提供重要的理論基礎和技術支持。

一、小波分析原理

小波是一個衰減的波形，它在有限的區域里存在(不為零)，且其均值為零。小波是尖銳變化而且是不規則的波形，因此用小波能更好地刻畫信號的局部特徵。小波變換是小波分析的核心，設測井曲線為f(h)，小波變換的演算法如下:

0

公式中ωf(i，j)為尺度i下刻度j處的小波系數;a為尺度參數，a=2－j;b=a·j。

採用低通和高通濾波器，求取低頻系數CA₁和高頻系數CD₁，然後再分解CA₁為CA₂和CD₂，再分解CA₂為CA₃和CD₃，如此類推(圖4-16)，然後再對低頻系數和高頻系數進行重構高頻信號D₃、D₂、D₁和低頻信號A₃，分別得到不同頻率的周期旋迴。

圖4-16一維小波變換三次分解與重構示意圖

二、小波分析在高解析度層序地層劃分中的意義

1.識別轉換基準面和剝蝕面

小波變換是傅里葉變換的發展，其實質是引入伸縮、平移思想，對不同頻率成分自動地選取時域和取樣步長，從而能夠聚焦到物體的任意微小細節，因此小波變換被譽為「數學顯微鏡」。利用小波變換能夠有效地檢測包含非常重要的信息的奇異點和不規則點的信號。在測井信號中，奇異點常常代表地層信息變化劇烈點，如具有特殊意義的剝蝕面和轉換基準面。因此利用小波變換的這一特點，可以很容易地識別測井曲線的劇變點以及劇變程度，分析出對於高解析度層序劃分具有特殊意義的剝蝕面和轉換基準面，從而確定層序邊界，如圖4-17所示為鄂爾多斯盆地A16井延長組與富縣組的分界線的位置為不整合界線，經小波變換後檢測包含明顯信息的奇異點信號，同樣在湖泛面的位置也檢測出不規則的信號。

圖4-17A16井小波分析的延長組與富縣組界面和旋迴周期性識別

2.識別不完整基準面旋迴周期

在地質過程中，存在眾多的侵蝕、沉積物過而不沉和欠補償沉積等作用，導致地質記錄和保存的不完整性，測井信號所反映出的周期性、旋迴性必然是殘缺不全的。而傳統的頻譜分析方法，如經典的傅里葉變換，是對信號在整個時間域上的周期性進行分析，而不能彌補這些不完整的旋迴記錄。而小波分析對於地層中發育的單獨的上升半旋迴、下降半旋迴和不對稱旋迴，可以實現對信號在時間域和頻率域局部化分析，又如圖4-17所示，在鄂爾多斯盆地A16井1860～2020m段的自然伽馬測井曲線中，直觀的分析曲線中以正旋迴為主，逆旋迴難以觀察，經小波分析後，可以明顯的分析出該井段的測井曲線包含2.5個較高頻旋迴，由此來看，小波分析為基準面旋迴的識別提供了技術支持。

3.識別測井曲線的隱含周期

高解析度等時地層對比的關鍵是識別地層記錄中這些代表多級次基準面旋迴的地層旋迴，在基準面旋迴的識別中，測井序列是迄今為止所能獲得的解析度最高、連續性最好的地質數據，其中蘊藏著豐富的地質信息。不同的測井數據在不同程度上記錄著地質演化的歷史，從不同側面反映著地層形成演化的條件和影響因素，如海平面變化、古環境、古地理、古氣候信息及其變化情況等，是了解地質過程的最好工具。但由於地質過程的多時間多尺度特徵和各種串級過程，地層沉積序列實際上是各種地質周期沉積響應的疊加，再加上不確定因素和局部因素引起的隨機波動的干擾，從測井曲線上很難直觀地從測井信號中判讀各種隱含周期。

圖4-18疊加信號的小波分析分解圖

小波分析技術可以把測井信號分解在任意精度的不同頻帶內(圖4-18)，根據感興趣的信號頻帶范圍，把信號在一定尺度上分解，從而提取相應頻帶的信息，得到相應的地質周期(余繼峰等，2003)，對信號中的低頻成分，採用寬的時間窗，得到低的頻率解析度;對信號中的高頻成分，採用窄的時間窗，得到高的頻率解析度(Daubechies，1991)。由此可以利用小波分析的這種自適應特徵對測井曲線進行多尺度分析，選取信號中代表地質長周期的低頻部分來確定大的層序地層旋迴，中等頻率用來確定中等的地層層序，高頻部分則代表短周期的旋迴，可以用於小層的精細對比和劃分。如圖4-18中的A曲線是一個由3個正弦函數f(t)=sin2πt+sin3πt+sin4πt組成的疊加信號，經過小波分析後能夠將其分解為3個獨立的信號，如圖4-18中的曲線B、C和D，由此可見，小波分析可以用於測井曲線中隱含的旋迴周期的識別。

三、頻譜分析原理

1.頻譜分析基本原理

快速傅里葉變換法(FFT)是頻譜分析技術研究周期性最為常用的一種統計分析方法，主要通過對一復合波系進行傅里葉變換，將其分解成若干振幅和相位不同的間諧波，並找出其中振幅最大的波，即該復合波中的主要頻率。

傅里葉變換函數常以連續函數得出，若變換函數為x(t)，則傅里葉變換由下式:

0

式中:t為時間;f代表頻率。

在實際過程中，采樣都是離散和有限的，如果是連續信號，在應用計算機處理時也需要進行截斷並離散化，因此在實際數據處理時，都採用離散傅里葉變換，其變換公式如下:

假定有一段N項離散時間序列x_m，其離散傅里葉變換為

0

其中X_k稱為頻譜值，k=0，1，2，…，N－1。

快速傅里葉變換(FFT)是計算離散傅里葉變換(DFT)的一種快速演算法，能迅速提高DFT的運算速度。FFT的演算法種類很多，基於FFT演算法中的頻率抽選法(Cara，1982)是其中的常用的一種。

當k為偶數時，公式為

0

k=0，1，…，A－1;A=N/2，此處，k已用2k所代替。

當k為奇數時，用2k+1代替k，可以得到:

0

k=0，1，…，A－1;A=N/2。

由於頻譜值X_k沒有直觀意義，習慣上總把其轉化成能量，構成直觀的能量-頻率圖，對應於頻率f_k的諧振能量為

P_k=(real(X_k))²+(imag(X_k))²

其中，k=0，1，2，…，N－1。

2.頻譜分析的Matlab實現

對於大多數地質工作者來說，要進行上述復雜的頻譜分析數學計算，可能會有一定的難度。由Mathwork公司推出的Matlab軟體集數字分析、矩陣運算、信號處理和圖像處理、顯示於一體，構成了一個方便靈活的、界面友好的用戶環境。在這個環境下，對所求解的問題用戶只需簡單地列出數學表達式，其結果便以數值或圖形的方式顯示出來。以下是應用Matlab實現復合信號頻譜分析進行信號分解的一個實例，以期為Matlab在地質中的應用起到拋磚引玉的作用。

利用Matlab中的快速傅里葉變換函數fft()實現頻譜分析的實例如下:

0

對上述信號signal，利用Matlab進行快速傅里葉變換:

R=fft(signal，N)

式中:R為測井信號的快速傅里葉變換頻譜值序列，fft為快速傅里葉變換函數，N是signal的長度。

頻譜值R對應於頻率f_k的諧振能量為

P_f(k)=R_k*conj(R_k)

式中:P_f(k)為頻率f_k的諧振能量值，R_k頻率f_k的快速傅里葉變換(FFT)值，conj(R_k)為R_k的共軛復數。由於所求功率為模數，頻譜曲線左右對稱，因此僅取其中的一半進行計算，所以k=0，1，2…Round(N/2)，Round(N/2)是對N/2進行取整。

0

運行結果如圖4-19所示。

圖4-19復合正弦曲線的Matlab頻譜分析

上述實例可以看出，使用Matlab進行數據處理十分方便與靈活。Matlab為用戶提供了大量的功能函數，可以為研究人員避免大量重復性的數學運算、而把更多的精力集中到專業的方法研究中提供了便利。

對於測井曲線頻譜分析，是測井曲線的數字化處理，這些與測井相關的問題正是Matlab很方便解決地。

3.頻譜分析在高解析度層序地層劃分中的意義

(1)天文周期的確定。在六個級次的基準面旋迴級次劃分方案中所強調的主控制因素各不相同，以陸相盆地為例，低頻長周期的旋迴如超長期和長期旋迴主要受構造作用控制，而高頻短周期旋迴如中期、短期和超短期基準面旋迴則分別受天文因素的偏心率長周期、偏心率短周期和歲差周期所控制(鄭榮才等，2001)。因此，如何從復雜的地層信息中識別出保存有天文周期的信息是進行高頻基準面旋迴分析的重要內容，而測井曲線的頻譜分析可以從地層中獲取包含不同級次的天文周期信息，從而更有效地進行高頻周期旋迴控制因素的分析。

(2)估算地層沉積速率。通過頻譜分析所得到的天文周期旋迴，在地質歷史時期，組成米蘭科維奇旋迴的偏心率周期和歲差周期是穩定的，因此在已知沉積周期厚度和周期持續的時間，就可以得到相應的沉積速率。

⑶ matlab怎樣進行頻譜分析

姓名：張猛

【嵌牛導讀】：如何對一個信號畫出頻譜並進行分析，從頻譜中得到有用的信息

引用：http://blog.sina.com.cn/s/blog_a07f4fe301013gj3.html

【嵌牛鼻子】：matlab  fft  頻域

【嵌牛提問】：如何畫頻譜，對頻譜如何分析

【嵌牛正文】

圖像的頻率是表徵圖像中灰度變化劇烈程度的指標，是灰度在平面空間上的梯度。如：大面積的沙漠在圖像中是一片灰度變化緩慢的區域，對應的頻率值很低；而對於地表屬性變換劇烈的邊緣區域在圖像中是一片灰度變化劇烈的區域，對應的頻率值較高。傅立葉變換在實際中有非常明顯的物理意義，設f是一個能量有限的模擬信號，則其傅立葉變換就表示f的譜。從純粹的數學意義上看，傅立葉變換是將一個函數轉換為一系列周期函數來處理的。從物理效果看，傅立葉變換是將圖像從空間域轉換到頻率域，其逆變換是將圖像從頻率域轉換到空間域。換句話說，傅立葉變換的物理意義是將圖像的灰度分布函數變換為圖像的頻率分布函數，傅立葉逆變換是將圖像的頻率分布函數變換為灰度分布函數。

這樣通過觀察傅立葉變換後的頻譜圖，也叫功率圖，我們首先就可以看出，圖像的能量分布，如果頻譜圖中暗的點數更多，那麼實際圖像是比較柔和的（因為各點與鄰域差異都不大，梯度相對較小），反之，如果頻譜圖中亮的點數多，那麼實際圖像一定是尖銳的，邊界分明且邊界兩邊像素差異較大的。對頻譜移頻到原點以後，可以看出圖像的頻率分布是以原點為圓心，對稱分布的。將頻譜移頻到圓心除了可以清晰地看出圖像頻率分布以外，還有一個好處，它可以分離出有周期性規律的干擾信號，比如正弦干擾，一副帶有正弦干擾，移頻到原點的頻譜圖上可以看出除了中心以外還存在以某一點為中心，對稱分布的亮點集合，這個集合就是干擾噪音產生的，這時可以很直觀的通過在該位置放置帶阻濾波器消除干擾。另外我還想說明以下幾點：

1、圖像經過二維傅立葉變換後，其變換系數矩陣表明：

若變換矩陣Fn原點設在中心，其頻譜能量集中分布在變換系數短陣的中心附近(圖中陰影區)。若所用的二維傅立葉變換矩陣Fn的原點設在左上角，那麼圖像信號能量將集中在系數矩陣的四個角上。這是由二維傅立葉變換本身性質決定的。同時也表明一股圖像能量集中低頻區域。

2 、變換之後的圖像在原點平移之前四角是低頻，最亮，平移之後中間部分是低頻，最亮，亮度大說明低頻的能量大（幅角比較大）。

從計算機處理精度上就不難理解，一個長度為N的信號，最多隻能有N/2+1個不同頻率，再多的頻率就超過了計算機所能所處理的精度范圍）

X[]數組又分兩種，一種是表示餘弦波的不同頻率幅度值：Re X[]，另一種是表示正弦波的不同頻率幅度值：Im X[]，Re是實數(Real)的意思，Im是虛數(Imagine)的意思，採用復數的表示方法把正餘弦波組合起來進行表示，但這里我們不考慮復數的其它作用，只記住是一種組合方法而已，目的是為了便於表達（在後面我們會知道，復數形式的傅立葉變換長度是N，而不是N/2+1）。

用Matlab實現快速傅立葉變換

FFT是離散傅立葉變換的快速演算法，可以將一個信號變換到頻域。有些信號在時域上是很難看出什麼特徵的，但是如果變換到頻域之後，就很容易看出特徵了。這就是很多信號分析採用FFT變換的原因。另外，FFT可以將一個信號的頻譜提取出來，這在頻譜分析方面也是經常用的。

雖然很多人都知道FFT是什麼，可以用來做什麼，怎麼去做，但是卻不知道FFT之後的結果是什意思、如何決定要使用多少點來做FFT。

現在就根據實際經驗來說說FFT結果的具體物理意義。一個模擬信號，經過ADC采樣之後，就變成了數字信號。采樣定理告訴我們，采樣頻率要大於信號頻率的兩倍，這些我就不在此啰嗦了。

采樣得到的數字信號，就可以做FFT變換了。N個采樣點，經過FFT之後，就可以得到N個點的FFT結果。為了方便進行FFT運算，通常N取2的整數次方。

假設采樣頻率為Fs，信號頻率F，采樣點數為N。那麼FFT之後結果就是一個為N點的復數。每一個點就對應著一個頻率點。這個點的模值，就是該頻率值下的幅度特性。具體跟原始信號的幅度有什麼關系呢？假設原始信號的峰值為A，那麼FFT的結果的每個點（除了第一個點直流分量之外）的模值就是A的N/2倍。而第一個點就是直流分量，它的模值就是直流分量的N倍。而每個點的相位呢，就是在該頻率下的信號的相位。第一個點表示直流分量（即0Hz），而最後一個點N的再下一個點（實際上這個點是不存在的，這里是假設的第N+1個點，也可以看做是將第一個點分做兩半分，另一半移到最後）則表示采樣頻率Fs，這中間被N-1個點平均分成N等份，每個點的頻率依次增加。例如某點n所表示的頻率為：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到頻率為Fs/N，如果采樣頻率Fs為1024Hz，采樣點數為1024點，則可以分辨到1Hz。1024Hz的采樣率采樣1024點，剛好是1秒，也就是說，采樣1秒時間的信號並做FFT，則結果可以分析到1Hz，如果采樣2秒時間的信號並做FFT，則結果可以分析到0.5Hz。如果要提高頻率分辨力，則必須增加采樣點數，也即采樣時間。頻率解析度和采樣時間是倒數關系。

假設FFT之後某點n用復數a+bi表示，那麼這個復數的模就是An=根號a*a+b*b，相位就是Pn=atan2(b,a)。根據以上的結果，就可以計算出n點（n≠1，且n<=N/2）對應的信號的表達式為：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。對於n=1點的信號，是直流分量，幅度即為A1/N。由於FFT結果的對稱性，通常我們只使用前半部分的結果，即小於采樣頻率一半的結果。

下面以一個實際的信號來做說明。假設我們有一個信號，它含有2V的直流分量，頻率為50Hz、相位為-30度、幅度為3V的交流信號，以及一個頻率（f0）為75Hz、相位為90度、幅度為1.5V的交流信號。用數學表達式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)。式中cos參數為弧度，所以-30度和90度要分別換算成弧度。我們以256Hz的采樣率對這個信號進行采樣，總共采樣256點。按照我們上面的分析，Fn=(n-1)*Fs/N，我們可以知道，每兩個點之間的間距就是1Hz，第n個點的頻率就是n-1。我們的信號有3個頻率：0Hz、50Hz、75Hz，應該分別在第1個點、第51個點、第76個點上出現峰值，其它各點應該接近0。實際情況如何呢？我們來看看FFT的結果的模值如圖所示。

從圖中我們可以看到，在第1點、第51點、和第76點附近有比較大的值。我們分別將這三個點附近的數據拿上來細看：1點：512+0i2點：-2.6195E-14 - 1.4162E-13i3點：-2.8586E-14 - 1.1898E-13i50點：-6.2076E-13 - 2.1713E-12i51點：332.55 - 192i52點：-1.6707E-12 - 1.5241E-12i75點：-2.2199E-13 -1.0076E-12i76點：3.4315E-12 + 192i77點：-3.0263E-14 +7.5609E-13i很明顯，1點、51點、76點的值都比較大，它附近的點值都很小，可以認為是0，即在那些頻率點上的信號幅度為0。接著，我們來計算各點的幅度值。分別計算這三個點的模值，結果如下：1點：51251點：38476點：192按照公式，可以計算出直流分量為：512/N=512/256=2；50Hz信號的幅度為：384/(N/2)=384/(256/2)=3；75Hz信號的幅度為192/(N/2)=192/(256/2)=1.5。可見，從頻譜分析出來的幅度是正確的。然後再來計算相位信息。直流信號沒有相位可言，不用管它。先計算50Hz信號的相位，atan2(-192, 332.55)=-0.5236,結果是弧度，換算為角度就是180*(-0.5236)/pi=-30.0001。再計算75Hz信號的相位，atan2(192, 3.4315E-12)=1.5708弧度，換算成角度就是180*1.5708/pi=90.0002。可見，相位也是對的。根據FFT結果以及上面的分析計算，我們就可以寫出信號的表達式了，它就是我們開始提供的信號。

總結：假設采樣頻率為Fs，采樣點數為N，做FFT之後，某一點n（n從1開始）表示的頻率為：Fn=(n-1)*Fs/N；該點的模值除以N/2就是對應該頻率下的信號的幅度（對於直流信號是除以N）；該點的相位即是對應該頻率下的信號的相位。相位的計算可用函數atan2(b,a)計算。atan2(b,a)是求坐標為(a,b)點的角度值，范圍從-pi到pi。要精確到xHz，則需要采樣長度為1/x秒的信號，並做FFT。要提高頻率解析度，就需要增加采樣點數，這在一些實際的應用中是不現實的，需要在較短的時間內完成分析。解決這個問題的方法有頻率細分法，比較簡單的方法是采樣比較短時間的信號，然後在後面補充一定數量的0，使其長度達到需要的點數，再做FFT，這在一定程度上能夠提高頻率分辨力。具體的頻率細分法可參考相關文獻。

附貼上上述例子的matlab程序：

Matlab的例子（一）

t=0:1/256:1;%采樣步長

y= 2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180);

N=length(t); %樣點個數

plot(t,y);

fs=256;%采樣頻率

df=fs/(N-1);%解析度

f=(0:N-1)*df;%其中每點的頻率

Y=fft(y)/N*2;%真實的幅值

%Y=fftshift(Y);

figure(2)

plot(f,abs(Y));

由於以上程序是結合傅里葉演算法轉換得到的對稱圖，而常用的只需要一半就可以了。對應的程序如下：

t=0:1/256:1;%采樣步長

y= 2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180);

N=length(t); %樣點個數

plot(t,y);

fs=256;%采樣頻率

df=fs/(N-1);%解析度

f=(0:N-1)*df;%其中每點的頻率

Y=fft(y(1:N))/N*2;%真實的幅值

%Y=fftshift(Y);

figure(2)

plot(f(1:N/2),abs(Y(1:N/2)));