matlab代码实现MCMC波达角估计

MCMC（Markov Chain Monte Carlo）波达角估计是一种基于概率的估计方法，可以用于估计信号在空间中的传播方向，以下是一个MATLAB实现的示例代码：

% 设置参数
M = 10000;      % MCMC样本数
N = 100;        % 数据点数
T = 1;          % 采样时间
L = T*N;        % 信号长度
theta_true = 45;% 真实波达角
sigma = 0.1;    % 噪声标准差

% 生成数据
s = sin(2*pi*1/T*(0:L-1));       % 信号源
s = s-mean(s);
s = s/norm(s);
a = exp(1j*(0:180)*pi/180);      % 阵列响应
x = real(a(theta_true+1)'*s)'+sigma*randn(L,1); % 观测信号

% MCMC波达角估计
theta = zeros(M,1);  % 存储MCMC样本
theta(1) = randi(180);% 初始值随机取
for i = 2:M
    theta_prop = mod(theta(i-1)+randn*10,180); % 提议分布
    a_prop = exp(1j*(0:180)*pi/180*theta_prop);% 提议分布对应的阵列响应
    x_prop = real(a_prop'*s)'+sigma*randn(L,1);% 提议分布对应的观测信号
    % 计算接受率
    alpha = min(1,exp(-norm(x-x_prop)^2/2/sigma^2));
    % 生成接受/拒绝样本
    if rand<alpha
        theta(i) = theta_prop;
    else
        theta(i) = theta(i-1);
    end
end

% 结果展示
figure;
subplot(2,1,1);
plot(x);
xlabel('采样点');
ylabel('幅度');
title('观测信号');
subplot(2,1,2);
histogram(theta,0:180);
xlabel('波达角');
ylabel('频数');
title('波达角估计直方图');

该代码中，首先生成一个正弦信号源，并通过阵列响应模拟信号在不同波达角下的接收情况，再加上高斯噪声生成观测信号。然后，使用MCMC算法对波达角进行估计，其中提议分布为均匀分布，接受率计算使用高斯噪声的概率密度函数。最后，使用直方图展示波达角的估计结果。