連著三天8:00上課，我真的是痩不了了呀村刨。
一告抄、迭代法收斂理論
1.結(jié)論：對(duì)角占優(yōu)（對(duì)角線元素大）似乎比較有利于收斂，且GS法收斂的更快嵌牺！
代碼驗(yàn)證一下打洼，用一個(gè)已知可以求解的例子，分別用雅克比和高斯賽德?tīng)柗椒ǖ蠼饽娲狻Ａ钇鋵?duì)角元素逐漸增大募疮，繪制其迭代次數(shù)隨增大數(shù)值的變化曲線。
def Jocobi(A,b,initial,delta,isPrint=False):
'''
@author:zengwei
輸入：A是系數(shù)矩陣僻弹，N階方陣
b是N*1列向量
initial是解的初始值阿浓，N*1大小
delta是我給定的迭代值與真實(shí)值之間的誤差
輸出：迭代后的求解結(jié)果
'''
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) # 對(duì)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 # 因?yàn)榍懊嬗辛艘淮蔚? 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('理論上的最大迭代次數(shù)為：%d次' %(int(k)+1))
print("實(shí)際上的最大迭代次數(shù)為：%d次" %times)
return X,times
else:
print('Sorry,不可收斂。')
print('譜半徑為：',np.max(np.abs(lamda)))
except:
print('對(duì)角矩陣D沒(méi)有逆矩陣蹋绽！')
def Seidel(A,b,initial,delta,isPrint=False):
'''
@author:zengwei
輸入：A是系數(shù)矩陣搔扁，N階方陣
b是N*1列向量
initial是解的初始值，N*1大小
delta是我給定的迭代值與真實(shí)值之間的誤差
輸出：迭代后的求解結(jié)果
'''
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 # 因?yàn)榍懊嬗辛艘淮蔚? 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('理論上的最大迭代次數(shù)為：%d次' %(int(k)+1))
print("實(shí)際上的最大迭代次數(shù)為：%d次" %times)
return X,times
else:
print('Sorry,不可收斂蟋字。')
print('譜半徑為：',np.max(np.abs(lamda)))
except:
print('矩陣[D-L]沒(méi)有逆矩陣稿蹲！')
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
理論上的最大迭代次數(shù)為：12次
實(shí)際上的最大迭代次數(shù)為：13次
array([[2.99999933],
[2.00000577],
[1.00000419]])
'''
Seidel(A,b,initial,delta,isPrint=True)
'''
譜半徑為： 0.13055824196677338
理論上的最大迭代次數(shù)為：6次
實(shí)際上的最大迭代次數(shù)為：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('對(duì)角元素加%d的迭代次數(shù)變化圖'%S)
plt.legend(['Jocobi','Seidel'])
plt.xlabel('對(duì)角元素每次加%d'%S)
plt.ylabel('迭代次數(shù)')
plt.savefig('result.png')
plt.show()
迭代速度對(duì)比圖.png
從測(cè)試來(lái)看鹊奖，結(jié)論得以驗(yàn)證苛聘。
2.用迭代法計(jì)算希爾伯特矩陣
生成希爾伯特矩陣：
def Hilbert(N):
return 1. / (np.arange(1, N+1) + np.arange(0, N)[:, np.newaxis])
生成一個(gè)6階的希爾伯特矩陣，并且解均為1忠聚。
當(dāng)我用雅克比迭代法進(jìn)行計(jì)算時(shí)设哗，會(huì)顯示不可收斂，因?yàn)樽V半徑大于1了两蟀。當(dāng)用高斯賽德?tīng)柕〞r(shí)网梢，會(huì)出現(xiàn)第一次迭代的結(jié)果。我注意到這時(shí)這個(gè)值會(huì)變得很大赂毯，且理論上的最大迭代次數(shù)
也變得很大战虏。因?yàn)榧鹪祝瑢?xiě)SOR函數(shù)時(shí)先不用這些收斂條件，先粗糙地設(shè)置最大迭代次數(shù)烦感，同時(shí)參考書(shū)后發(fā)現(xiàn)SOR對(duì)對(duì)稱正定的三對(duì)角矩陣能取到最優(yōu)的松弛因子
巡社。一個(gè)不太成熟的SOR代碼：
def SOR(A,b,initial,delta,w,maxTimes,isw=False):
'''
@author:zengwei
輸入：A是系數(shù)矩陣，N階方陣
b是N*1列向量
initial是解的初始值手趣，N*1大小
delta是我給定的迭代值與真實(shí)值之間的誤差
w:松弛因子
isw:對(duì)于對(duì)稱正定的三對(duì)角矩陣可以求最優(yōu)的w
輸出：迭代后的求解結(jié)果
'''
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]沒(méi)有逆矩陣！')
我用SOR迭代法計(jì)算上述6階希爾伯特矩陣時(shí)可以得到結(jié)果绿渣，并分別試了松弛因子為1.0,1.25,1.5時(shí)得到的結(jié)果朝群。我沒(méi)法用肉眼分辨哪個(gè)結(jié)果更好，因?yàn)槲沂褂玫蹬c真實(shí)值的均方誤差來(lái)進(jìn)行比較中符。
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時(shí)潜圃，數(shù)值解與真實(shí)值之間的均方誤差為:',MES(SOR(A,b,initial,delta,1,100),yTrue))
print('w=1.25時(shí)，數(shù)值解與真實(shí)值之間的均方誤差為:',MES(SOR(A,b,initial,delta,1.25,100),yTrue))
print('w=1.5時(shí)舟茶，數(shù)值解與真實(shí)值之間的均方誤差為:',MES(SOR(A,b,initial,delta,1.5,100),yTrue))
'''
w=1時(shí)，數(shù)值解與真實(shí)值之間的均方誤差為: 0.012836171242905989
w=1.25時(shí)堵第，數(shù)值解與真實(shí)值之間的均方誤差為: 0.02127700958940544
w=1.5時(shí)吧凉，數(shù)值解與真實(shí)值之間的均方誤差為: 0.03818250349102156
'''
二、自相關(guān)
自相關(guān)在工程上有非常多的應(yīng)用踏志，典型應(yīng)用是從隨機(jī)信號(hào)中找出周期信號(hào)阀捅。噪聲的自相關(guān)會(huì)在0時(shí)刻取到最大值，隨后衰減到0针余。
參考《數(shù)字信號(hào)處理》（胡廣書(shū)）編寫(xiě)自相關(guān)函數(shù)有兩種思路：一種是按照定義計(jì)算饲鄙，可以轉(zhuǎn)換為矩陣運(yùn)算；另一種是先補(bǔ)N個(gè)0圆雁，再做FFT忍级，然后對(duì)N分之一幅度平方做IFFT。代碼如下：
def myAutocorrelation(SigA):
'''
@author:zengwei
自相關(guān)函數(shù),屬于有偏估計(jì)伪朽，對(duì)應(yīng)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()函數(shù)
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]
得到的結(jié)果與matlab中的xcorr(SigA,'biased')結(jié)果一致轴咱，但看到更多情況下用xcorr(SigA,'unbiased')。matlab說(shuō)明文檔詳細(xì)說(shuō)明了其中的區(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('輸入信號(hào)')
plt.plot(n,x)
plt.subplot(212)
plt.plot(cor)
plt.title('自相關(guān)結(jié)果')
plt.savefig('自相關(guān).png')
plt.show()
自相關(guān)示例.png
對(duì)于如何利用自相關(guān)結(jié)果找周期信號(hào)伶授，可以參見(jiàn)matlab說(shuō)明文檔進(jìn)一步了解：https://ww2.mathworks.cn/help/signal/ug/find-periodicity-using-autocorrelation.html
三谁鳍、自相關(guān)法求功率譜估計(jì)
用Python做功率譜估計(jì)感覺(jué)怪怪的，因?yàn)閙atlab有現(xiàn)成的函數(shù)做各種功率譜估計(jì)。我這里用信號(hào)自相關(guān)函數(shù)的FFT變換做為功率譜估計(jì)阵苇。如下所示：
# 生成原始輸入信號(hào)
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)
# 自相關(guān)法功率譜估計(jì)
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('輸入信號(hào)')
plt.tight_layout()
plt.plot(x)
ax=plt.subplot(312)
ax.set_title('自相關(guān)法功率譜估計(jì)P_BT')
plt.plot(P_BT[:num_fft//2])
plt.tight_layout()
ax=plt.subplot(313)
ax.set_title('自相關(guān)法功率譜估計(jì)20*log10(P_BT)')
plt.plot(20*np.log10(P_BT[:num_fft//2]))
plt.tight_layout()
plt.savefig('P_BT.png')
plt.show()
自相關(guān)功率譜估計(jì).png
接下來(lái)試試對(duì)自相關(guān)函數(shù)加窗函數(shù)，其余操作不變灯抛，并且窗函數(shù)加在自相關(guān)函數(shù)整個(gè)上面妇押。然后看看有啥變化需了。
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)
自相關(guān)功率譜估計(jì)加窗.png
放在一副圖上比較.png
不加窗，加矩形窗橡疼，加漢明窗援所，加漢寧窗，其方差分別為：579.6754欣除，579.6754住拭，422.8173，413.1546历帚。感官上滔岳，有一丟丟改善。應(yīng)該比不上分段加窗強(qiáng)挽牢。
四谱煤、直接法功率譜估計(jì)的平滑改進(jìn)
【10月12日批注】我這里分段時(shí)沒(méi)有考慮到頻率分辨率的問(wèn)題，可以移步下一篇
依舊對(duì)上述相同的信號(hào)進(jìn)行處理禽拔。
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
感覺(jué)可能有錯(cuò)誤的地方。明天繼續(xù)......
