Files
josh-blog-source/source/_posts/2024-02-01-algo-MVDR.md
T
2024-02-01 17:14:17 +08:00

7.8 KiB
Raw Blame History

title, mathjax, comments, copyright, toc, abbrlink, date, tags, categories
title mathjax comments copyright toc abbrlink date tags categories
MVDR 算法仿真 true true true
enable number wrap expand_all max_depth
true true true true 4
86ade7bc 2024-02-01 10:09:39
Beamforming
波束成形
MVDR
误码率分析
Josh 的学习笔记
算法仿真

  MVDR (Minimum Variance Distortionless Response,最小方差无失真响应)算法是 Capon 于 1969 年提出的经典波束形成算法,其解决的是无失真约束下输出噪声加干扰功率最小化问题。

阵列模型

  首先对阵列接收信号进行建模。假设 N 元阵列天线在 t 时刻的阵列接收信号和阵列输出信号分别为 \mathbf{x}(t) 和 $y(t)$,则对于 K 个入射信号,有

$$\begin{equation} \begin{cases} \mathbf{x}(t) = \mathbf{A}(\boldsymbol{\theta})\mathbf{s}(t) + \mathbf{n}(t) \ y(t) = \mathbf{w}^\mathrm{H}\mathbf{x}(t) \end{cases},\quad t = 1,\cdots,T \end{equation}$$

其中

  • 阵列接收信号 \mathbf{x}(t) = [x_1(t),\cdots,x_N(t)]^\mathrm{T}\in\mathbb{C}^{N\times 1}

  • 噪声 \mathbf{n}(t) = [n_1(t),\cdots,n_N(t)]^\mathrm{T}\in\mathbb{C}^{N\times 1}

  • 入射信号 \mathbf{s}(t) = [s_1(t),\cdots,s_K(t)]^\mathrm{T}\in\mathbb{C}^{K\times 1}

  • 阵列流形矩阵 $\mathbf{A}(\boldsymbol{\theta}) = [\mathbf{a}(\theta_1),\cdots,\mathbf{a}(\theta_K)]\in \mathbb{C}^{N\times K}$,导向矢量 $\mathbf{a}(\theta_i)\in\mathbb{C}^{N\times 1}$,K 个信号来向集合 \boldsymbol{\theta} = \{\theta_i\}_{i=1}^K

  • 波束形成复加权 \mathbf{w} = [w_1,\cdots,w_N]^\mathrm{T}

MVDR 算法简述

  假设第 1 个信号是期望信号(desired signal),令其导向矢量 $\mathbf{a}_\mathrm{d} = \mathbf{a}(\theta_1)$,则阵列接收信号可以进一步写为

$$\begin{equation} \mathbf{x}(t) = \underset{\mathbf{x}\mathrm{d}(t)}{\underbrace{\mathbf{a}\mathrm{d}s_1(t)}} + \underset{\mathbf{x}\mathrm{i+n}(t)}{\underbrace{\sum{k=2}^{K}\mathbf{a}(\theta_k)s_k(t) + \mathbf{n}(t)}} \end{equation}$$

至此,可以给出 MVDR 问题的基本表述

$$\begin{equation} \begin{aligned} &\underset{\mathbf{w}}{\text{minimize}} &&\mathbf{w}^\mathrm{H}\mathbf{R}{\mathrm{i+n}}\mathbf{w} \ &\text{subject to} && \mathbf{w}^\mathrm{H}\mathbf{a}\mathrm{d} = 1 \end{aligned} \end{equation}$$

其闭式解可以由 Lagrange 乘子法得到

$$\begin{equation} \mathbf{w} = \alpha\mathbf{R}{\mathrm{i+n}}^{-1}\mathbf{a}\mathrm{d} \end{equation}$$

其中

  • 噪声加干扰协方差矩阵 \mathbf{R}_{\mathrm{i+n}} = \displaystyle\frac{1}{T} \sum_{t=1}^{T}\mathbf{x}_\mathrm{i+n}(t)\mathbf{x}_\mathrm{i+n}^\mathrm{H}(t)

  • 归一化因子(使得复加权对于期望信号的增益为 1) \alpha = \displaystyle\frac{1}{\mathbf{a}_\mathrm{d}^\mathrm{H}\mathbf{R}_{\mathrm{i+n}}^{-1}\mathbf{a}_\mathrm{d}}

MVDR 算法仿真

理想空域谱

理想 Capon 谱{width=800px}

{% note primary MATLAB Code %}

%% High-Resolution Frequency-Wavenumber Spectrum Analysis
% J. Capon, “High-Resolution Frequency-Wavenumber Spectrum Analysis,” _Proceedings
% of the IEEE_, vol. 57, no. 8, pp. 1408–1418, Aug. 1969, doi: <https://doi.org/10.1109/PROC.1969.7278
% 10.1109/PROC.1969.7278>.
%
% H. L. Van Trees, _Optimum Array Processing: Part IV of Detection, Estimation,
% and Modulation Theory_. New York, NY: Wiley, 2002.

clear
close all
%% Parameters

rng(253)
% Array Parameters

% elements in array
elemNum = 12;

% spacing wrt wavelength
elemSpac = 0.5;

% element locations
elemLoc = 0:elemSpac:(elemNum-1)*elemSpac;
% Signal Parameters

% signal number
sigNum = 6;

% direction of arrival
sigAng = rand(1,sigNum)*180;

% signal to noise ratio
snr = [30,70];

% sample number
sampleNum = 9000;

% parameter match check
if (length(snr) ~= sigNum) | (length(sigAng) ~= sigNum)
  disp('WARN: Signal number and parameter number NOT match.')
end
%% Impinging Waves/Signals
% Transmitted Signals

% wave generation, here are some sets of uncorrelated white noise
txWave = randn(sigNum,sampleNum);

% transmit power control based on given snr
if length(snr) == 1
  txWave = db2mag(snr)*txWave;
elseif length(snr) >= sigNum
  if length(snr) > sigNum
    disp('WARN: Extra SNRs are ignored.')
  end
  for i = 1:sigNum
    txWave(i,:) = db2mag(snr(i))*txWave(i,:);
  end
elseif length(snr) < sigNum
  disp('WARN: The last SNR is used for the rest of the signals.')
  for i = 1:sigNum
    if i <= length(snr)
      txWave(i,:) = db2mag(snr(i))*txWave(i,:);
    else
      txWave(i,:) = db2mag(snr(end))*txWave(i,:);
    end
  end
end
% Array Manifold Matrix
% For each signal, the array manifold vector is
%
% $$\mathbf{a}\left(\theta \right)=\exp \left(-j\frac{2\pi \mathbf{d}}{\lambda
% \;}\cos \theta \;\right)$$
%
% in which the element locations are
%
% $$\mathbf{d}={\left\lbrack 0,\frac{\lambda }{2},\cdots ,\left(N-1\right)\frac{\lambda
% }{2}\right\rbrack }^{\mathrm{T}}$$

if length(sigAng) > sigNum
  disp('WARN: Extra angles are ignored.')
elseif length(sigAng) < sigNum
  disp('ERROR: Angle number is less than signal number.')
  return
end

arrayMF = exp(-1i*pi*elemLoc'*cos(sigAng(1:sigNum)/180*pi));

if sigNum > 1
  % signal & interference manifold vector
  sigMF = arrayMF(:,1);
  intMF = arrayMF(:,2:end);

  % impinging signal & interference
  rxSig = sigMF*txWave(1,:);
  rxIpn = intMF*txWave(2:end,:) + randn(elemNum,sampleNum);
else
  % signal manifold vector
  sigMF = arrayMF;

  % impinging signal
  rxSig = sigMF*txWave;
  rxIpn = randn(elemNum,sampleNum);
end
% Received Wave

rxWave = rxSig + randn(elemNum,sampleNum);
%% Array Processing
% The MVDR problem is to minimize output interference-plus-noise power under
% distortionless contraint of target signal
%
% $$\begin{array}{l}\underset{\mathbf{w}}{\mathrm{minimize}\;\;} \;{\mathbf{w}}^{\mathrm{H}}
% {\mathbf{R}}_{\mathrm{i}+\mathrm{n}} \mathbf{w}\\\mathrm{subecto}\;\mathrm{to}\;\;{\mathbf{w}}^{\mathrm{H}}
% \mathbf{a}=1\end{array}$$
%
% The problem can be solve using the Lagrange multiplier method, and the solution
% is
%
% $$\mathbf{w}=\alpha {\mathbf{R}}_{\mathrm{i}+\mathrm{n}}^{-1} \mathbf{a}$$
%
% where $\alpha \;$is the normalization factor that ensures the power of array
% input and output desired signal remains unchanged
%
% $$\alpha =\frac{1}{{\mathbf{a}}^{\mathrm{H}} {\mathbf{R}}_{\mathrm{i}+\mathrm{n}}
% \mathbf{a}}$$

% interference-plus-noise covariance matrix
ipnCovar = rxIpn*rxIpn';

% MVDR weight vector
weight = ipnCovar\sigMF/(sigMF'*(ipnCovar\sigMF));
%      = inv(ipnCovar)*sigMF/(sigMF'*inv(ipnCovar)*sigMF);

% output SINR
sinrOut = pow2db(var(weight'*rxSig)/var(weight'*rxIpn));
%% Spectrum & Figure

% scan grid
scanGrid = 0:0.1:180;

% spectrum
spec = zeros(size(scanGrid));

for i = 1:length(scanGrid)
  % steering vector
  scanMF = exp(-1i*pi*elemLoc'*cos(scanGrid(i)/180*pi));

  % MVDR spectrum
  spec(i) = weight'*scanMF;
end

% normalized spectrum
specNorm = pow2db(abs(spec));

% spectrum & signal location indicator
plot(scanGrid,specNorm)
hold on
for i = 1:sigNum
  if i == 1
    plot([sigAng(i),sigAng(i)],[max(specNorm),min(specNorm)],'-.r');
  else
    plot([sigAng(i),sigAng(i)],[max(specNorm),min(specNorm)],'-.g');
  end
end
hold off

axis tight
title('Ideal Capon Spectrum')
xlabel('Degree/°')
ylabel('Magnitude/dB')
grid on
grid minor

{% endnote %}

理想误码率

失配输出 SINR

失配误码率

Reference

  1. J. Capon, “High-Resolution Frequency-Wavenumber Spectrum Analysis,” Proceedings of the IEEE, vol. 57, no. 8, pp. 1408–1418, Aug. 1969, doi: 10.1109/PROC.1969.7278.

  2. H. L. Van Trees, Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory. New York, NY: Wiley, 2002.