采用前后向空间平滑（FBSS）的DOA估计算法（python）_python doa估计

相干信号会导致信源协方差矩阵秩亏，采用空间平滑算法可以实现降维处理，但是会牺牲阵列的孔径。

这里先采用前后向空间平滑算法解相干，再使用MUSCI算法进行DOA估计。

需要注意的是，在空间平滑之后，新的协方差矩阵的大小是子阵大小，不是原来的协方差矩阵大小。

import numpy as np
import scipy.signal as ss
import scipy.linalg as LA
import matplotlib.pyplot as plt
 
def awgn(x, snr):
    spower = np.sum((np.abs(x) ** 2)) / x.size
    x = x + np.sqrt(spower / snr) * (np.random.randn(x.shape[0], x.shape[1]) + 1j * np.random.randn(x.shape[0], x.shape[1]))
    return x
 
derad = np.pi / 180
radeg = 180 / np.pi
d = np.arange(0,3.5,0.5)
theta = np.array([10,30,60]).reshape(1,-1)
n=500
snr = 10
iwave = theta.size
A = np.exp(-1j*2*np.pi*d.reshape(-1,1)@np.sin(theta*derad))
 
def coherent_signal_Rxx(d, theta, n, snr):
    A = np.exp(-1j * 2 * np.pi * d.reshape(-1, 1) @ np.sin(theta * derad))
    S = np.random.randn(iwave-1, n)
    S = np.concatenate((S[0,:].reshape(1,-1),S), axis=0)
    X = A @ S
    X = awgn(X, snr)
    Rxx = X @ (X.conj().T) / n
    return Rxx
 
def  FBSS(Rxx, sub_num, iwave):   ##subarray element number
    K = Rxx.shape[0]
    N = K-sub_num+1
    J = np.flip(np.eye(K),axis=0)
    crfb = (Rxx + J@(Rxx.T)@J)/2
    crs = np.zeros([sub_num,sub_num])
    for i in range(N):
        crs = crs + crfb[i:i+sub_num,i:i+sub_num]
        # crs = crs + Rxx[i:i + sub_num, i:i + sub_num] #FBB
 
    crs = crs/N
 
    D, EV = LA.eig(crs)
    index = np.argsort(D)  # ascending order index
    EN = EV[:, index][:, 0:sub_num - iwave]  #
 
    Angles = np.linspace(-np.pi / 2, np.pi / 2, 360)
    numAngles = Angles.size
    SP = np.empty(numAngles, dtype=complex)
 
    for i in range(numAngles):
        a = np.exp(-1j * 2 * np.pi * d.reshape(-1, 1)[0:sub_num] * np.sin(Angles[i]))
        SP[i] = ((a.conj().T @ a) / (a.conj().T @ EN @ EN.conj().T @ a))[0, 0]
 
    SP = np.abs(SP)
    SPmax = np.max(SP)
    SP = 10 * np.log10(SP / SPmax)
    x = Angles * radeg
 
    return x, SP
 
co_Rxx = coherent_signal_Rxx(d=d,theta=theta,n=n,snr=snr)
x, SP = FBSS(Rxx=co_Rxx,sub_num=6,iwave=iwave)
plt.plot(x, SP)
plt.show()