發表文章

DSP筆記6_功率頻譜密度 Power Spectral Density(PSD)_Periodogram(週期圖)最原始簡單的作法

圖片
事實上,在上一篇我們完全不依賴scipy手刻的PSD計算的python程式 它的演算法流程就是所謂的Periodogram(週期圖)最原始簡單的作法。 periodogram(週期圖)是一種用來估計訊號功率頻譜密度 使用單一區段全部資料來計算 PSD 雖然速度快,但具有對雜訊敏感且變異數較大的限制。 顯示訊號的功率如何分布在不同頻率上。 常用於signal processing、vibration analysis、PPG/EEG/ECG 分析 在數學上,對於一個長度為  N  的離散時間訊號  x [ n ] : 會用以下數學式子表示週期圖功率頻譜密度 推導過程 Step1.計算離散傅立葉轉換(實際應用場景中,可用快速傅立葉轉換更有效率做出來。) Step2.計算功率譜,也就是計算「振幅(幅度大小)平方」 Step3.除以樣本數N做正規化 原始periodogram功率頻譜密度 提供 每個樣本的功率估計值 。 原始 periodogram就是 最初版本、最簡單的 PSD 估計方法。 在實務上,通常使用 FFT 來有效率地進行計算。 原始periodogram功率頻譜密度特性 變異數 : 高;即使增加 \(N\),變異數也不會下降,因此結果通常會比較多雜訊干擾。 偏差(Bias)、頻譜洩漏(Spectral Leakage): 在有限長度訊號中的不連續現象,會使功率擴散到鄰近的頻率。 解析度(Resolution) Δ f = f s N ;較長的訊號可以提升解析度。 一致性(Consistency): 原始 periodogram是不一致的估計量 原始periodogram功率頻譜密度限制 High variance : 峰值會產生波動 Spectral leakage : 有限長度資料會造成頻譜能量向周圍頻率擴散。 Poor consistency : N→∞時,原始 periodogram 不會收斂。 Trade-off between bias and resolution : 加窗(Windowing)可以降低頻譜洩漏,但會使頻譜峰值變寬。 常見改善方法 加窗(Windowing): Hann、Hamming、Blackman 視窗可降低 頻譜洩漏(spectral leakage) 。 Welch 方法(Welch'...

DSP筆記5_功率頻譜密度 Power Spectral Density(PSD)_從時域訊號到幅度大小的平方再除以資料總長度

圖片
功率頻譜密度(Power Spectral Density, PSD) 用來描述訊號在不同頻率下的功率分布情形。 功率頻譜密度(Power Spectral Density, PSD) S x ( f ) S_x(f) 適用於 功率訊號(power signals) 以及 廣義平穩程序(wide-sense stationary processes) 。 定義, 透過自相關函數(Autocorrelation)根據( Wiener–Khinchin theorem,維納–辛欽定理 )可得到如下式子 另一種定義 (此公式十分重要) : 總功率與性質(Total Power and Properties) 換言之,總功率(Total Power)可以由功率頻譜密度 S x ( f ) S_x(f) 對所有頻率積分後得到。 把給定區間的功率譜密度全部加總起來 練習程式ver1.建立一個典型的10Hz正弦波訊號 import numpy as np import matplotlib.pyplot as plt # Sinusoid parameters fs = 100 # Sampling freq (Hz)每秒取樣 100 次。 t = np.arange( 0 , 1 , 1 /fs)#[0. 0.01 0.02 ... 0.98 0.99],每次增加 0.01。 f0 = 10 # signal frequency # Sinusoidal signal x = np.sin( 2 * np.pi * f0 * t) # Plot plt.figure(figsize=( 8 , 4 )) plt.plot(t, x) plt.title( "10 Hz Sinusoid" ) plt.xlabel( "Time (s)" ) plt.ylabel( "Amplitude" ) plt.grid( True ) plt.show() 現在產生的是: x ( t ) = sin ⁡ ( 2 π f 0 t ) x(t)=\sin(2\pi f_0t) 其中 f 0 = 10  Hz f_0=10\text{ Hz} 意思是一秒鐘會完成...

DSP筆記4_能量頻譜密度 Energy Spectral Density(ESD)_從時域訊號到幅度大小的平方

圖片
  能量頻譜密度 Energy Spectral Density(ESD) 用於描述訊號在不同頻率下能量分布情形。 用  S E ( f )  來表示,X(f)為原先時間域x(t)的訊號做傅立葉轉換的結果。 能量譜密度的總能量關係滿足如下式子,就稱作帕塞瓦爾定里Parseval's theorem)。 總能量會等同於對時域跟頻域函數從負無窮到正無窮取平方的積分 所以 練習程式ver1.矩形脈衝Energy Signal import numpy as np import matplotlib.pyplot as plt fs = 1000 #Sampling Frequency T = 1 # Duration (seconds) t = np.linspace( 0 , T, int (T*fs), endpoint= False ) #建立一個全為0的訊號 x = np.zeros( len (t)) x[ int ( 0.4 *fs): int ( 0.6 *fs)] = 1 # 0.4 ~ 0.6 秒設定為 1 # Plot plt.figure(figsize=( 8 , 4 )) plt.plot(t, x) plt.title( "Rectangular Pulse" ) plt.xlabel( "Time (s)" ) plt.ylabel( "Amplitude" ) plt.grid( True ) plt.show() 運行效果 可以先得知目前此時間訊號的函數區間定義如下 x ( t ) = { 1 , 0.4 ≤ t < 0.6 0 , 其他時間 x(t)= \begin{cases} 1, & 0.4\leq t<0.6\\ 0, & \text{其他時間} \end{cases} 也就是一個寬度 0.6 − 0.4 = 0.2  秒 0.6-0.4=0.2\text{ 秒} 的矩形脈衝,是一個典型的 Energy Signal。。 它的能量為: E = ∫ − ∞ ∞ ∣ x ( t ) ∣ 2 d t E=\int_{-\infty}^{\infty}|x(t)|^2dt ...

Stochastic(Random) Processes筆記1_Random Variables

圖片
當產生訊號的過程過於複雜,以致無法建立合理的確定性模型(deterministic model)時,通常會使用隨機程序(random processes)來建立訊號模型。 10-6.FOURIER ANALYSIS OF STATIONARY RANDOM SIGNALS: THE PERIODOGRAM Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd ed https://www.brainkart.com/article/Random-Processes_6476/ Random:adj. 隨機的 在沒有特定方法或沒有經過有意識決策的情況下所做、完成或發生的。 在統計學(Statistics)中,指由隨機機率(chance)所支配或涉及相等機率的情況。 Stochastic:adj. 隨機的/隨機性的 具有隨機的機率分布(random probability distribution)或模式,可以透過統計方法進行分析,但可能無法被精確預測。 Process:n. 過程 為了達成某個特定目標而採取的一連串動作或步驟。 隨機過程  Stochastic(Random) Processes 為什麼要引入隨機過程(random process)? 為了方便對一個系統進行分析。 在一個系統中,例如通訊系統(communication system),可以考慮兩種模型: 確定性模型(Deterministic model) 對其隨時間變化的(精確)行為不存在不確定性。 隨機模型/隨機過程模型(Random or Stochastic model) 對其隨時間變化(time-dependent)的行為存在不確定性,但對其統計行為(statistical behavior)則具有一定的確定性。 隨機變數 Random Variables 在機率空間 ( S , F , P ) (S,\mathcal{F},P) 上的一個隨機變數,其中: F \mathcal{F} 是一個 σ-field, 對於每一個 x ∈ R x\in\mathbb{R} , { ζ ∈ S : x ( ζ ) ≤ x } 都必須是一個 事件(even...

DSP筆記3_Power Spectrum Estimation_功率譜及功率譜密度(Periodogram 與 Welch’s Method)

圖片
當訊號可以由多個正弦波的總和,或由週期性脈衝序列(periodic pulse train)激發的線性系統(linear system)來建模時,訊號有限長度分段(finite-length segments)的傅立葉轉換,可以很方便且自然地利用傅立葉轉換(Fourier transforms)、加窗(windowing)與線性系統理論(linear system theory)來加以解釋。 然而,對於較類似雜訊(noiselike)的訊號,通常會將其建模為隨機訊號(random signals),其。 日常生活中充斥的訊號,常見線上心電圖、語音、雷達等等,絕大部分都是非平穩訊號(nonstationary signals)。 截自 Fundamentals of Sound and Vibration 功率譜(Power Spectrum) 功率譜描述的訊號的功率或能量如何分布在各頻率上,是一種頻域中的表示方法。 可觀測出訊號是由那些頻率的訊號所組成以及它的能量頻寬、分布。 常應用在訊號處理、通訊、震動分析、雜訊描述。 比方我們有一個x(t)訊號 ,因為要將訊號從時間域轉換到頻率域。 經過傅立葉轉換後的式子可表示如下 常見轉換工具有傅立葉級數跟傅立葉轉換,如果是 固定規則(週期性訊號)的訊號 採用傅立葉級數, 不規則訊號(非週期性訊號) 就採用傅立葉轉換。則其Power Spectrum就可用對X(f)取絕對值平方表示,它表示各個頻率成分對訊號功率的貢獻程度。 訊號分類 根據訊號在「無限時間範圍內的總能量」與「長時間平均功率」來分類,可細分成兩種類: 能量訊號(Energy signals) 總能量是有限的,平均功率為0。 功率訊號(Power signals) 總能量是無限的,平均功率是有限的。 Power Spectral Density of Random Signals 通常在實際應用中所遇到的大多數訊號,其未來的變化無法被精確得知。 我們只能針對這些變化做出機率性的敘述。 用來描述這類訊號的數學工具,是由一組可能實現值(realizations)所構成的隨機序列(random sequence),其中每一個實現值都有其對應的發生機率。 然而,在所有可能的實現值中,實驗者通常只能觀察到訊號的一個實現值,因此無法直接照搬針對確定性訊號所建立的定義來套用。因為...

DSP筆記2_傅立葉轉換於頻譜分析的應用筆記

圖片
  關於「譜」 https://www.youtube.com/watch?v=F2Yl7N2uQWs 關於譜的觀念可以從光學領域追朔,以前學過當一束光經過三稜鏡,會根據光的波長(週期)被分解為紅橙黃綠藍....等色彩,也稱作(可見光)光譜。 對於訊號而言,頻譜分析也就是所謂借助傅立葉變化,類似三稜鏡一般將訊號轉換到頻率域,藉此去解析出構成原始訊號不同頻率成分。而描述頻率分量的曲線就稱作訊號的頻譜。 如圖中所示,原始波形訊號經由傅立葉分解出三條不同訊號波長的偕波集合。而此時若我們從這張圖的時間軸的朝向去看,波長的倒數也就是頻率,可去看不同頻率上的幅值。 傅立葉這位法國數學家,當時提出「任何時間函數都可分解為若干正弦訊號的總和」。 正弦函數是構成訊號的描述基本單元 以下這三個參數通常就能用於區分不同正弦訊號的Identifier,更可以理解為訊號的資訊載體。 A 振幅,反映訊號的強度(顯著程度) ω 頻率,衡量單位時間內訊號均勻重複的次數 Φ 相位,衡量時間延遲或時間偏移,可提供我們訊號於全局時間軸上出現時刻。 每個參數皆能乘載反映訊號產生原本質特徵的識別資訊 像是以下這個週期性波形,僅依靠肉眼觀察。我們無法很明確得知其數學模型。 無法得知該波形的準確頻率。若只依靠時域資訊,我們對複雜模式分析能力有限。 這也是為何需要將其分解為正弦訊號的原因。 因為我們人類對正弦訊號的特性比較熟悉。 傅立葉變換的核心是根據頻率,將原始訊號分解為若干正弦訊號的總和。 比方下面左側原始波形,經過傅立葉轉換分解後,得知可以用4個正弦訊號合成這個波形。 所以可以檢測到四個頻率成分 若要把這塊訊息用更簡潔呈現方式,可以透過繪製一個新的圖表。 橫軸以頻率而非時間作為自變數,縱軸以振幅來表示。 第一個正弦訊號在頻率對應位置是1,振幅為3。 第二個正弦訊號在頻率對應位置是7,振幅為1.5。 其他依此類推。 而相位參數則可以獨立呈現在另一張圖中 傅立葉函數以時域訊號作為輸入,輸出包含幅度和相位資訊的頻率函數。 當我們有能力拆解出頻譜資訊後可以做捨麼用途? 在工業領域中,可以想像遇到類似區分齒輪箱振動訊號原始波形,若從原始時域訊號難以發現問題,但從頻譜則可以清楚看到主要振動頻率,從而延伸做像是按傳動關係進一步推定振動源或故障異常分析等應用。 要從給定時域訊號中提取傅立葉變化就兩步驟 Step1....