MVDR最小方差无畸形相应波束形成器 Python 实现

# -*- coding: utf-8 -*-
# author：半斤地瓜烧
# time: 2021年11月15日
import numpy as np
# 对12阵元均匀线阵进行仿真，频率f = 10 G H z f=10\rm{GHz}f=10GHz，期望信号角度为ϕ d = 1 0 ∘ \phi_d=10^\circϕ
import matplotlib.pyplot as plt

M = 12  #
c = 3e8  # 光速
f = 10e9  # 频率
_l = c / f  # 波长
d = _l / 2  # 阵元间隔 ，应该是半波振子（望指教）

angle = [10, 30]
phi = [i * np.pi / 180 for i in angle]  # 两个信号的 角度
Kt = 1  # 目标信号数目
snr = 20  # 信噪比
inr = 10  # 干噪比 干扰信号/噪声信号
angle_sample_num = 3601  # 角度采样数

angle_array = [i * d for i in range(M)]  # 10，30°之间的阵列角度

# 计算振元增益函数
A_gain = []
for i in phi:
    temp_gain = []
    for m in angle_array:
        gain = np.exp(1j * 2 * np.pi * m * np.sin(i) / _l)
        gain = np.round(gain, 6)
        temp_gain.append(gain)
    A_gain.append(temp_gain)
A_gain = np.matrix(A_gain).T  # 增益函数的维度为 （12， 2 ）

#  构造信号
object_signal = []  # 1 * 1024 目标信号
raw_signal = []  # 2 * 1024 总信号
signal_data_length = 1024
Kj = len(angle) - Kt  # 只有两个信号，一个是干扰信号，一个是期望信号，因此干扰信号个数是 总信号 - 干扰信号
for i in range(signal_data_length):
    s = np.sqrt(10 ** (snr / 10)) * np.exp(1j * 2 * np.pi * f * i / (2 * signal_data_length))  # 其中snr 为信噪比
    object_signal.append(np.round(s, 6))

raw_signal.append(object_signal)  # 将期望信号加在信号里面

# sqrt(10^(inr/10)/2) * (randn(1) + 1j*randn(1)) * exp(1j * 2 * pi * f * i/(2*Ts))
interfering_signal = []  # 定义干扰信号
for i in range(signal_data_length):
    s = np.sqrt(pow(10, (inr / 10) / 2)) * (np.random.randn(1)[0] + 1j * np.random.randn(1)[0]) * \
        np.exp(1j * 2 * np.pi * f * i / 2 * signal_data_length)  # 干扰信号
    interfering_signal.append(np.round(s, 6))
    # type(np.random.randn(1)) 是一个 numpy.ndarray 类型,取其中的数出来要 np.random.randn(1)[0]

raw_signal.append(interfering_signal)
raw_signal = np.array(raw_signal)  # 原信号的大小 (2,1024)
# object_signal = np.array(object_signal)  #


# 接收数据向量
noise_signal = (np.random.randn(M, signal_data_length) + 1j * np.random.randn(M, signal_data_length)) / np.sqrt(
    2)  # shape = (12,1024)
noise_signal = np.matrix(noise_signal)
# n = (randn(M,Ts)+1j*randn(M,Ts))/sqrt(2);

# x = A * s + n;
receive_signal = np.matrix(A_gain) * np.matrix(raw_signal) + noise_signal  # (12,1024)
# 接收信号 = 天线增益函数 * 原生信号（也携带噪声） + 噪声信号 ,注意这里用矩阵的计算,计算的维度关系为  (12,2).(2,1024) + (12,1024)

# MVDR算法
# J = A(:,Kt+1:K) * s(Kt+1:K,:);
J = A_gain[:, 1] * raw_signal[1, :]
Rnj = (1 / signal_data_length) * (J * J.T + noise_signal * noise_signal.T)
# Rnj = 1/Ts * (J*J'+n*n');
Rpriv = np.linalg.pinv(Rnj)  # 伪逆 pinv(X)=(X.T*X)^-1* X.T ，容易发现 pinv(X)X=Ipinv(X)X=I

# w = Rpriv*A(:,1)/(A(:,1)'*Rpriv*A(:,1));

weight = (Rpriv * A_gain[:, 0]) / (A_gain[:, 0].T * Rpriv * A_gain[:, 0])  # (12, 1)

theta = np.linspace(-np.pi, np.pi, angle_sample_num)
y = []
for i in theta:
    a = []
    for m in angle_array:  # angle_array 是阵列的角度位置
        a.append([np.exp(1j * 2 * np.pi * m * np.sin(i) / _l)])
    a = np.array(a)
    # print(weight.shape) # 这个详细计算不懂看看shape
    # print(weight.T.shape)
    # print(a.shape)
    y.append(weight.T * a)
y = np.array(y)
  # 这里就是正常的dB 计算公式
y = list(y.squeeze())

ydB_list = []  # 计算dB
max_abs = abs(max((y)))
for i in y:

    ydB_list.append(20 * np.log10(abs(i) / max_abs))

plt.plot(theta * 180 / np.pi, ydB_list)

plt.show()