josh/josh-blog-source
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

update several posts, add matlab code to figs

Browse Source
This commit is contained in:
josh
2023-06-20 19:17:57 +08:00
parent 67b99a2560
commit 85e2b1bb85
5 changed files with 727 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files
_config.next.yml
+1 -1
View File
@@ -531,7 +531,7 @@ pdf:
# Mermaid tag
mermaid:
enable: false
enable: true
# Available themes: default | dark | forest | neutral
theme:
light: forest
_config.yml
+2
View File
@@ -55,6 +55,8 @@ highlight:
tab_replace: ' '
wrap: true
hljs: false
exclude_languages:
- mermaid
prismjs:
enable: false
preprocess: true
source/_posts/2023-04-20-josh-oap-part-1-3.md
+292
View File
@@ -75,6 +75,51 @@ $$\begin{equation}
![图 1-3-2 用 dB 表示 $\left|\boldsymbol{\varUpsilon}(\psi)\right|$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-20-josh-oap-part-1-3/2023-04-20-josh-oap-part-1-3-020-FrequencyWavenumberResponseFunctionInDB.svg){width=1000px}
{% note primary MATLAB Code %}
``` matlab fig2_1516.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.15 & 2.16
% Beampattern of uniformly weighted linear array
% Xiaomin Lu
% Updated 1/5/99
% Last updated by K. Bell 7/22/01, 10/4/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N = 11; % 11 elements
% Note, Eq. (2.92) is in fact the standard Dirichlet function
% beam = diric(psi, N);
% if you don't have the signal processing toolbox, this does the same thing
psi = (-5:1/400:5)*pi;
beam = zeros(1,length(psi(:)));
y=sin(.5*psi);
i=find(abs(y)>1e-12); % set where x is not divisible by 2 pi
j=1:length(psi(:));
j(i)=[]; % complement set
beam(i)=sin((N/2)*psi(i))./(N*y(i));
beam(j)=sign(cos(psi(j)*((N+1)/2)));
plot(psi/pi, beam)
grid
xlabel('$\psi/\pi$','Interpreter','latex')
ylabel('Frequency wavenumber response function')
axis([-5 5 -0.4 1])
set(gcf,'Position',[0 0 1000 600])
figure
beam = 20*log10(abs(beam));
plot(psi/pi,beam);
xlabel('$\psi/\pi$','Interpreter','latex')
ylabel('Frequency wavenumber response function (dB)')
axis([-5 5 -25 0])
grid
set(gcf,'Position',[0 0 1000 600])
```
{% endnote %}
若不指定 $\boldsymbol{w}$ 的类型，则 $\boldsymbol{\varUpsilon}_\psi(\psi)$ 是复数，所以相位也应该画出来。但是，均匀加权线阵具有对称性，因此得到的频率-波数响应是实函数。
  也可以用 $k_z$ 来表示频率-波数响应
@@ -121,10 +166,138 @@ $$\begin{equation}
![图 1-3-3 在极坐标系中画出 $B_\theta(\theta)$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-20-josh-oap-part-1-3/2023-04-20-josh-oap-part-1-3-030-PloarPlotOfBeamPatternInAngleSpace.svg){width=800px}
{% note primary MATLAB Code %}
``` matlab fig2_17.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.17
% Polar plot of B (Theta)
% Lillian Xu 1/5/99
% Updated by K. Bell 6/25/01
% Functions called: polardb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N = 11; % Elements in array
d = 0.5; % spacing wrt wavelength
beamwidth = 2/(N*d); % null-to-null (LL BW is half this)
D=d*(-(N-1)/2:1:(N-1)/2); % element locations
ang = pi*(-1:0.001:1);
u = cos(ang);
n2=size(u,2);
AS = exp(-1i*2*pi*cos(90/180*pi)*D'); % BP points to 90 (broadside)
Au = exp(-1i*2*pi*D'*u);
B = real(AS'*Au)/N;
G = 20*log10((abs(B)));
h=polardb(ang,G,-40);
```
{% endnote %}
<a id="fig.1-3-4"></a>
![图 1-3-4 $d = \lambda/2, N=10$，线阵的 $\left|\boldsymbol{\varUpsilon}_\psi(\psi)\right|$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-20-josh-oap-part-1-3/2023-04-20-josh-oap-part-1-3-040-AbsOfFWRFForALinearArray.svg){width=1000px}
{% note primary MATLAB Code %}
``` matlab fig2_18.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.18
% beam patterns in different spaces
% Lillian Xiaolan Xu
% Last updated 09/07/2000
% updated by K. Bell 7/22/01, 10/4/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N = 10;
n = (-(N-1)/2:(N-1)/2).';
psi = pi*(-3:0.001:3);
w = 1/N*[ones(N,1)];
d = 1/2;
vv = exp(1i*2*d*n*psi);
B = (abs(w'*vv));
u1=-1;
u2=3;
y1=0;
y2=1;
subplot(4,1,1)
plot(-psi/d,B);
hold on
axis([u1*pi/d u2*pi/d y1 y2])
set(gca,'XTick',[-pi/d 0 pi/d 2*pi/d 3*pi/d])
set(gca,'YTick',[0 0.5 1])
set(gca,'XTickLabel',{'$-\pi/d$';'$0$';'$\pi/d$';'$2\pi/d$';'$3\pi/d$'}, ...
'TickLabelInterpreter','latex')
grid on
point=0.8;
point_up=0.84;
point_down=0.76;
plot([-pi/d pi/d],[point point],'--')
plot([-1*pi/d -0.9*pi/d],[point point_up])
plot([-1*pi/d -0.9*pi/d],[point point_down])
plot([0.9*pi/d pi/d],[point_up point])
plot([0.9*pi/d pi/d],[point_down point])
text(-0.3*pi/d,1.2,'Visible region')
plot([pi/d 3*pi/d],[point point],'--')
plot([pi/d 1.1*pi/d],[point point_up])
plot([pi/d 1.1*pi/d],[point point_down])
plot([2.9*pi/d 3*pi/d],[point_up point])
plot([2.9*pi/d 3*pi/d],[point_down point])
text(1.7*pi/d,1.2,'Virtual region')
legend({'$k_z$-space'},'Interpreter','latex')
subplot(4,1,2)
plot(psi,B);
hold on
axis([u1*pi u2*pi y1 y2])
set(gca,'XTick',[-pi 0 pi 2*pi 3*pi])
set(gca,'YTick',[0 0.5 1])
set(gca,'XTickLabel',{'$-\pi$';'$0$';'$\pi$';'$2\pi$';'$3\pi$'}, ...
'TickLabelInterpreter','latex')
grid on
legend({'$\psi$-space'},'Interpreter','latex')
subplot(4,1,3)
plot(psi/pi,B);
hold on
axis([u1 u2 y1 y2])
set(gca,'XTick',[-1 0 1 2 3])
set(gca,'YTick',[0 0.5 1])
grid on
legend({'$u$-space'},'Interpreter','latex')
theta = -180:0.1:360;
vv = exp(1i*2*d*n*pi*cos(theta/180*pi));
B = (abs(w'*vv));
subplot(4,1,4)
plot(theta,B);
hold on
axis([-180 180 y1 y2])
set(gca,'XTick',[-180 -90 0 90 180])
set(gca,'YTick',[0 0.5 1])
set(gca,'XTickLabel',{'$180^\circ$';'$90^\circ$';'$0^\circ$';'$-90^\circ$'; ...
'$-180^\circ$'},'TickLabelInterpreter','latex')
grid on
legend({'$\theta$-space'},'Interpreter','latex')
set(gcf,'Position',[0 0 800 800])
```
{% endnote %}
{% note info %}
虽然以上例子比较简单，但是可以用来展示线性阵列的一些重要特征。
@@ -151,6 +324,76 @@ $$\begin{equation}
![图 1-3-5 波束方向图的主波束](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-20-josh-oap-part-1-3/2023-04-20-josh-oap-part-1-3-050-MainLobeOfBeamPattern.svg){width=700px}
{% note primary MATLAB Code %}
``` matlab fig2_19.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.19
% Main lobe of beam pattern
% Kristine Bell
% Last updated 6/4/01, 10/4/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% bp_sect - plots conventional bp with sector
clear all
% close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Uniform Linear Array
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 10; % Elements in array
d = 0.5; % spacing wrt wavelength
beamwidth = 2/(N*d); % null-to-null (LL BW is half this)
D=d*(-(N-1)/2:1:(N-1)/2); % element locations
u = -0.45:0.01:0.45;
AS = exp(-1i*2*pi*D'*0); % BP points to 0
Au = exp(-1i*2*pi*D'*u);
B = real(AS'*Au)/N; % BP
h=plot(u,B,'-');
set(h,'LineWidth',1.5)
hold on
% HPBW
I=find(B>=0.707);
plot([u(min(I)-1) u(max(I)+1)],[B(min(I)-1) B(max(I)+1)],'-')
plot(u(min(I)-1)*[1 1],B(min(I)-1)*[1 1]+0.03*[-1 1],'-')
plot(u(max(I)+1)*[1 1],B(max(I)+1)*[1 1]+0.03*[-1 1],'-')
plot([0.03 0.13],[0.72 0.845],'-')
text(0.155,0.84,'$\Delta u_1 = $HPBW','Interpreter','latex');
% BW-NN
plot([-0.2 0.2],[-0.3 -0.3],'-')
plot(-0.2*[1 1],-0.3*[1 1]+0.03*[-1 1],'-')
plot(0.2*[1 1],-0.3*[1 1]+0.03*[-1 1],'-')
text(-0.05,-0.35,'$\Delta u_2 = BW_{NN}$','Interpreter','latex');
% axes
plot([-1 1 ],[0 0],'-')
plot([0 0],[-0 1.1],'-')
% tick marks
plot(-0.4*[1 1],0.03*[-1 1])
plot(-0.2*[1 1],0.03*[-1 1])
plot(0.2*[1 1],0.03*[-1 1])
plot(0.4*[1 1],0.03*[-1 1])
yy = -0.1;
%tick labels
text(-0.45,yy,'$-2\frac{\lambda}{ND}$','Interpreter','latex');
text(-0.21,yy,'$-\frac{\lambda}{ND}$','Interpreter','latex');
text(-0.005,yy,'0','Interpreter','latex');
text(0.18,yy,'$\frac{\lambda}{ND}$','Interpreter','latex');
text(0.38,yy,'$2\frac{\lambda}{ND}$','Interpreter','latex');
text(0.01,1.1,'$B(u)$','Interpreter','latex')
hold off
axis([-0.42 0.42 -0.4 1.2])
set(gca,'Visible','off')
```
{% endnote %}
&emsp;&emsp;3dB 波束宽度是波束宽度的一个度量。定义为对应 $\left| B_u(u) \right|^2 = 0.5$ 或 $\left| B_u(u) \right| = 1/\sqrt{2}$ 的点。可以通过令波束方向图 $B_u(u)$（式 $\eqref{BeamPatternOfUWLAInDirectionCosineDomain}$）等于 $1/\sqrt{2}$，来求出 $u$ 空间的半功率点。对于 $N\geqslant 10$，通过解下面的方程可以得到一个很好的近似：
$$\begin{equation}
@@ -295,6 +538,55 @@ $$\begin{equation}
![图 1-3-6 阵元间距对波束方向图的影响：(a) $d=\lambda/4$；(b) $d=\lambda/2$；(c) $d=\lambda$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-20-josh-oap-part-1-3/2023-04-20-josh-oap-part-1-3-060-%20EffectOfElementSpacingOnBeamPattern.svg){width=1000px}
{% note primary MATLAB Code %}
``` matlab fig2_20.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.20
% Effect of element spacing on beam pattern
% Xin Zhang
% Lillian Xu updated 09/2000
% updated by K. Bell 7/22/2001, 10/4/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N = 10;
n = (-(N-1)/2:(N-1)/2).';
psi = pi*(-3:0.001:3);
w = 1/N*[ones(N,1)];
d = [1/4 1/2 1];
figure
for k =1:length(d)
vv = exp(1i*2*d(k)*n*psi);
B(k,:) = 20*log10(abs(w'*vv));
subplot(3,1,k)
plot(psi/pi,B(k,:));
hold on
plot([-1 -1],[-25 0],'--');
plot([1 1],[-25 0],'--');
axis([-3 3 -25 5])
set(gca,'YTick',[-25 -20 -15 -10 -5 0])
grid on
xlabel('$u$','Interpreter','latex')
if k == 1
plot([-1 1],[2.5 2.5])
plot([-1 -0.9],[2.5 0.5])
plot([-1 -0.9],[2.5 4.5])
plot([0.9 1],[0.5 2.5])
plot([0.9 1],[4.5 2.5])
text(-0.4,7,'Visible region')
end
hold off
end
set(gcf,'Position',[0 0 800 600])
```
{% endnote %}
&emsp;&emsp;在[图 1-3-6](#fig.1-3-6) 中，针对多种 $d/\lambda$ 的值，画出了 $\left| \boldsymbol{\varUpsilon}_u(u) \right|$。图中说明了“栅瓣”这个重要的概念，即和主波束一样高的波瓣。栅瓣发生在当波束方向图 $B_\psi(\psi)$（式 $\eqref{BeamPatternOfUWLAInWaveNumberDomain}$）中的分子和分母均为零的时候。栅瓣出现的间隔为
$$\begin{equation}
source/_posts/2023-04-22-josh-oap-part-1-4.md
+161
View File
@@ -130,6 +130,48 @@ $$\begin{equation}
![图 1-4-2 阵列调向对栅瓣的影响：$N=10$，(a) $d=2\lambda/3, \bar\theta = 30^\circ$；(b) $d=\lambda/2, \bar\theta = 90^\circ$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-22-josh-oap-part-1-4/2023-04-22-josh-oap-part-1-4-020-%20EffectOfSteeringOnTheGratingLobes.svg){width=1000px}
{% note primary MATLAB Code%}
``` matlab fig2_22.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.22
% Effect of element spacing on beam pattern
% Lillian Xu
% Last updated by K. Bell 7/22/01, 10/4/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N=10;
n=(-(N-1)/2:(N-1)/2)';
w=ones(N,1)/N;
u=-3:0.01:3;
d=[2/3 1/2];
theta=[30 90];
figure
s = ['(a)';'(b)'];
for m=1:2
beam=(w.*exp(1i*2*pi*n*d(m)*sin(theta(m)/180*pi)))' ...
*exp(1i*2*pi*n*d(m)*u);
subplot(2,1,m);
plot(u,20*log10(abs(beam)));
hold on
axis([-3 3 -25 5])
ylabel('Beam pattern (dB)')
xlabel('$u$','Interpreter','latex')
text(0,-32,s(m,:),'HorizontalAlignment','center')
% title(['N=',num2str(N),', d=',num2str(d(m),2),' lambda, theta-Bar=', ...
% num2str(theta(m)),' degrees'])
grid on
end
set(gcf,'Position',[0 0 800 600])
```
{% endnote %}
&emsp;&emsp;通常，为了避免栅瓣进入可视区域，需要有
$$\begin{equation}
@@ -161,6 +203,38 @@ $$\begin{equation}
![图 1-4-3 10 阵元均匀阵列（$d = \lambda/2$）扫描到 30°（和正侧向夹角为 60°）时的波束方向图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-22-josh-oap-part-1-4/2023-04-22-josh-oap-part-1-4-030-BeamPatternFor10ElementUniformArrayScannedTo30Degree.svg){width=800px}
{% note primary MATLAB Code %}
``` matlab fig2_23.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.23
% Beam pattern for 10-element uniform array (d=lambda/2)
% scanned to 30 degrees (60 degrees from broadside)
% Xin Zhang 1/20/99
% Last updated by K. Bell 6/25/01
% Functions called: polardb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N = 10;
n = (-(N-1)/2:(N-1)/2).';
theta = pi*(-1:0.001:1);
u = cos(theta);
d = 1/2;
vv = exp(1i*2*pi*d*n*u);
theta_T = 30/180*pi;
w = 1/N*ones(N,1).*exp(1i*n*pi*cos(theta_T));
B = w'*vv;
B = 10*log10(abs(B).^2);
figure
h=polardb(theta,B,-40);
hold off
```
{% endnote %}
&emsp;&emsp;为了研究在 $\theta$ 空间的半功率带宽情况，我们利用式 $\eqref{SteeredBeamPatternInAngleDomain}$ 和 [$N \geqslant 10$ 时半功率带宽的估算式][]，可以得到在 $u$ 空间的左右半功率点分别为
$$\begin{equation}
@@ -196,6 +270,61 @@ $$\begin{equation} \label{HPBWOfUWLA}
![图 1-4-4 HPBW 和阵列调向角的关系：标准线阵，采用均匀加权](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-22-josh-oap-part-1-4/2023-04-22-josh-oap-part-1-4-040-HPBWVersusSteeringAngleStandardUWLA.svg){width=1000px}
{% note primary MATLAB Code %}
``` matlab fig2_24.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.24
% HPBW versus steering angle
% Lillian Xu
% Modified by Xin Zhang
% Last updated by L. Xu 11/30/00, K. Bell 7/22/01, 10/4/01, 10/17/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
close all;
Nd_set=[1000:-5:700 700:-0.5:200 200:-0.05:20 15:-0.005:1];
theta_deg=[2.5 5 10 20 30 45 90];
%theta_deg=[2.5 5 10 20 30 45 60 90];
theta_set=theta_deg*pi/180;
figure
for num=1:size(theta_set,2)
theta=theta_set(num);
bw=0;
for num1=1:size(Nd_set,2)
Nd=Nd_set(num1);
r=cos(theta)+0.443/Nd;
if abs(r)<=1
bw(num1)=acos(r)-acos(cos(theta)-0.443/Nd);
end
end
loglog(Nd_set(1:size(bw,2)),abs(bw*180/pi));
hold on
rightpoint=abs(bw(1)*180/pi);
text(1040,rightpoint,[num2str(theta_deg(num)),'$^\circ$'],Interpreter='latex'); % used to be 7
end
axis([1 1000 0.04 100])
xlabel('$Nd/\lambda$','Interpreter','latex')
ylabel('3-dB beamwidth in degrees')
scanlimit=-acos(1-2*0.443./Nd_set);
loglog(Nd_set,abs(scanlimit*180/pi),'--')
Endfire = 2*acos(1-0.433./Nd_set);
loglog(Nd_set,abs(Endfire*180/pi))
rightpoint=abs(Endfire(1)*180/pi);
text(1040,rightpoint,'$\theta=0^\circ$',Interpreter='latex'); % used to be 7
%title('HPBW versus steering angle: standard linear array with uniform weighting')
text(12,35,'Endfire')
text(12,14.4,'Scan limit')
text(12,1.8,'Broadside')
grid on;
set(gcf,'Position',[0 0 800 500])
```
{% endnote %}
&emsp;&emsp;当波束调向到正侧向的附近（$\bar{\theta}_\mathrm{T}$ 的值很小）且 $Nd \gg \lambda$ 时，$\theta_\mathrm{H}$ 很小，利用小角度展开形式得到
$$\begin{equation} \label{HPBWWithSmallAngle}
@@ -220,6 +349,38 @@ $$\begin{equation} \label{HPBWOfEndfireUWLA}
![图 1-4-6 一个标准 10 阵元线阵。采用均匀幅度加权时在端射时的波束方向图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-22-josh-oap-part-1-4/2023-04-22-josh-oap-part-1-4-060-BeamPatternOfStandard10ElementUWLAAtEndfire.svg){width=800px}
{% note primary MATLAB Code %}
``` matlab fig2_26.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 2.26
% Beam pattern of a standard 10-element linear array
% with uniform amplitude weighting at end-fire
% Xin Zhang 1/20/99
% Last updated by K. Bell 6/25/01
% Functions called: polardb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N = 10;
n = (-(N-1)/2:(N-1)/2).';
theta = pi*(-1:0.001:1);
u = cos(theta);
d = 1/2;
vv = exp(1i*2*pi*d*n*u);
theta_T = 0/180*pi;
w = 1/N*ones(N,1).*exp(1i*n*pi*cos(theta_T));
B = w'*vv;
B = 10*log10(abs(B).^2);
figure
h=polardb(theta,B,-40);
hold off
```
{% endnote %}
可以把式 $\eqref{HPBWOfEndfireUWLA}$ 写成
$$\begin{equation}
source/_posts/2023-06-19-josh-oap-part-2-1.md
+271 -4
View File
@@ -21,7 +21,7 @@ abbrlink: 61f0754c
date: 2023-06-19 16:25:26
---
&emsp;&emsp;下面主要研究对每个传感器输出进行加权的技术，以获得具有理想性质的频率-波数响应和波束方向图。加权的过程在一些文献中也称为锐化 (tapering) 或赋形 (shading) 。许多阵列加权技术等价于时间序列谱分析中的加窗或锐化 (Taper) 技术。
&emsp;&emsp;谱加权（Spectral Weighting）技术主要研究对如何每个传感器输出进行加权，以获得具有理想性质的频率-波数响应和波束方向图。加权的过程在一些文献中也称为锐化 (tapering) 或赋形 (shading) 。许多阵列加权技术等价于时间序列谱分析中的加窗或锐化 (Taper) 技术。
&emsp;&emsp;对于一个阵元间距小于或等于 $\lambda/2$ 的线阵，由于不需要考虑栅瓣和空域混迭，可以直接应用连续函数的傅里叶变换理论对权重和频率-波数响应之间存在的傅里叶变换关系进行分析。本文关注线性阵列，讨论如何在降低旁瓣的同时使主瓣宽度的增长最小（最小化方向性损失），将尝试不同的权值矢量，并分析它们的性能。
@@ -111,6 +111,50 @@ $$\begin{equation}
![图 2-1-1 均匀加权和 Cosine 加权的波束方向图: $N = 11$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-06-19-josh-oap-part-2-1/2023-06-19-josh-oap-part-2-1-010-BeamPatternForUniformAndCosineWeighting.svg){width=700px}
{% note primary MATLAB Code %}
``` matlab fig3_01.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 3.1
% Cosine weighting: Array beam patterns
% Kristine Bell
% Modified by Xin Zhang 3/17/99
% Lillian 04/16/2001, K. Bell 7/20/01, K. Bell 9/30/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
M=11;
d=0.5; % sensor spacing wrt wavelength
D = (-(M-1)/2:1:(M-1)/2)*d; % sensor positions in wavelengths
% weights, normalized so that w(0)=1
W_unf = ones(1,M);
W_cos = cos(pi*D*2/M);
% Beampatterns
u = 0:0.001:1;
A = exp(-1i*2*pi*D.'*u);
G_unf = W_unf*A;
G_unf = G_unf/(max(abs(G_unf)));
G_cos = W_cos*A;
G_cos = G_cos/(max(abs(G_cos)));
figure
% array only
clf
plot(u,20*log10(abs(G_unf)),'--')
hold on
plot(u,20*log10(abs(G_cos)),'-')
hold off
axis([0 1 -80 0])
%title([int2str(M) ' element array'])
h=legend('Uniform','Cosine');
xlabel('$u$','Interpreter','latex')
ylabel('Beam pattern (dB)')
```
{% endnote %}
# 3. 升 Cosine 加权
&emsp;&emsp;可以把矩形均匀加权和 Cosine 加权组合起来，获得两种加权方法一些理想性质的组合。对应的阵列加权为
@@ -129,7 +173,7 @@ $$\begin{equation}
&emsp;&emsp;在[图 2.1.2](#fig.2-1-2) 中给出了 $p = 0.31, 0.17, 0$ 时的升 Cosine 加权的波束方向图。当 $p$ 减小时，第一旁瓣的高度减小，主波束宽度增加。波束方向图参数如下：
| $p$ | HPBW | $BW_{NN}$ | 第一旁瓣高度 | $D_N$[^1] |
| $p$ | HPBW | $BW_{NN}$ | 第一旁瓣高度 | $D_N$ |
| :--: | :--: | :--: | :--: | :--: |
| 0.31 | $1.03 \displaystyle\frac{2}{N}$ | $2.50 \displaystyle\frac{2}{N}$ | $-20.0 \text{dB}$ | $0.928$ |
| 0.17 | $1.09 \displaystyle\frac{2}{N}$ | $3.70 \displaystyle\frac{2}{N}$ | $-22.0 \text{dB}$ | $0.886$ |
@@ -139,6 +183,52 @@ $$\begin{equation}
![图 2-1-2 升 Cosine 加权的波束方向图: $N = 11$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-06-19-josh-oap-part-2-1/2023-06-19-josh-oap-part-2-1-020-BeamPatternForRaisedCosineWeighting.svg){width=700px}
{% note primary MATLAB Code %}
``` matlab fig3_02.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 3.2
% Raised cosine weighting function
% Kristine Bell
% Modified by Xin Zhang 3/17/99
% Lillian Xu 04/16/2001, K. Bell 7/20/01, 9/30/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
M=11;
d=0.5; % sensor spacing wrt wavelength
D = (-(M-1)/2:1:(M-1)/2)*d; % sensor positions in wavelengths
% weights, normalized so that w(0)=1
p=[0;0.17;0.31;1];
W_rcos = p*ones(1,M)+(ones(4,1)-p)*cos(pi*D*2/M);
% Beampatterns
u = 0:0.001:1;
c=p+(ones(4,1)-p)*2/pi;
A = exp(-1i*2*pi*D.'*u);
G_rcos = W_rcos*A;
for m=1:4
G_rcos(m,:) = G_rcos(m,:)/(max(abs(G_rcos(m,:))));
end
figure
clf
plot(u,20*log10(abs(G_rcos(3,:))),'-.')
hold on
plot(u,20*log10(abs(G_rcos(2,:))),'--')
plot(u,20*log10(abs(G_rcos(1,:))),'-')
hold off
axis([0 1 -80 0])
h=legend({'$p=0.31$','$p=0.17$','$p=0$'},'Interpreter','latex');
xlabel('$u$','Interpreter','latex')
ylabel('Beam pattern (dB)')
```
{% endnote %}
&emsp;&emsp;至此，我们已经可以使 HPBW 变窄，并使得第一旁瓣比均匀分布的情况要低得多。
# 4. Cosine$^m$ 加权
@@ -157,7 +247,7 @@ $$\begin{equation}
&emsp;&emsp;$m=2$ 时的加权有时称为是 Hann 加权[^2]。同样，波束方向图可以利用 Cosine 函数的指数形式来计算。波束方向图在[图 2.1.3](#fig.2-1-3) 中给出。当 $m$ 增加时，旁瓣减小，但是主波束变宽。波束方向图的参数为:
| $m$ | HPBW | $BW_{NN}$ | 第一旁瓣高度 | $D_N$[^1] |
| $m$ | HPBW | $BW_{NN}$ | 第一旁瓣高度 | $D_N$ |
| :--: | :--: | :--: | :--: | :--: |
| 2 | $1.44 \displaystyle\frac{2}{N}$ | $4 \displaystyle\frac{2}{N}$ | $-31.4 \text{dB}$ | $0.667$ |
| 3 | $1.66 \displaystyle\frac{2}{N}$ | $5 \displaystyle\frac{2}{N}$ | $-39.4 \text{dB}$ | $0.576$ |
@@ -167,6 +257,56 @@ $$\begin{equation}
![图 2-1-3 Cosine$^m$ 加权的波束方向图: $N = 11$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-06-19-josh-oap-part-2-1/2023-06-19-josh-oap-part-2-1-030-BeamPatternForCosineMWeighting.svg){width=700px}
{% note primary MATLAB Code %}
``` matlab fig3_03.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 3.3
% Cosine^m weighting: m = 2,3,and 4
% Kristine Bell
% Modified by Xin Zhang 3/17/99
% Lillian Xu 04/16/2001, K. Bell 7/20/01, 9/30/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
M=11;
d=0.5; % sensor spacing wrt wavelength
D = (-(M-1)/2:1:(M-1)/2)*d; % sensor positions in wavelengths
% weights, normalized so that w(0)=1
W_cos2 = cos(pi*D*2/M).^2;
W_cos3 = cos(pi*D*2/M).^3;
W_cos4 = cos(pi*D*2/M).^4;
% Beampatterns
u = 0:0.001:1;
A = exp(-1i*2*pi*D.'*u);
G_cos2 = W_cos2*A;
G_cos2 = G_cos2/(max(abs(G_cos2)));
G_cos3 = W_cos3*A;
G_cos3 = G_cos3/(max(abs(G_cos3)));
G_cos4 = W_cos4*A;
G_cos4 = G_cos4/(max(abs(G_cos4)));
figure
clf
plot(u,20*log10(abs(G_cos2)),'--')
hold on
plot(u,20*log10(abs(G_cos3)),'-')
plot(u,20*log10(abs(G_cos4)),'-.')
hold off
axis([0 1 -80 0])
h=legend({'$m=2$','$m=3$','$m=4$'},'Interpreter','latex');
set(h,'Fontsize',12)
xlabel('$u$','Interpreter','latex')
ylabel('Beam pattern (dB)')
```
{% endnote %}
# 5. 升 Cosine 平方加权
&emsp;&emsp;升 Cosine 平方系列的加权为
@@ -200,12 +340,139 @@ $$\begin{equation}
B_u(u) = 0.54 \frac{\sin\left(\displaystyle\frac{N\pi u}{2}\right)}{\sin\left(\displaystyle\frac{\pi u}{2}\right)} + 0.23 \left[ \frac{\sin\left( \displaystyle\frac{N\pi}{2} \left( u - \frac{2}{N} \right) \right)}{\sin\left( \displaystyle\frac{\pi}{2} \left( u - \frac{2}{N} \right) \right)} + \frac{\sin\left( \displaystyle\frac{N\pi}{2} \left( u + \frac{2}{N} \right) \right)}{\sin\left( \displaystyle\frac{\pi}{2} \left( u + \frac{2}{N} \right) \right)} \right]
\end{equation}$$
&emsp;&emsp;[图 2.1.4](#fig.2-1-4) 给出了 Hamming 加权、Hann 加权和 Blackman-Harris 加权的波束方向图。第一个零点的位置在 $u = 4/N$，第一个非零旁瓣的高度为 $-39.5\text{dB}$。加权有步进不连续性，使得旁瓣的衰减渐近为 $O(1/n)$。但是，第一个旁瓣被完全对消掉了，剩余的旁瓣都很低，尽管它们的衰减比较慢。而且，注意到 Hamming 窗的波束宽度小于 Hann 加权的波束宽度，且第一旁瓣更低。这是一般变化趋势的一个例外情况。
&emsp;&emsp;[图 2.1.4](#fig.2-1-4) 给出了 Hamming 加权、Hann 加权和 Blackman-Harris 加权的波束方向图。第一个零点的位置在 $u = 4/N$，第一个非零旁瓣的高度为 $-39.5\text{dB}$。由于加权具有步进不连续性（step discontinuity），旁瓣的衰减渐近为 $O(1/n)$。但是，第一个旁瓣被完全对消掉了，尽管剩余的旁瓣衰减得比较慢，它们的幅度都不是很高。而且注意到 Hamming 窗的波束宽度小于 Hann 加权的波束宽度，且第一旁瓣更低。
<a id="fig.2-1-4"></a>
![图 2-1-4 Hann 加权、Hamming 加权和 Blackman-Harris 加权的波束方向图: $N = 11$](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-06-19-josh-oap-part-2-1/2023-06-19-josh-oap-part-2-1-040-BeamPatternForHannHammingAndBlackmanHarrisWeighting.svg){width=700px}
{% note primary MATLAB Code %}
``` matlab fig3_04.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figure 3.4
% Beam patterns for Hann, Hamming, and Blackman-Harris
% Kristine Bell
% Modified by Xin Zhang 3/17/99
% Lillian Xu 04/16/2001, K. Bell 7/20/01, 9/30/01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hann, Hamming, Blackman-Harris
% M-element array, d=lambda/2
clear all
close all
M=11;
d=0.5; % sensor spacing wrt wavelength
D = (-(M-1)/2:1:(M-1)/2)*d; % sensor positions in wavelengths
% weights, normalized so that w(0)=1
W_Hann = 0.5*(ones(1,M)+cos(2*pi*D*2/M));
W_Hamm = 0.54*ones(1,M)+0.46*cos(2*pi*D*2/M);
W_BH = 0.42*ones(1,M)+0.5*cos(2*pi*D*2/M)+0.08*cos(4*pi*D*2/M);
% Beampatterns
u = 0:0.001:1;
A = exp(-1i*2*pi*D.'*u);
G_Hann = W_Hann*A;
G_Hann = G_Hann/(max(abs(G_Hann)));
G_Hamm = W_Hamm*A;
G_Hamm = G_Hamm/(max(abs(G_Hamm)));
G_BH = W_BH*A;
G_BH = G_BH/(max(abs(G_BH)));
figure
clf
plot(u,20*log10(abs(G_Hann)),'--')
hold on
plot(u,20*log10(abs(G_Hamm)),'-')
plot(u,20*log10(abs(G_BH)),'-.')
hold off
axis([0 1 -80 0])
h=legend('Hann','Hamming','Blackman-Harris');
xlabel('$u$','Interpreter','latex')
ylabel('Beam pattern (dB)')
```
{% endnote %}
# 7. Blackman-Harris 加权
&emsp;&emsp;Blackman-Harris 加权把加权过程扩展到了更高阶谐波的情况，并在前两个旁瓣的峰值位置上放置了零点。加权函数为
$$\begin{equation}
w(\tilde{n}) = 0.42 + 0.5 \cos\left( \frac{2\pi \tilde{n}}{N} \right) + 0.08 \cos\left( \frac{4\pi \tilde{n}}{N} \right), \quad \tilde{n} = -\frac{N-1}{2}, \cdots, \frac{N-1}{2}
\end{equation}$$
对应的波束方向图是多个常规波束方向图的和，已在[图 2.1.4](#fig.2-1-4) 中给出
$$\begin{equation}
\begin{aligned}
B_u(u) =& 0.42 \frac{\sin\left( \displaystyle\frac{N\pi u}{2} \right)}{\sin\left( \displaystyle\frac{\pi u}{2} \right)} + 0.25 \left[ \frac{\sin\left( \displaystyle\frac{N\pi}{2} \left( u - \frac{2}{N} \right) \right)}{\sin\left( \displaystyle\frac{\pi}{2} \left( u - \frac{2}{N} \right) \right)} + \frac{\sin\left( \displaystyle\frac{N\pi}{2} \left( u + \frac{2}{N} \right) \right)}{\sin\left( \displaystyle\frac{\pi}{2} \left( u + \frac{2}{N} \right) \right)} \right] \\
&+ 0.04 \left[ \frac{\sin\left( \displaystyle\frac{N\pi}{2} \left( u - \frac{2}{N} \right) \right)}{\sin\left( \displaystyle\frac{\pi}{2} \left( u - \frac{2}{N} \right) \right)} + \frac{\sin\left( \displaystyle\frac{N\pi}{4} \left( u + \frac{4}{N} \right) \right)}{\sin\left( \displaystyle\frac{\pi}{4} \left( u + \frac{4}{N} \right) \right)} \right]
\end{aligned}
\end{equation}$$
波束方向图的参数为
| 加权方式 | HPBW | $BW_{NN}$ | 第一旁瓣高度 | $D_N$ |
| :--: | :--: | :--: | :--: | :--: |
| Hann | $1.44 \displaystyle\frac{2}{N}$ | $4.0 \displaystyle\frac{2}{N}$ | $-31.4 \text{dB}$ | $0.667$ |
| Hamming | $1.31 \displaystyle\frac{2}{N}$ | $4.0 \displaystyle\frac{2}{N}$ | $-39.5 \text{dB}$ | $0.730$ |
| Blackman-Harris | $1.65 \displaystyle\frac{2}{N}$ | $6.0 \displaystyle\frac{2}{N}$ | $-56.6 \text{dB}$ | $0.577$ |
至此，所有的加权方法都是基于各种正弦函数进行的。现在看看其他类型的阵列加权技术，试试能否更好在波束宽度和旁瓣之间的获得折中。
# 8. 长球函数 (Prolate Spheroidal Functions)
&emsp;&emsp;在设计加权方案时，我们的目标是使得一个给定角度范围内的能量占总能量的比例最大，也即最大化如下比值
$$\begin{equation}
\alpha = \frac{\displaystyle\iint_{\varOmega_1}\left| B(\theta,\varphi) \right|^2 \sin\theta\,\mathrm{d}\theta\,\mathrm{d}\varphi}{\displaystyle\int_0^\pi \int_0^{2\pi}\left| B(\theta,\varphi) \right|^2 \sin\theta\,\mathrm{d}\varphi\,\mathrm{d}\theta}
\end{equation}$$
其中 $\varOmega_1$ 是主波束周围的一个区域。对于标准线性阵列，可以写成
$$\begin{equation}
\alpha = \frac{\displaystyle\int_{-\psi_0}^{\psi_0}\left| B(\psi) \right|^2\,\mathrm{d}\psi}{\displaystyle\int_{-\pi}^{\pi}\left| B(\psi) \right|^2\,\mathrm{d}\psi}
\end{equation}$$
其分子和分母部分可分别写为
$$\begin{equation}
\begin{aligned}
&\alpha_N &&\xlongequal{B(\psi) = \boldsymbol{w}^\mathrm{H}\boldsymbol{v}_\psi(\psi)} \int_{-\psi_0}^{\psi_0}\boldsymbol{w}^\mathrm{H}\boldsymbol{v}_\psi(\psi) \boldsymbol{v}_\psi^\mathrm{H}(\psi)\boldsymbol{w}\,\mathrm{d}\psi \\
& &&= \boldsymbol{w}^\mathrm{H} \left[ \int_{-\psi_0}^{\psi_0} \boldsymbol{v}_\psi(\psi) \boldsymbol{v}_\psi^\mathrm{H}(\psi) \,\mathrm{d}\psi \right]\boldsymbol{w} \\
& &&= \boldsymbol{w}^\mathrm{H}\boldsymbol{A}\boldsymbol{w} \\
&\alpha_D &&\xlongequal{B(\psi) = \boldsymbol{w}^\mathrm{H}\boldsymbol{v}_\psi(\psi)} \int_{-\pi}^{\pi}\boldsymbol{w}^\mathrm{H}\boldsymbol{v}_\psi(\psi) \boldsymbol{v} \\
& &&= \boldsymbol{w}^\mathrm{H}\boldsymbol{A}\boldsymbol{w}
\end{aligned}
\end{equation}$$
$$\begin{equation}
\end{equation}$$
其中
$$\begin{equation}
\begin{aligned}
\boldsymbol{A} &\triangleq \int_{-\psi_0}^{\psi_0} \boldsymbol{v}_\psi(\psi) \boldsymbol{v}_\psi^\mathrm{H}(\psi) \,\mathrm{d}\psi \\
\boldsymbol{B} &\triangleq \int_{-\pi}^{\pi} \boldsymbol{v}_\psi(\psi) \boldsymbol{v}_\psi^\mathrm{H}(\psi) \,\mathrm{d}\psi = 2\pi \boldsymbol{I}
\end{aligned}
\end{equation}$$
$$\begin{equation}
[\boldsymbol{A}]_{mn} = \int_{-\psi_0}^{\psi_0} e^{jm\psi} e^{-jn\psi} \,\mathrm{d}\psi = \frac{2\sin ((m-n)\psi_0)}{m-n} = 2\psi_0 \,\mathrm{sinc}((m-n)\psi_0)
\end{equation}$$
因此
$$\begin{equation}
\alpha = \frac{\boldsymbol{w}^\mathrm{H}\boldsymbol{A}\boldsymbol{w}}{2\pi \boldsymbol{w}^\mathrm{H}\boldsymbol{w}}
\end{equation}$$
[^1]: 参数 $D_N$ 是阵列的归一化方向性。它是相对于一个均匀加权阵列的方向性进行归一化的。对于标准线阵, $D_N$ 也是归一化的白噪声增益。