最基本的麦克风阵列延时累加和算法，DelayAndSum beamformer

最近从 phased.TimeDelayBeamformer 扒了它的代码，简单改造作为最基本的麦克风阵列算法学习

refer to https://www.mathworks.com/help/phased/ref/phased.timedelaybeamformer-system-object.html

一、延时累加和beamformer

可以说是最基本的麦克风阵列算法，

 

 方法很简单，就是通过平移每个通道数据到合适的位置，然后累计调整系数即可。

但问题是什么？

一般采集到的麦克风数据都是数字化的，采样点的间隔是1/fs 如果要平移的时间量刚好是N个sample点，那么刚好是有对应的 采样值的 如果平移量不是，那么得考虑下面两种方案：

1、时域上做插值，根据周围的采样点来插值插出来你想要的某些平移量的数据

2、转换到频域上操作

频域上的 延时，只要对信号 相位做调整即可，

 

具体操作：通过fft转到频域，执行相位变换，再ifft回到时域

 二、具体的matlab代码

clear all
close all
clc;

% phased.TimeDelayBeamformer

fs = 8000;

%入射角度，或者叫阵列想指向哪里，旋转角和相位角，其中相位角是离开xy平面到z轴的角度
incidentAngle = [-50;30];

%阵列采集数据
load rx

%麦克风单元数

M=11;

%阵列拓扑结构，其他任意阵列形状也可

ArrayPos = [0    0    0    0    0    0    0    0    0    0    0;
-0.2    -0.16    -0.12    -0.08    -0.04    0    0.04    0.08    0.12    0.16    0.2;
0    0    0    0    0    0    0    0    0    0    0];

%旋转角和俯仰角，

azang = incidentAngle(1,:);
elang = incidentAngle(2,:);

%导向矢量
incidentdir = [-cosd(elang).*cosd(azang);...
    -cosd(elang).*sind(azang);...
    -sind(elang)];

%延时时间

delay = ArrayPos.'*incidentdir/c;
pSampleRate = fs;

% x_steered = delayseq(rx,-delay,pSampleRate);

% 初始化
x=rx;
inputLength = size(x,1);

%刚好是“反“移动，对阵列本身的校准，这样就达到延迟一致了
dt = -delay;
dt = dt *fs;

%取整数延时部分（sample点数）
delayInt = round(dt);    % Integer delays in samples

%分数延迟部分（sample点数，单位是个）
delayFrac = dt - delayInt;

%找到最大延迟点，数据长度因为延迟而膨胀了
maxLength = inputLength+max(0, max(delayInt)); % maximum sequence length
% Define upperbound
maxLengthLimit = maxLength;

%找到合适的fft size

nfftLimit = 2^nextpow2(2*inputLength);
output = zeros(maxLengthLimit,length(dt));

% Perform delay operation，对每个延时通道
for colI=1:length(dt)

    endIdx = inputLength+delayInt(colI);
    % no operation if delayed or advanced out of scope
    if (endIdx <= 0) || (endIdx > 2*inputLength)
        continue;
    end
    tmpx = x(:,colI);
    

%如果有分数个延时点的要求，转换到频域
    if delayFrac(colI)    
        nfft = 2^nextpow2(inputLength + max(0, delayInt(colI)));
        assert(nfft <= nfftLimit);
        binStart = floor(nfft/2);
        % Notice the FFT bins must belong to [-pi, pi].
        fftBin = 2*pi*ifftshift(((0:nfft-1)-binStart).')/nfft;
        if isscalar(tmpx) && nfft==1
            tmpxd = tmpx;
        else%核心代码
            tmpxd = fft(tmpx,nfft);
            tmpxd = ifft(tmpxd(:).*exp(-1i*dt(colI)*fftBin));
        end
        %以下认为正向延迟的头几个点不要，置零

        %你可以认为从负时间轴（t<0) 的地方带(移）不过来数据，都是零

        if delayInt(colI) >= 0
            orgStart = delayInt(colI) + 1;
            newStart = delayInt(colI) + 1;
        else
            orgStart = 1;
            newStart = 1;
        end
        orgEnd = endIdx;
        if isreal(x)
           output(newStart:endIdx,colI) = real(tmpxd(orgStart:orgEnd));
        else
           output(newStart:endIdx,colI) = tmpxd(orgStart:orgEnd);
        end
    else
        % Integer sampling shifting
        if delayInt(colI) >= 0
            orgStart = 1;
            newStart = delayInt(colI) + 1;
        else
            orgStart = 1-delayInt(colI);
            newStart = 1;
        end
        orgEnd = inputLength;
        output(newStart:endIdx,colI) = tmpx(orgStart:orgEnd);
    end
    
end

% output has the same length as input
x_steered = output(1:inputLength, :);
w  = ones(M,1)/M;% 归一化权重
y = x_steered*conj(w);%相乘作和

plot(t,rx(:,6),'r:',t,y)
xlabel('Time (sec)')
ylabel('Amplitude')
legend('Original','Beamformed')