用Python編程實作功率譜估計的平滑改進
周期圖法求解信號的功率譜
周期圖法求解信號功率譜的基本思想是：將隨機信號x(n)的N點觀測資料視為能量有限信號，對該信號求解傅里葉變換，求解變換結果幅值的平方，并且除以N，
這一部分的代碼如下：
#sig為所給信號，N為信號的長度
Pper=np.zeros(N)
xfft=np.fft.fft(sig)
for i in np.arange(N):
Pper[i]=np.abs(xfft[i])*np.abs(xfft[i])/N
試驗時所用的信號為：
N=200
sig=np.random.randn(N,)
程式運行結果為：
對周期圖法計算出的功率譜進行平滑改進
在本實驗中，利用了Welch法對周期圖法進行了改進，其基本思想是：對原始信號x(n)進行分段，并且允許每一段的資料有部分交疊；此外Welch法還允許選擇不同的窗函式來進行分段，本文利用了最簡單的窗函式來進行分段，
在利用Welch法時，首先要計算信號的分段數目，計算信號分段數的程式為：
def segL(N,M, overlap): #N：信號長度；M:分段長度；overlap:重疊長度；回傳值為分段數目sn，
d=M-overlap
n=np.int32(np.ceil((N-overlap)/d))
return n
在計算完分段數后，就可以構造一個由所有段信號組成的sn*M的xi矩陣，對該矩陣的每一行再進行加窗操作，窗函式的寬度為該矩陣的列數，
加窗函式的程式為：
for i in np.arange(sn):
xi=xi_snM[i] #xi_snM為求出的sn*M的xi矩陣
for j in np.arange(M):
xid[j]=xi[j]*d[j] #d為所用的窗函式
xifft[i]=np.fft.fft(xid) #計算加窗后信號的傅里葉變換
到此為止，計算出了snM矩陣的每一行信號的傅里葉變換，其中xifft也是一個snM的矩陣，利用Welch法計算功率譜時，需要對xifft進行變換，變換的目的主要是讓xifft的各行之間按照分段時的原則進行每段功率譜的相加求和，本文用的M=50，overlap=25，計算出的每小段開始時的索引號為：
計算完成的xifft和加窗后并且進行索引號變換的的xifftwin的資料如下圖所示：
原始信號和各小段信號組合后的新信號的圖形如下圖所示：
最后就算xifftwin每一列幅值的平方和，并且除以MUsn，就可以得到改進后的功率譜，其中U為歸一化因子，計算的程式為：
UWinwd=0
for n in np.arange(Winwd): #Winwd為窗函式的寬度，此處Winwd=M，
UWinwd=d[n]*d[n]+UWinwd #d為所選用的窗函式，本文選用的為矩形窗，計算出的歸一化因子為1，
U=UWinwd/Winwd
計算改進后的功率譜的程式為：
Pperwin=np.zeros(Nex)
for j in np.arange(Nex):
for i in np.arange(sn):
Pperwin[j]=np.abs(xifftwin[i,j])*np.abs(xifftwin[i,j])+Pperwin[j]
Pperwin[j]=Pperwin[j]/(M*U*sn)
最終得到的結果為：
如果本文有什么錯誤，還請讀者諒解，并且能夠幫忙指正，感激不盡！！！
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/181313.html
標籤：AI
上一篇：2020-10-18周總結