第25個冬天

  Photoplethysmography (PPG) 光體積變化描記法是一種用於監測多種生命徵象的光學技術，例如脈搏率、呼吸率及血氧飽和度，由於其非侵入性的特點而被廣受採用。 於早期研究中可得知，這些變化也可透過遠距方式測量，進而形成遠程光體積變化描記法（rPPG）的方法。於學術研究中，僅需低成本的RGB攝影機就能夠做實驗。而rPPG技術受限受測者不能有劇烈運動。當受測者只要出現顯著動作，現有演算法便會失效。 在2013年時候Gerard de Haan釋出一篇論文是探討「基於色度的 rPPG 強健脈搏率」， 提出新的 rPPG 技術比過去盲源分離方法(BSS-based methods)來得更穩健，其在訊雜比（SNR）與抗動作干擾能力方面，均優於所有先前的方法。他們不是單純從 RGB 中找週期，而是從「皮膚反射光的物理模型」推導出色度方法。論文的實驗包含 117 位受試者，並另外用腳踏車與踏步器測試動作干擾。 CHROM 到底想解決什麼問題？ 假設你對一張臉拍攝影片，每一幀都取臉部皮膚 ROI 的平均顏色，就會得到三條隨時間變動的訊號： \[ R(t),\quad G(t),\quad B(t) \] 心臟跳動時，皮膚中的血液量會微微改變，因此吸收、反射的光也會跟著改變。攝影機雖然看不到明顯的「臉變色」，但 RGB 數值其實會有非常微小的週期變化。 論文一開始將某個顏色通道 \(C\) 寫成： \[ C_i=I_i^C(\rho_{dc}^{C}+\rho_i^C) \] 其中 \(C\) 可以是 R、G、B。 可以理解成 :  攝影機看到的顏色=照到臉上的光x皮膚反射回來的顏色。 而皮膚反射又可以想成： 皮膚原本的顏色+心跳造成的微小變化 也就是說，我們真正要找到的是那個 非常小的心跳變化 。 但攝影機同時還會看到燈光亮暗、臉部移動、反光等大得多的變化。 第一步：先把 R、G、B 各自正規化 如果今天整張臉突然變亮，R、G、B 都可能一起變大。 因此作者先把每個顏色通道除以該時間區間的平均值： \[ R_n=\frac{R}{\mu(R)} \] \[ G_n=\frac{G}{\mu(G)} \] \[ B_n=\frac{B}{\mu(B)} \] 這樣就可以降低不同膚色、不同曝光亮度造成的絕對數值差異。 例如某一段影片： \[ \...

事實上，在上一篇我們完全不依賴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'...

功率頻譜密度（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} 意思是一秒鐘會完成...

  能量頻譜密度 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 ...

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

方法覆寫Override 程式1-原本子類別方法呼叫 class Parent : def myMethod ( self ): print ( "父類別方法" ) class Child (Parent): def myMethod ( self ): print ( "子類別方法" ) c = Child() c.myMethod() 程式2-方法覆寫Override super(子類別,self).同名方法 class Parent : def myMethod ( self ): print ( "父類別方法" ) class Child (Parent): def myMethod ( self ): print ( "子類別方法" ) super (Child, self ).myMethod() c = Child() c.myMethod() 方法多載Overloading 在python程式語言中，加減乘除四則運算子皆可用於預設變數、list的資料型別，但是如果是兩個自訂物件要進行運算，就勢必要重新定義這些運算子，不然就會發生不如預期的錯誤。 除了算術運算以外，比如要針對物件之間的比較運算、邏輯運算也都需要進行overloading。 固定覆寫寫法要留意，是以下這些制式化方法名稱 以下是一個簡單的二維點座標案例 TypeError: unsupported operand type(s) for +: 'Point' and 'Point' 另外像是要打印點座標直接print也會直接印出記憶體位置 程式1-As If class Point : def __init__ ( self ,x= 0 ,y= 0 ): self .x = x self .y = y p1 = Point( 4 , 3 ) p2 = Point( 2 , 1 ) #錯誤示範 #print(p2) #print(p1+p2) 程式2-To Be class Point : ...

  一個簡單的Employee物件示範 __init__用途是定義物件剛建立後，預計初始化的流程 類別內的method和一般python函數定義不同，必須包含參數self，並且務必放在第一個參數位置 在單一個main.py程式定義物件並使用物件 class Employee : empCount = 0 #定義物件剛建立後，預計初始化的流程。 #self就是物件自己本身 def __init__ ( self , name, salary): self .name = name self .yoursalary = salary Employee.empCount += 1 #和一般python函數定義不同，類別內的method必須包含參數self，並且務必放在第一個參數位置 def displayTotalEmpCount ( self ): print ( "Total Employee Count: %d" % self .empCount) def displayEmpInfo ( self ): print ( "Name:" , self .name, "Salary:" , self .yoursalary) emp1 = Employee( "Mike" , 42000 ) emp1.displayEmpInfo() emp1.displayTotalEmpCount() emp2 = Employee( "Michelle" , 50000 ) emp2.displayEmpInfo() emp1.displayTotalEmpCount() emp2.displayTotalEmpCount() 可看到empCount會是被各個獨立的物件實體，所共用的屬性，類似以前C#、Java的static修飾。 練習程式2-獨立物件實體成員屬性額外定義與刪除 class Employee : empCount = 0 #定義物件剛建立後，預計初始化的流程。 #self就是物...