update MVDR simulation

source/_posts/2024-02-01-algo-MVDR.md

@@ -82,13 +82,25 @@ $$\begin{equation}
 
 # MVDR 算法仿真
 
-## 理想空域谱
+  使用 ULA（Uniform Linear Array，均匀线阵）对 MVDR 算法进行仿真，阵元间距为半波长。
 
-![理想 Capon 谱](../images/post/2024-02-01-algo-MVDR/2024-02-01-algo-MVDR-010-IdealCaponSpectrum.svg){width=800px}
+## 理想情况
+
+  在理想情况下，各阵元的位置与预期一致，各阵元的复加权器是没有误差的，目标信号来向与期望方向一致，且信号和干扰加噪声可以分离。
+
+### 空域谱
+
+  Capon 空域谱表示 MVDR 复加权在整个空域范围内的响应。假设 1 个信号和 6 个干扰同时入射，信噪比 $\mathrm{SNR}=30\mathrm{dB}$，干噪比 $\mathrm{INR} = 70\mathrm{dB}$，信号和干扰来向在 $[0^\circ,180^\circ]$ 内服从均匀分布，使用 1000 个快拍进行仿真。
+
+  得到的空域谱如[图 1](#fig.1)，其中期望来向用红色点划线表示，干扰来向用绿色点划线表示。仿真同时还计算了输出 SINR，用于衡量残余干扰的大小。
+
+<a id="fig.1"></a>
+
+![图 1 理想 Capon 谱](../images/post/2024-02-01-algo-MVDR/2024-02-01-algo-MVDR-010-IdealCaponSpectrum.svg){width=700px}
 
 {% note primary MATLAB Code %}
 
-``` matlab mvdr_idea.m
+``` matlab mvdr_spec.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: <https://doi.org/10.1109/PROC.1969.7278
@@ -124,7 +136,7 @@ sigAng = rand(1,sigNum)*180;
 snr = [30,70];
 
 % sample number
-sampleNum = 9000;
+sampleNum = 1000;
 
 % parameter match check
 if (length(snr) ~= sigNum) | (length(sigAng) ~= sigNum)
@@ -254,6 +266,8 @@ for i = 1:sigNum
 end
 hold off
 
+legend({'Capon Spec','Desired','Interference'},Location="southeast")
+
 axis tight
 title('Ideal Capon Spectrum')
 xlabel('Degree/°')
@@ -264,11 +278,159 @@ grid minor
 
 {% endnote %}
 
-## 理想误码率
+### 分辨率
 
-## 失配输出 SINR
+&emsp;&emsp;MVDR 算法的分辨率用来衡量 MVDR 算法对接近目标方向的干扰的抑制能力。假设 1 个信号和 1 个干扰同时入射，信噪比 $\mathrm{SNR}=0 \mathrm{dB}$，干噪比 $\mathrm{INR} = 70\mathrm{dB}$，假定目标来向为 90°，干扰和信号来向差在 0°~50° 范围内扫描，对阵元数 $N=4,8,12,16$ 四种情况进行仿真。
 
-## 失配误码率
+&emsp;&emsp;[图 2](#fig.2) 以 MVDR 算法的输出 SINR 作为指标,给出了随期望、干扰信号来向的角度差的变化曲线，其中输出 SINR 归一化到 0 dB，便于不同阵元数情况的对比。点划线给出了低于理想 SINR 3 dB 的值，将点划线和实线的交点作为 MVDR 算法的分辨率。
+
+<a id="fig.2"></a>
+
+![图 2 Capon 分辨率](../images/post/2024-02-01-algo-MVDR/2024-02-01-algo-MVDR-020-IdealCaponResolution.svg){width=700px}
+
+{% note primary MATLAB Code %}
+
+``` matlab mvdr_resolution.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: <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
+
+rng(253)
+
+%% Scan Configuration
+elemNumArr = [4,8,12,16];
+intAngArr = 0.1:0.1:50;
+
+sinrOut = zeros(length(elemNumArr),length(intAngArr));
+
+for elemNumIdx = 1:length(elemNumArr)
+for intAngIdx = 1:length(intAngArr)
+
+%% Array Parameters
+
+% elements in array
+elemNum = elemNumArr(elemNumIdx);
+
+% spacing wrt wavelength
+elemSpac = 0.5;
+
+% element locations
+elemLoc = 0:elemSpac:(elemNum-1)*elemSpac;
+
+% Signal Parameters
+
+% signal number
+sigNum = 2;
+
+% direction of arrival
+desAng = 90;
+sigAng = [desAng,desAng+intAngArr(intAngIdx)];
+
+% signal to noise ratio
+snr = [0,70];
+
+% sample number
+sampleNum = 1024;
+
+%% Impinging Waves/Signal
+
+% 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
+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));
+
+% 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);
+
+%% AWGN Channel & Received Wave
+rxWave = rxSig + randn(elemNum,sampleNum);
+
+%% Array Processing
+
+% interference-plus-noise covariance matrix
+ipnCovar = rxIpn*rxIpn';
+
+% MVDR weight vector
+weight = ipnCovar\sigMF/(sigMF'*(ipnCovar\sigMF));
+%      = inv(ipnCovar)*sigMF/(sigMF'*inv(ipnCovar)*sigMF);
+
+% nomarlized output SINR
+sinrIdeal = snr(1)+pow2db(elemNum);
+sinrOut(elemNumIdx,intAngIdx) = pow2db(var(weight'*rxSig)/ ...
+  var(weight'*rxIpn)) - sinrIdeal;
+
+end
+
+plot(intAngArr,smoothdata(sinrOut(elemNumIdx,:)),DisplayName=['$N=', ...
+  num2str(elemNumArr(elemNumIdx)),'$'])
+hold on
+
+end
+
+%% figures
+plot([min(intAngArr),max(intAngArr)],[-3,-3],HandleVisibility="off")
+hold off
+
+title('MVDR Output SINR v.s. Angular Difference')
+xlabel('Angular Difference between Desigred and Interference Signal (°)')
+ylabel('Normalized Output SINR (dB)')
+ylim([-5 max(smoothdata(sinrOut,2),[],'all')+1])
+
+legend(Interpreter="latex",Location="southeast");
+
+grid on
+```
+
+{% endnote %}
+
+### 理想误码率
+
+## 失配情况
+
+### 分辨率
+
+### 误码率
 
 # Reference