空域的匹配滤波

对于一个固定的阵型的阵列来说，理想远场情况下，声源处发出声音到MIC收到声音信号之间存在延迟，如果阵元位置位于相同的位置，则阵元与阵元之间的延时差就等于零，此时如果将声源按照固定的距离为半径绕阵元位置一周，并扫描各个方向上信号的累加值，可以发现阵列输出是各向同性的。

如果将阵元在一条直线上均匀散开，间距为d，同样将声源按照固定的距离为半径，绕阵列中心点扫描一周，因为延迟差的存在，各阵元在特定方向的累加值就与信号的输入角度有关，可以猜测信号与阵列平面垂直（定义此时入射角为0度，角度按逆时针递增）时候，阵列输出最大。

如果考虑声音在空气中传播时候的振幅衰减，扫描初始波束形状如下图：

从波束扫描的图形来看，可以看出当声源来自于垂直于阵列平面方向的时候，阵列的输出是最高的。

如果我们在阵元的各路信号中，对阵元的相位响应进行一定角度的调整，就可以将阵列的响应波束按照我们期望的方向改变。如我们将阵元的输出按照入射角45度的延时向量进行调整，可以得到以下波束：

如果在形成波束的时候，考虑声音传播过程中的幅值衰减，得到波束图形如下：

验证代码：

clc
clear 
close all
format long

%% 阵列的信号入射角度与输出之间的关系
% 阵列平面上五个阵元，间距0.2m
% 一个信号源S1 = 1*exp（a1 + jwt）,距离阵列中心点距离3.4m
% 阵列的输出波束形状与输入角度（信号延迟）有关，与输入信号和阵列之间的距离（声压衰减速度）有关
% 所以我们对声源信号模型进行近场、远场分类，分情况讨论

speedOfAcoustic = 340;
distanceOfElement = 0.2;
frequency = 1000;
w = 2*pi*frequency;
a1 = 0;
length = 3.4;

% 取两个周期的信号
t = 1/(32*frequency):1/(32*frequency):64/(32*frequency);
S1 = exp(1j*w*t);
figure(1)
plot(real(S1))

A = zeros(360,64);

% 计算信号到每个阵元的延迟
% 信源与阵列中心点连线，与阵列平面法线之间的角度关系为 dgree1(逆时针,使用角度，求解时转换为弧度)
% 信源到阵列平面的距离h/length = cos(degree1弧度)
% 信源到法线的距离d/length = sin(degree1弧度)

degreeSwitch = 45;
    h = length * cos(degreeSwitch*pi/180);
    d = length * sin(degreeSwitch*pi/180);
    
    length1 = ((2*distanceOfElement - d)^2 + h^2)^0.5;
    length2 = ((distanceOfElement - d)^2 + h^2)^0.5;
    length3 = ((d)^2 + h^2)^0.5;
    length4 = ((-1*distanceOfElement - d)^2 + h^2)^0.5;
    length5 = ((-2*distanceOfElement - d)^2 + h^2)^0.5;
    
    % 模拟权值向量w对相位以及波束指向的调整
    t11 = length1/speedOfAcoustic;
    t22 = length2/speedOfAcoustic;
    t33 = length3/speedOfAcoustic;
    t44 = length4/speedOfAcoustic;
    t55 = length5/speedOfAcoustic;
    
    lengthTotal = zeros(5,360);

for degree1 = 1:360
    h = length * cos(degree1*pi/180);
    d = length * sin(degree1*pi/180);
    
    length1 = ((2*distanceOfElement - d)^2 + h^2)^0.5;
    lengthTotal(1,degree1) = length1;
    length2 = ((distanceOfElement - d)^2 + h^2)^0.5;
    lengthTotal(2,degree1) = length2;
    length3 = ((d)^2 + h^2)^0.5;
    lengthTotal(3,degree1) = length3;
    length4 = ((-1*distanceOfElement - d)^2 + h^2)^0.5;
    lengthTotal(4,degree1) = length4;
    length5 = ((-2*distanceOfElement - d)^2 + h^2)^0.5;
    lengthTotal(5,degree1) = length5;
    
    t1 = length1/speedOfAcoustic;
    t2 = length2/speedOfAcoustic;
    t3 = length3/speedOfAcoustic;
    t4 = length4/speedOfAcoustic;
    t5 = length5/speedOfAcoustic;
    
%     A(degree1,:) = exp(1j*(w*(t - t1))) + exp(1j*(w*(t - t2))) + exp(1j*(w*(t - t3))) ...
%         + exp(1j*(w*(t - t4))) + exp(1j*(w*(t - t5)));                    
%     A(degree1,:) = (length3/length1)^2*exp(1j*(w*(t - t1))) + (length3/length2)^2*exp(1j*(w*(t - t2))) + (length3/length3)^2*exp(1j*(w*(t - t3))) ...
%         + (length3/length4)^2*exp(1j*(w*(t - t4))) + (length3/length5)^2*exp(1j*(w*(t - t5)));
%     A(degree1,:) = exp(1j*(w*(t - t1 + t11))) + exp(1j*(w*(t - t2 + t22))) + exp(1j*(w*(t - t3 + t33))) ...
%           + exp(1j*(w*(t - t4 + t44))) + exp(1j*(w*(t - t5 + t55)));
    A(degree1,:) = (length3/length1)^2*exp(1j*(w*(t - t1 + t11))) + (length3/length2)^2*exp(1j*(w*(t - t2 + t22))) + (length3/length3)^2*exp(1j*(w*(t - t3 + t33))) ...
          + (length3/length4)^2*exp(1j*(w*(t - t4 + t44))) + (length3/length5)^2*exp(1j*(w*(t - t5 + t55)));
end

absA = abs(A(:,1));
gain = 20*log10(absA/max(absA));
figure(2)
polarplot((1:360)*pi/180,gain')
rlim([min(gain),0])
title('调整45度波束形状')
figure(3)
plot((1:360)*pi/180,gain')
ylim([min(gain),0])
figure(4)
hold on
plot(lengthTotal(1,:))
plot(lengthTotal(2,:))
plot(lengthTotal(3,:))
plot(lengthTotal(4,:))
plot(lengthTotal(5,:))
hold off
legend('1','2','3','4','5')