矩阵束估计二维频率的MATLAB代码

以下是一份MATLAB代码，用于实现矩阵束估计二维频率：

function [f_hat, A_hat] = matrix_pencil(X, K)

% X为观测数据矩阵，K为估计频率个数
[N, M] = size(X);

% 构造基矩阵
J1 = repmat((0:N-1)', 1, M);
J2 = repmat(0:M-1, N, 1);
J = J1 .* J2;

% 计算信号子空间
[U, ~, ~] = svd(X);
Un = U(:, K+1:end);

% 计算矩阵P和矩阵Q
P = Un * Un';
Q = X * P;

% 奇异值分解
[Uq, Sq, Vq] = svd(Q);

% 估计频率和振幅
f_hat = sort(angle(diag(Uq(:, 1:K)' * J * Vq(:, 1:K))), 'ascend');
A_hat = abs(Uq(:, 1:K)' * X);

end

该函数输入观测数据矩阵X和估计的频率个数K，并输出估计的频率f_hat和振幅A_hat。注意，该代码实现的是二维频率估计，因此观测数据矩阵X应为二维矩阵。