連著三天8:00上課，我真的是痩不了了呀。
一、迭代法收斂理論
1.結論：對角占優(yōu)（對角線元素大）似乎比較有利于收斂，且GS法收斂的更快！
代碼驗證一下，用一個已知可以求解的例子，分別用雅克比和高斯賽德爾方法迭代求解。令其對角元素逐漸增大，繪制其迭代次數隨增大數值的變化曲線。
def Jocobi(A,b,initial,delta,isPrint=False):
'''
@author:zengwei
輸入：A是系數矩陣，N階方陣
b是N*1列向量
initial是解的初始值，N*1大小
delta是我給定的迭代值與真實值之間的誤差
輸出：迭代后的求解結果
'''
D = np.diag(np.diag(A)) # 獲得D矩陣
L = -np.tril(A,-1) # 獲得L矩陣
U = -np.triu(A,1) # 獲得U矩陣
try:
d = np.linalg.inv(D) # 對D矩陣求逆
BJ = np.dot(d,L+U) # 迭代矩陣BJ
lamda,_ = np.linalg.eig(BJ)
if np.max(np.abs(lamda))<1: # 譜半徑小于1
f = np.dot(d,b)
X = np.dot(BJ,initial)+f # 初次的解
k = -np.log(delta)/-np.log(np.max(np.abs(lamda)))
BJnorm = np.linalg.norm(BJ)
times = 1 # 因為前面有了一次迭代
while (BJnorm**k/(1-BJnorm))*np.linalg.norm(X - initial)> delta:
initial = X
X = np.dot(BJ,initial) + f
times = times +1
if isPrint==True:
print('譜半徑為：',np.max(np.abs(lamda)))
print('理論上的最大迭代次數為：%d次' %(int(k)+1))
print("實際上的最大迭代次數為：%d次" %times)
return X,times
else:
print('Sorry,不可收斂。')
print('譜半徑為：',np.max(np.abs(lamda)))
except:
print('對角矩陣D沒有逆矩陣！')
def Seidel(A,b,initial,delta,isPrint=False):
'''
@author:zengwei
輸入：A是系數矩陣，N階方陣
b是N*1列向量
initial是解的初始值，N*1大小
delta是我給定的迭代值與真實值之間的誤差
輸出：迭代后的求解結果
'''
D = np.diag(np.diag(A))
L = -np.tril(A,-1)
U = -np.triu(A,1)
try:
d = np.linalg.inv( D - L ) # 這里是和雅克比不同的地方
BG = np.dot(d,U) # 迭代矩陣BG
lamda,_ = np.linalg.eig(BG)
if np.max(np.abs(lamda))<1: # 譜半徑小于1
f = np.dot(d,b)
X = np.dot( BG ,initial ) + f
k = -np.log(delta)/-np.log(np.max(np.abs(lamda)))
BGnorm = np.linalg.norm(BG)
times = 1 # 因為前面有了一次迭代
while (BGnorm**k/(1-BGnorm))*np.linalg.norm( X - initial ) > delta:
initial = X
X = np.dot( BG,initial )+f
times = times +1
if isPrint==True:
print('譜半徑為：',np.max(np.abs(lamda)))
print('理論上的最大迭代次數為：%d次' %(int(k)+1))
print("實際上的最大迭代次數為：%d次" %times)
return X,times
else:
print('Sorry,不可收斂。')
print('譜半徑為：',np.max(np.abs(lamda)))
except:
print('矩陣[D-L]沒有逆矩陣！')
A = np.array([[8,-3,2],[4,11,-1],[6,3,12]],dtype=float)
b = np.array([[20],[33],[36]],dtype=float)
initial = np.zeros((3,1))
delta = 10**(-5)
Jocobi(A,b,initial,delta,isPrint=True)
'''
譜半徑為： 0.35924985028455664
理論上的最大迭代次數為：12次
實際上的最大迭代次數為：13次
array([[2.99999933],
[2.00000577],
[1.00000419]])
'''
Seidel(A,b,initial,delta,isPrint=True)
'''
譜半徑為： 0.13055824196677338
理論上的最大迭代次數為：6次
實際上的最大迭代次數為：7次
array([[3.00000201],
[1.9999987 ],
[0.99999932]])
'''
def addII(A,s):
A = A.copy()
for i in np.arange(len(A)):
A[i,i] = A[i,i] + s
return A
def getTimeslist(N,S,JorGs):
'''
增加N次，每次增加s大小,JorGs做算法選擇
'''
timeslist = []
n = np.arange(S,N*S+S,S)
for s in n:
newA = addII(A,s)
if JorGs=='Jocobi':
_,times = Jocobi(newA,b,initial,delta)
if JorGs=='Seidel':
_,times = Seidel(newA,b,initial,delta)
timeslist.append(times)
timeslist = np.array(timeslist)
return n,timeslist
N = 50
S = 5
n1,Jtimeslist = getTimeslist(N,S,'Jocobi')
n2,GStimeslist = getTimeslist(N,S,'Seidel')
plt.figure(figsize=(12,5))
plt.scatter(n1,Jtimeslist)
plt.plot(n1,Jtimeslist)
plt.scatter(n2,GStimeslist)
plt.plot(n2,GStimeslist)
plt.title('對角元素加%d的迭代次數變化圖'%S)
plt.legend(['Jocobi','Seidel'])
plt.xlabel('對角元素每次加%d'%S)
plt.ylabel('迭代次數')
plt.savefig('result.png')
plt.show()
迭代速度對比圖.png
從測試來看，結論得以驗證。
2.用迭代法計算希爾伯特矩陣
生成希爾伯特矩陣：
def Hilbert(N):
return 1. / (np.arange(1, N+1) + np.arange(0, N)[:, np.newaxis])
生成一個6階的希爾伯特矩陣，并且解均為1。
當我用雅克比迭代法進行計算時，會顯示不可收斂，因為譜半徑大于1了。當用高斯賽德爾迭代法時，會出現第一次迭代的結果。我注意到這時這個值會變得很大，且理論上的最大迭代次數
也變得很大。因為，寫SOR函數時先不用這些收斂條件，先粗糙地設置最大迭代次數，同時參考書后發(fā)現SOR對對稱正定的三對角矩陣能取到最優(yōu)的松弛因子
。一個不太成熟的SOR代碼：
def SOR(A,b,initial,delta,w,maxTimes,isw=False):
'''
@author:zengwei
輸入：A是系數矩陣，N階方陣
b是N*1列向量
initial是解的初始值，N*1大小
delta是我給定的迭代值與真實值之間的誤差
w:松弛因子
isw:對于對稱正定的三對角矩陣可以求最優(yōu)的w
輸出：迭代后的求解結果
'''
D = np.diag(np.diag(A))
L = -np.tril(A,-1)
U = -np.triu(A,1)
if isw==True:
BJ = np.dot(np.linalg.inv(D),L+U)
lamdaBJ,_ = np.linalg.eig(BJ)
rouBJ = np.max(np.abs(lamdaBJ))
w = 2./(1+np.sqrt(1-rouBJ**2))
try:
d = np.linalg.inv(D - w*L)
Bw = np.dot(d,(1-w)*D+w*U) # 迭代矩陣Bw
lamda,_ = np.linalg.eig(Bw)
print('譜半徑為：',np.max(np.abs(lamda)))
if np.max(np.abs(lamda))<1:
f = np.dot(d,w*b)
X = np.dot( Bw ,initial ) + f
for i in np.arange(maxTimes):
initial = X
X = np.dot( Bw,initial )+f
return X
else:
print('Sorry,不可收斂。')
print('譜半徑為：',np.max(np.abs(lamda)))
except:
print('矩陣[D-wL]沒有逆矩陣！')
我用SOR迭代法計算上述6階希爾伯特矩陣時可以得到結果，并分別試了松弛因子為1.0,1.25,1.5時得到的結果。我沒法用肉眼分辨哪個結果更好，因為我使用迭代值與真實值的均方誤差來進行比較。
N = 6
A = Hilbert(N)
x = np.ones((N,1))
b = np.dot(A,x)
initial = np.zeros((N,1))
delta = 10**(-5)
def MES(y,yTrue):
return np.sum([i**2 for i in (y-yTrue)])/len(y)
yTrue = 1
print('w=1時，數值解與真實值之間的均方誤差為:',MES(SOR(A,b,initial,delta,1,100),yTrue))
print('w=1.25時，數值解與真實值之間的均方誤差為:',MES(SOR(A,b,initial,delta,1.25,100),yTrue))
print('w=1.5時，數值解與真實值之間的均方誤差為:',MES(SOR(A,b,initial,delta,1.5,100),yTrue))
'''
w=1時，數值解與真實值之間的均方誤差為: 0.012836171242905989
w=1.25時，數值解與真實值之間的均方誤差為: 0.02127700958940544
w=1.5時，數值解與真實值之間的均方誤差為: 0.03818250349102156
'''
二、自相關
自相關在工程上有非常多的應用，典型應用是從隨機信號中找出周期信號。噪聲的自相關會在0時刻取到最大值，隨后衰減到0。
參考《數字信號處理》（胡廣書）編寫自相關函數有兩種思路：一種是按照定義計算，可以轉換為矩陣運算；另一種是先補N個0，再做FFT，然后對N分之一幅度平方做IFFT。代碼如下：
def myAutocorrelation(SigA):
'''
@author:zengwei
自相關函數,屬于有偏估計，對應matlab的xcorr(SigA,'biased')
https://ww2.mathworks.cn/help/matlab/ref/xcorr.html
'''
N = len(SigA)
SigB = np.zeros(2*N-1)
SigB[0:N] = SigA
L = len(SigB)
SigC = np.array(SigB.tolist()*(L+1))
Matrix = np.zeros((L,L))
start = N-1
for i in np.arange(L):
end = start + L # 還可以用np.roll()函數
Matrix[i,:] = SigC[start:end]
start = end - 1
return np.dot(Matrix,SigB)/N
def selfAutocorrelation(signal):
N = len(signal)
signal2N = np.hstack((signal,np.zeros(N)))
X = np.fft.fft(signal2N)
x = (np.abs(X)**2)/(N)
rm = np.fft.ifft(x)
return np.roll(rm,N-1)[0:-1]
得到的結果與matlab中的xcorr(SigA,'biased')結果一致，但看到更多情況下用xcorr(SigA,'unbiased')。matlab說明文檔詳細說明了其中的區(qū)別：https://ww2.mathworks.cn/help/matlab/ref/xcorr.html 。
Fs = 500
n = np.arange(0,1,1/Fs)
noise = np.random.randn(len(n))
x = 2*np.cos(2*np.pi*5*n)+noise
cor = myAutocorrelation(x)
plt.figure(figsize=(12,5))
plt.subplot(211)
plt.title('輸入信號')
plt.plot(n,x)
plt.subplot(212)
plt.plot(cor)
plt.title('自相關結果')
plt.savefig('自相關.png')
plt.show()
自相關示例.png
對于如何利用自相關結果找周期信號，可以參見matlab說明文檔進一步了解：https://ww2.mathworks.cn/help/signal/ug/find-periodicity-using-autocorrelation.html
三、自相關法求功率譜估計
用Python做功率譜估計感覺怪怪的，因為matlab有現成的函數做各種功率譜估計。我這里用信號自相關函數的FFT變換做為功率譜估計。如下所示：
# 生成原始輸入信號
fs = 1000
t = np.arange(0, 1, 1/fs)
f = 100
x = 2*np.cos(2*np.pi*f*t) + np.random.randn(t.size)
# 自相關法功率譜估計
num_fft = 1024
cor_x = myAutocorrelation(x)
cor_X = np.fft.fft(cor_x, num_fft)
P_BT = np.abs(cor_X)
plt.figure(figsize=(12, 8))
ax=plt.subplot(311)
ax.set_title('輸入信號')
plt.tight_layout()
plt.plot(x)
ax=plt.subplot(312)
ax.set_title('自相關法功率譜估計P_BT')
plt.plot(P_BT[:num_fft//2])
plt.tight_layout()
ax=plt.subplot(313)
ax.set_title('自相關法功率譜估計20*log10(P_BT)')
plt.plot(20*np.log10(P_BT[:num_fft//2]))
plt.tight_layout()
plt.savefig('P_BT.png')
plt.show()
自相關功率譜估計.png
接下來試試對自相關函數加窗函數，其余操作不變，并且窗函數加在自相關函數整個上面。然后看看有啥變化。
def Hamming(N):
return np.array([0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1)) for n in range(N)])
def Hanning(N):
return np.array([0.5 - 0.5 * np.cos(2 * np.pi * n / (N - 1)) for n in range(N)])
def Rect(N):
return np.ones(N)
自相關功率譜估計加窗.png
放在一副圖上比較.png
不加窗，加矩形窗，加漢明窗，加漢寧窗，其方差分別為：579.6754，579.6754，422.8173，413.1546。感官上，有一丟丟改善。應該比不上分段加窗強。
四、直接法功率譜估計的平滑改進
【10月12日批注】我這里分段時沒有考慮到頻率分辨率的問題，可以移步下一篇
依舊對上述相同的信號進行處理。
1.Bartlett法
def Bartlett(signal,L,num_fft):
M = len(signal)//L
signal = signal.reshape(L,M)
P_BT = []
for i in signal:
cor_I = np.fft.fft(i,num_fft)
P_BTi = np.abs(cor_I)**2/M
P_BT.append(P_BTi.tolist())
P_BT = np.sum(np.array(P_BT),axis=0)/L
return np.abs(P_BT)
BP_BT = Bartlett(x,100,1024)
plt.figure(figsize=(15,5))
plt.title('Bartlett平滑')
plt.plot(np.log10(BP_BT[:num_fft//2]))
plt.savefig('Bartlett.png')
plt.show()
Bartlett.png
2.Welch法
這里用的是漢寧窗。
def Welch(sig,M,num_fft):
L = int((len(sig)-M/2)//(M/2))
signal = np.zeros((L,M))
for m in np.arange(L):
signal[m,:] = x[int(m/2):int(m/2)+M]
d2 = Hanning(M)
U = np.sum(d2**2)/M
P_BT = []
for i in signal:
cor_I = np.fft.fft(i*d2,num_fft)
P_BTi = np.abs(cor_I)**2/(M*U)
P_BT.append(P_BTi.tolist())
P_BT = np.sum(np.array(P_BT),axis=0)/L
return np.abs(P_BT)
WP_BT = Welch(x,10,1024)
plt.figure(figsize=(15,5))
plt.title('Welch平滑')
plt.plot(np.log10(WP_BT[:num_fft//2]))
plt.savefig('Welch.png')
plt.show()
Welch.png
感覺可能有錯誤的地方。明天繼續(xù)......
