matlab室内麦克风阵列仿真

一、整体信号模型

M 个麦克风，K 个声源（含噪源），第 m 通道输出：

\[x_m(t) = \sum_{k=1}^{K} s_k(t) * h_{mk}(t) + n_{\text{sensor},m}(t) \]

• (\(h_{mk}\))：声源 k → 麦克风 m 的 RIR（含混响+直达）
• 噪声可以是独立传感噪（各通道不相关），也可以是空间相关噪

---

二、房间 + 阵列 + 声源几何定义

%% room_mic_setup.m
function [room, array, src] = room_mic_setup()
% 房间 [Lx Ly Lz] (m)
room.L = [6, 5, 3];

% 反射系数 -> RT60（简化 Sabine: RT60 ≈ 0.161V / (Σα_i S_i)）
room.alpha = 0.3;       % 墙面平均吸声系数 (0.1~0.6，住宅~0.3)
room.c = 343;           % 声速
room.fs = 16000;

% === 阵列：4 元线性阵，沿 x 轴 ===
array.M = 4;
array.d = 0.04;         % 间距 4 cm（< λ/2 @ 4kHz）
array.pos = zeros(4,3);
for m = 1:4
    array.pos(m,:) = [(m-1)*array.d - 1.5*array.d + 2, 2, 1.2];
    % 阵列中心 (2,2,1.2)m，离地 1.2m
end

% === 目标声源（说话人）===
src(1).pos = [3, 3.5, 1.6];   % 目标说话人
src(1).type = 'point';

% === 噪声源 1：风扇，指定角度（远场近似，方位 30°）===
src(2).pos = [];               % 远场用 DOA 代替
src(2).az = deg2rad(30);       % 方位角（x 轴正向为 0°）
src(2).el = deg2rad(0);        % 俯仰
src(2).type = 'farfield';      % 平面波，只加延时+幅度

% === 噪声源 2：角落点源（扩散场近似用多点）===
src(3).pos = [5.5, 0.5, 1.0];
src(3).type = 'point';
end

想要环形阵 / 方形阵 / 随机扰动阵，把 array.pos 换掉即可，后面 RIR 和 DOA 延时都只认坐标。

---

三、RIR 生成

没 Audio Toolbox 也能跑，有 Toolbox 的话直接换 roomImpulseResponse 更快。

%% ism_rir.m — 图像法 RIR（矩形房，单音源→单麦克风）
function h = ism_rir(src_pos, mic_pos, room, max_order)
% src_pos, mic_pos: [x y z]
% room.L = [Lx Ly Lz], room.c, room.alpha, room.fs
% max_order: 镜像阶数（3~5 够家用，>7 很慢）

if nargin < 4, max_order = 4; end

L = room.L; c = room.c; fs = room.fs;
alpha = room.alpha;   % 简化：所有墙同 α

% 镜像坐标枚举（1D 笛卡尔积）
lims = -max_order:max_order;
[idx_x, idx_y, idx_z] = meshgrid(lims, lims, lims);
idx_x = idx_x(:); idx_y = idx_y(:); idx_z = idx_z(:);

h = zeros(ceil(0.5*fs)+1, 1);  % RIR 最长 0.5s

for k = 1:length(idx_x)
    ox = idx_x(k); oy = idx_y(k); oz = idx_z(k);
    
    % 镜像源位置
    mx = src_pos(1) + 2*ox*L(1);
    if ox < 0, mx = -mx; end   % 修正符号（标准 ISM）
    my = src_pos(2) + 2*oy*L(2);
    if oy < 0, my = -my; end
    mz = src_pos(3) + 2*oz*L(3);
    if oz < 0, mz = -mz; end
    
    % 距离
    r = norm([mx, my, mz] - mic_pos);
    tau = r / c;
    n = round(tau * fs) + 1;
    if n > length(h), continue; end
    
    % 反射衰减：(alpha)^(|ox|+|oy|+|oz|)
    refl = alpha^(abs(ox)+abs(oy)+abs(oz));
    
    % 球面扩散 + 镜像衰减
    h(n) = h(n) + refl / (4*pi*r);
end

% 归一化 + 加微小传感器噪
h = h / max(abs(h)) * 0.8;
h = h + 1e-6*randn(size(h));   % 数值稳定
end

严格 ISM 要对每面墙用反射系数 R=(Z2-Z1)/(Z2+Z1)，上面用统一 alpha 是住宅级的合理近似。要更准可换成各墙独立 alpha_x+, alpha_x-, ...。

---

四、远场噪声 DOA → 各通道延时

如果噪声是远场平面波（风扇、空调、窗外车流），不给位置只给 (az, el)，按阵列几何算延时：

%% farfield_delay.m
function [delays, att] = farfield_delay(array, az, el)
% az: 方位角，x 轴正向=0°, 逆时针 +
% el: 俯仰角
% 返回: delays (samples), att (幅度，球面衰减 ≈1 远场)

M = array.M;
u = [cos(el)*cos(az), cos(el)*sin(az), sin(el)]';  % 传播方向单位矢
c = 343; fs = array.fs;

delays = zeros(M,1);
att   = ones(M,1);
for m = 1:M
    % 波程差 = -r_m·u（负号：波从 u 方向来）
    d_tau = -array.pos(m,:) * u / c;
    delays(m) = d_tau * fs;   % 可正可负，小数样点
end
% 以首麦为参考偏移到非负
delays = delays - min(delays);
end

小数延时用 fracDelay 或 interp1 实现：

function y = apply_frac_delay(x, delay_samples)
% x: 输入, delay_samples: 标量或向量(每通道)
% 用 FIR 分数延时（简单版用 FFT 相位偏移更稳）
Nf = 64;  b = fir1(Nf, 0.4);  % 低通
% 这里给 FFT 版（更快）
N = length(x);
X = fft(x);
f = (0:N-1)/N*2*pi;
H = exp(-1j*f*delay_samples);
y = real(ifft(X .* H));
end

---

五、主仿真脚本（混响 + 任意角度噪声）

%% main_indoor_mic_sim.m
clear; clc; close all;

fs = 16000;
[room, array, src] = room_mic_setup();
room.fs = fs;

%% ========== 1. 读干净语音 ==========
% 用内置示例或你自己 wav
[clean, fs_orig] = audioread('speech.wav');
if fs_orig ~= fs, clean = resample(clean, fs, fs_orig); end
clean = clean(:,1);                    % 单通道
clean = clean / max(abs(clean)) * 0.7; % 归一
N = length(clean);

%% ========== 2. 目标源 → 各麦 RIR 卷积 ==========
M = array.M;
x_clean = zeros(N, M);
for m = 1:M
    h = ism_rir(src(1).pos, array.pos(m,:), room, 4);
    x_clean(:,m) = conv(clean, h, 'same');
end

%% ========== 3. 加噪声（两种模式演示）==========
% --- 噪声 a：远场 30° 风扇（平面波，走 DOA 延时）---
noise_wave = 0.15 * randn(N,1);   % 白噪近似风扇底噪
[delays_a, ~] = farfield_delay(array, src(2).az, src(2).el);
x_noise_a = zeros(N, M);
for m = 1:M
    x_noise_a(:,m) = apply_frac_delay(noise_wave, delays_a(m));
end

% --- 噪声 b：角落点源（走 RIR，空间相关强）---
x_noise_b = zeros(N, M);
for m = 1:M
    hb = ism_rir(src(3).pos, array.pos(m,:), room, 3);
    nb = 0.2 * randn(N,1);       % 点源激励
    x_noise_b(:,m) = conv(nb, hb, 'same');
end

% --- 传感器自噪 ---
sensor_noise = 0.01 * randn(N, M);

%% ========== 4. 按 SNR 混合 ==========
% 以第 1 麦干净信号能量为基准
E_s = mean(x_clean(:,1).^2);
target_SNR_db = 10;   % 总噪（a+b+传感器）
E_n_total = E_s / (10^(target_SNR_db/10));
% 归一化各噪贡献
E_a = mean(x_noise_a(:).^2);
E_b = mean(x_noise_b(:).^2);
E_sens = mean(sensor_noise(:).^2);
scale = sqrt(E_n_total / (E_a+E_b+E_sens));

x_total = x_clean + scale*(x_noise_a + x_noise_b + sensor_noise);

%% ========== 5. 听 + 画 ==========
soundsc(x_total(:,1), fs)   % 听第 1 麦

figure('Color','white','Position',[100 100 1000 500])
subplot(2,2,1)
plot(x_clean(:,1))
title('Clean (Mic1, reverberant)'); xlabel('Sample'); grid on

subplot(2,2,2)
plot(x_total(:,1))
title(sprintf('Noisy (SNR≈%d dB)', target_SNR_db)); xlabel('Sample'); grid on

subplot(2,2,3)
% 画阵列+源位置（俯视）
plot(array.pos(:,1), array.pos(:,2), 'bo', 'MarkerSize',10,'LineWidth',2)
hold on
plot(src(1).pos(1), src(1).pos(2), 'r^', 'MarkerSize',12,'LineWidth',2)
plot(src(3).pos(1), src(3).pos(2), 'ks', 'MarkerSize',10,'LineWidth',2)
% 远场噪用箭头标 DOA
quiver(3,4, cos(src(2).az), sin(src(2).az), 'r--','LineWidth',1.5)
legend('Array','Target Src','Point Noise','Far-field Noise DOA')
axis equal; grid on; title('Geometry (Top View)'); xlabel('x'); ylabel('y')

subplot(2,2,4)
% 第 1 麦频谱
Nfft = 1024;
[PX,f] = pwelch(x_total(:,1), hann(Nfft), Nfft/2, Nfft, fs);
plot(f, 10*log10(PX))
xlabel('Freq (Hz)'); ylabel('PSD (dB)'); title('Mic1 Spectrum'); grid on
xlim([0 4000])

参考代码 针对室内麦克风阵列仿真信号的产生，有混响，可以添加任意角度的噪声 www.youwenfan.com/contentcnw/82663.html

六、扩散场噪声

如果要模拟多人嘈杂/会议室扩散场，比"多点点源+大 α"更稳的做法是特征分解法生成空间相关噪：

%% diffuse_noise.m — 生成扩散场（各通道相关系数 sinc(2πf d / c)）
function nd = diffuse_noise(N, fs, array)
% nd: N×M
M = array.M;
nd = randn(N, M);
% FFT 域逐频施加空间相关
Nd = fft(nd, N);
f = (0:N-1)'/N*fs;
for k = 2:N/2+1
    % 理论扩散场协方差：R_ij = sinc(2πf(k) * norm(pi-pj)/c)
    R = zeros(M);
    for i = 1:M
        for j = 1:M
            d_ij = norm(array.pos(i,:)-array.pos(j,:));
            R(i,j) = sinc(2*pi*f(k)*d_ij/343);
        end
    end
    % 特征分解保证半正定
    [V,D] = eig((R+R')/2);
    L = real(diag(D));
    L(L<0) = 0;
    W = V * diag(sqrt(L)) * V';
    Nd(k,:) = Nd(k,:) * W;            % 单频点
    Nd(N-k+2,:) = conj(Nd(k,:));      % 共轭对称
end
nd = real(ifft(Nd, N));
end

调用时把 x_noise_a 换成 diffuse_noise(N,fs,array) 即可，声源定位论文里"diffuse babble"基本就这个做法。

---

七、和现有工具对照

工具	特点
本程序	透明、可改 DOA/几何/RIR 阶数、轻量
MATLAB roomImpulseResponse (Audio TB)	官方 ISM，支持球房间，快
PyRoomAcoustics (Python)	Python 圈标配，有 pyroomacoustics.simulate
image_source_method (FileExchange)	老牌，但只单通道 RIR
gpuRIR	GPU 加速，大规模阵列