diff --git a/source/_posts/2024-02-01-algo-MVDR.md b/source/_posts/2024-02-01-algo-MVDR.md new file mode 100644 index 0000000..6c3a944 --- /dev/null +++ b/source/_posts/2024-02-01-algo-MVDR.md @@ -0,0 +1,277 @@ +--- +title: MVDR 算法仿真 +mathjax: true +comments: true +copyright: true +toc: + enable: true + number: true + wrap: true + expand_all: true + max_depth: 4 +abbrlink: 86ade7bc +date: 2024-02-01 10:09:39 +tags: + - Beamforming + - 波束成形 + - MVDR + - 误码率分析 +categories: + - - 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 谱](../images/post/2024-02-01-algo-MVDR/2024-02-01-algo-MVDR-010-IdealCaponSpectrum.svg){width=800px} + +{% note primary MATLAB Code %} + +``` matlab mvdr_idea.m +%% 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: . +% +% 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](https://doi.org/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. diff --git a/source/images/post/2024-02-01-algo-MVDR/2024-02-01-algo-MVDR-010-IdealCaponSpectrum.svg b/source/images/post/2024-02-01-algo-MVDR/2024-02-01-algo-MVDR-010-IdealCaponSpectrum.svg new file mode 100644 index 0000000..a6dca1b --- /dev/null +++ b/source/images/post/2024-02-01-algo-MVDR/2024-02-01-algo-MVDR-010-IdealCaponSpectrum.svg @@ -0,0 +1,264 @@ + + +020406080100120140160180Degree/°-60-50-40-30-20-100Magnitude/dBIdeal Capon Spectrum