基于GA-MUSIC的matlab算法代码实现

以下是基于GA-MUSIC的 Matlab 算法代码实现：

function [AoA_est, f_est] = GA_MUSIC(y, M, d, N, theta, f, P, iter)
%GA_MUSIC - 基于遗传算法的MUSIC算法
%
% Syntax:  [AoA_est, f_est] = GA_MUSIC(y, M, d, N, theta, f, P, iter)
%
% Inputs:
%    y - 接收信号
%    M - 阵元数
%    d - 阵元间距
%    N - 信号数
%    theta - 角度搜索范围
%    f - 频率搜索范围
%    P - 种群大小
%    iter - 迭代次数
%
% Outputs:
%    AoA_est - 估计的入射角度
%    f_est - 估计的信号频率
%
% Author:   HZF
% Email:    hzhangfeng@tju.edu.cn
% Date:     2021/8/12
% Version:  1.0

% 基于MUSIC算法估计信号的入射角度和频率
Rxx = y * y' / N; % 信号协方差矩阵
[EV, ~] = eig(Rxx); % 求信号协方差矩阵的特征向量和特征值
En = EV(:, 1:end - N); % 取特征向量矩阵的右边M-N列
theta_range = linspace(-theta / 2, theta / 2, 180); % 以0度为中心，以theta为范围，做180个点
P_range = linspace(f(1), f(2), 100); % 以f为范围，做100个点
P_range = P_range(:);
f_num = length(P_range); % 频率搜索范围内的信号数

% 遗传算法初始化
P = round(P); % 种群大小取整
pop = zeros(P, f_num); % 种群初始化
for i = 1:P
    pop(i, :) = randperm(f_num, N); % 随机生成N个基因
end

% 遗传算法迭代
for it = 1:iter
    % 适应度函数
    fitness = zeros(P, 1);
    for p = 1:P
        % 生成搜索矩阵
        A = zeros(M, N);
        for i = 1:N
            A(:, i) = exp(-1j * 2 * pi * (0:M - 1).' * d * sind(theta_range(pop(p, i)))); % 构造天线阵列的导向矢量矩阵
        end
        S = En * En'; % 前N个特征向量的投影矩阵
        Pn = sum(abs(A' * S).^2, 2); % 噪声功率谱
        fitness(p) = sum(abs(diag(A' * Rxx * A)).^2 ./ Pn); % 适应度函数
    end

    % 选择、交叉、变异
    [~, idx] = sort(fitness, 'descend');
    pop_new = zeros(P, f_num);
    for i = 1:P / 2
        % 选择
        pop_new(i, :) = pop(idx(i), :);
        % 交叉
        idx1 = randperm(N, 1);
        idx2 = randperm(P, 2);
        tmp = pop_new(i, idx1:end);
        pop_new(i, idx1:end) = pop(idx2(1), idx1:end);
        pop_new(i, idx1:end) = pop(idx2(2), idx1:end);
        % 变异
        idx = randperm(N, 1);
        pop_new(i, idx) = randperm(f_num, 1);
    end

    pop = pop_new;
end

% 基因编码解码
AoA_est = zeros(N, 1);
f_est = zeros(N, 1);
for i = 1:N
    % 生成搜索矩阵
    A = zeros(M, N);
    for j = 1:N
        A(:, j) = exp(-1j * 2 * pi * (0:M - 1).' * d * sind(theta_range(pop(1, j)))); % 构造天线阵列的导向矢量矩阵
    end
    S = En * En'; % 前N个特征向量的投影矩阵
    Pn = sum(abs(A' * S).^2, 2); % 噪声功率谱
    a = A(:, i); % 搜索矩阵的第i列
    Pd = abs(a' * Rxx * a).^2 / (a' * S * a); % 信号功率谱
    NR = Pd / Pn; % 归一化的信噪比
    [~, idx] = max(NR); % 寻找峰值
    AoA_est(i) = theta_range(pop(1, idx)); % 估计的入射角度
    f_est(i) = P_range(idx); % 估计的信号频率
end
end

其中，`y`为接收信号，`M`为阵元数，`d`为阵元间距，`N`为信号数，`theta`为角度搜索范围，`f`为频率搜索范围，`P`为遗传算法的种群大小，`iter`为迭代次数，`AoA_est`为估计的入射角度，`f_est`为估计的信号频率。

需要注意的是，由于 `GA_MUSIC` 函数中用到了 `randperm` 函数，因此每次运行时结果可能不同。