MATLAB实现复数信号的盲源分离

function [Ae, Se] = jade_complex(X, n_sources)
% 复数JADE盲源分离算法
% 输入: X-混合信号矩阵, n_sources-源信号数量
% 输出: Ae-估计的混合矩阵, Se-估计的源信号

    [m, T] = size(X);
    if nargin < 2
        n_sources = m;
    end
    
    % 1. 数据中心化
    X_centered = X - mean(X, 2);
    
    % 2. 数据白化
    [X_white, W_white] = whiten_data_complex(X_centered, n_sources);
    
    % 3. 计算四阶累积量矩阵
    cumulant_matrices = compute_cumulant_matrices(X_white, n_sources);
    
    % 4. 联合近似对角化
    V = joint_diagonalization(cumulant_matrices, n_sources);
    
    % 5. 估计分离矩阵和源信号
    W = V' * W_white;
    Ae = pinv(W);
    Se = W * X_centered;
end

function [X_white, W_white] = whiten_data_complex(X, n_sources)
% 复数数据白化
    [m, T] = size(X);
    
    % 计算协方差矩阵
    C = (X * X') / T;
    
    % 特征值分解
    [V, D] = eig(C);
    [d, idx] = sort(real(diag(D)), 'descend');
    V = V(:, idx);
    
    % 保留主要成分
    V = V(:, 1:n_sources);
    d = d(1:n_sources);
    
    % 白化矩阵
    W_white = diag(1./sqrt(d)) * V';
    X_white = W_white * X;
end

function cumulant_matrices = compute_cumulant_matrices(X_white, n_sources)
% 计算四阶累积量矩阵
    [m, T] = size(X_white);
    cumulant_matrices = zeros(m, m, n_sources);
    
    for i = 1:n_sources
        % 选择参考信号
        y = X_white(i, :);
        
        % 计算四阶矩矩阵
        Q = (X_white .* conj(y)) * X_white' / T;
        M = Q - eye(m) - (y * y') / T * (X_white * X_white') / T;
        
        cumulant_matrices(:, :, i) = M;
    end
end

function V = joint_diagonalization(cumulant_matrices, n_sources)
% 联合近似对角化
    V = eye(n_sources);
    max_iter = 100;
    tolerance = 1e-10;
    
    for iter = 1:max_iter
        max_off = 0;
        
        for p = 1:n_sources-1
            for q = p+1:n_sources
                % 构建2x2子问题
                G = zeros(2, 2);
                
                for k = 1:size(cumulant_matrices, 3)
                    M = cumulant_matrices(:, :, k);
                    g11 = M(p, p) - M(q, q);
                    g12 = M(p, q) + M(q, p);
                    g22 = M(q, q) - M(p, p);
                    
                    G(1, 1) = G(1, 1) + abs(g11)^2;
                    G(1, 2) = G(1, 2) + g11' * g12;
                    G(2, 1) = G(2, 1) + g12' * g11;
                    G(2, 2) = G(2, 2) + abs(g12)^2 + abs(g22)^2;
                end
                
                % 特征值分解求最优旋转
                [U, D] = eig(G);
                [~, idx] = min(real(diag(D)));
                angle = U(:, idx);
                
                % Givens旋转矩阵
                Givens = eye(n_sources);
                Givens(p, p) = angle(1);
                Givens(p, q) = angle(2);
                Givens(q, p) = -conj(angle(2));
                Givens(q, q) = conj(angle(1));
                
                % 更新累积量矩阵和分离矩阵
                for k = 1:size(cumulant_matrices, 3)
                    cumulant_matrices(:, :, k) = Givens' * cumulant_matrices(:, :, k) * Givens;
                end
                
                V = V * Givens;
                
                % 计算非对角元素
                off_diag = 0;
                for k = 1:size(cumulant_matrices, 3)
                    M = cumulant_matrices(:, :, k);
                    off_diag = off_diag + sum(sum(abs(M - diag(diag(M))).^2));
                end
                max_off = max(max_off, off_diag);
            end
        end
        
        if max_off < tolerance
            break;
        end
    end
end

function [W, S] = fastica_complex(X, n_sources, max_iter)
% 复数FastICA算法实现
% 输入: X-混合信号, n_sources-源数, max_iter-最大迭代次数

    if nargin < 3
        max_iter = 200;
    end
    if nargin < 2
        n_sources = size(X, 1);
    end
    
    [m, T] = size(X);
    
    % 中心化
    X_centered = X - mean(X, 2);
    
    % 白化
    [X_white, W_white] = whiten_data_complex(X_centered, n_sources);
    
    % 初始化分离矩阵
    W = randn(n_sources, n_sources) + 1j * randn(n_sources, n_sources);
    W = orth(W')';
    
    % 复数非线性函数
    g = @(x) tanh(real(x)) + 1j * tanh(imag(x));
    g_prime = @(x) (1 - tanh(real(x)).^2) + 1j * (1 - tanh(imag(x)).^2);
    
    % ICA迭代
    for iter = 1:max_iter
        W_old = W;
        
        for i = 1:n_sources
            w = W(i, :);
            
            % 期望计算
            wx = w * X_white;
            g_wx = g(wx);
            g_prime_wx = g_prime(wx);
            
            term1 = mean(X_white .* conj(g_wx), 2);
            term2 = mean(g_prime_wx) * w';
            
            w_new = term1 - term2;
            
            % 归一化
            w_new = w_new / norm(w_new);
            
            % 去相关
            if i > 1
                w_new = w_new - W(1:i-1, :)' * (W(1:i-1, :) * w_new);
                w_new = w_new / norm(w_new);
            end
            
            W(i, :) = w_new';
        end
        
        % 检查收敛
        if norm(abs(W * W_old') - eye(n_sources), 'fro') < 1e-6
            break;
        end
    end
    
    % 分离信号
    W = W * W_white;
    S = W * X_centered;
end

function [S, A] = fajd_tvar_bss(X, n_sources, order, n_basis)
% 基于时变自回归模型的FAJD盲源分离算法
% 输入: X-混合信号, n_sources-源数, order-模型阶数, n_basis-基函数数量

    [m, T] = size(X);
    
    % 白化预处理
    [X_white, W_white] = whiten_data_complex(X, n_sources);
    
    % TVAR模型系数估计
    coefficient_matrices = estimate_tvar_coefficients(X_white, order, n_basis);
    
    % 快速近似联合对角化
    V = fajd(coefficient_matrices);
    
    % 信号分离
    W = V' * W_white;
    A = pinv(W);
    S = W * X;
end

function coefficient_matrices = estimate_tvar_coefficients(X, order, n_basis)
% 估计时变自回归模型系数
    [m, T] = size(X);
    t = linspace(0, 1, T);
    
    % 基函数（Legendre多项式）
    basis_functions = zeros(n_basis, T);
    for i = 1:n_basis
        basis_functions(i, :) = legendre_poly(i-1, t);
    end
    
    coefficient_matrices = zeros(m, m, n_basis * order);
    
    % 递推最小二乘法估计系数
    for i = 1:m
        for j = 1:m
            % 构建回归矩阵
            Phi = zeros(T - order, order * n_basis);
            for p = 1:order
                for k = 1:n_basis
                    col_idx = (p-1)*n_basis + k;
                    Phi(:, col_idx) = basis_functions(k, order+1:T) .* ...
                                     X(j, 1:T-order)';
                end
            end
            
            % 最小二乘估计
            y = X(i, order+1:T)';
            theta = pinv(Phi) * y;
            
            % 存储系数矩阵
            for p = 1:order
                for k = 1:n_basis
                    idx = (p-1)*n_basis + k;
                    coefficient_matrices(i, j, idx) = theta(idx);
                end
            end
        end
    end
end

% 复数盲源分离完整示例
clear; close all; clc;

%% 1. 生成测试信号
T = 2000; % 样本点数
t = linspace(0, 1, T);

% 源信号1: 复数QPSK信号
f1 = 10;
s1 = exp(1j * 2*pi*f1*t + 1j*pi/4*(round(4*t)));

% 源信号2: 复数调频信号
f2 = 5 + 3*sin(2*pi*t);
s2 = exp(1j * 2*pi*f2.*t);

% 源信号3: 复数噪声信号
s3 = (randn(1, T) + 1j*randn(1, T)) / sqrt(2);

S = [s1; s2; s3];

%% 2. 混合信号
m = 3; % 传感器数量
A = (randn(m, 3) + 1j*randn(m, 3)) / sqrt(2); % 随机混合矩阵
X = A * S + 0.01*(randn(m, T) + 1j*randn(m, T)); % 添加噪声

%% 3. 使用JADE算法进行盲源分离
[A_jade, S_jade] = jade_complex(X, 3);

%% 4. 使用FastICA算法进行盲源分离
[W_fastica, S_fastica] = fastica_complex(X, 3);

%% 5. 结果评估和显示
figure('Position', [100, 100, 1200, 800]);

% 原始信号
subplot(4, 3, 1); plot(real(S(1, :))); title('源信号1 (实部)');
subplot(4, 3, 2); plot(real(S(2, :))); title('源信号2 (实部)');
subplot(4, 3, 3); plot(real(S(3, :))); title('源信号3 (实部)');

% 混合信号
subplot(4, 3, 4); plot(real(X(1, :))); title('混合信号1 (实部)');
subplot(4, 3, 5); plot(real(X(2, :))); title('混合信号2 (实部)');
subplot(4, 3, 6); plot(real(X(3, :))); title('混合信号3 (实部)');

% JADE分离结果
subplot(4, 3, 7); plot(real(S_jade(1, :))); title('JADE分离1 (实部)');
subplot(4, 3, 8); plot(real(S_jade(2, :))); title('JADE分离2 (实部)');
subplot(4, 3, 9); plot(real(S_jade(3, :))); title('JADE分离3 (实部)');

% FastICA分离结果
subplot(4, 3, 10); plot(real(S_fastica(1, :))); title('FastICA分离1 (实部)');
subplot(4, 3, 11); plot(real(S_fastica(2, :))); title('FastICA分离2 (实部)');
subplot(4, 3, 12); plot(real(S_fastica(3, :))); title('FastICA分离3 (实部)');

%% 6. 性能评估
fprintf('性能评估:\n');

% JADE性能
G_jade = A_jade * A;
G_jade = G_jade ./ max(abs(G_jade), [], 2);
fprintf('JADE - 全局矩阵的非对角元素能量: %.6f\n', ...
    norm(G_jade - diag(diag(G_jade)), 'fro'));

% FastICA性能  
G_fastica = pinv(W_fastica) * A;
G_fastica = G_fastica ./ max(abs(G_fastica), [], 2);
fprintf('FastICA - 全局矩阵的非对角元素能量: %.6f\n', ...
    norm(G_fastica - diag(diag(G_fastica)), 'fro'));

%% 7. 星座图显示（适用于通信信号）
figure('Position', [100, 100, 1000, 600]);

subplot(2, 3, 1);
plot(real(S(1, 1:100:end)), imag(S(1, 1:100:end)), '.');
title('源信号1星座图'); axis equal;

subplot(2, 3, 2);
plot(real(X(1, 1:100:end)), imag(X(1, 1:100:end)), '.');
title('混合信号1星座图'); axis equal;

subplot(2, 3, 3);
plot(real(S_jade(1, 1:100:end)), imag(S_jade(1, 1:100:end)), '.');
title('JADE分离信号1星座图'); axis equal;

subplot(2, 3, 4);
plot(real(S(2, 1:100:end)), imag(S(2, 1:100:end)), '.');
title('源信号2星座图'); axis equal;

subplot(2, 3, 5);
plot(real(X(2, 1:100:end)), imag(X(2, 1:100:end)), '.');
title('混合信号2星座图'); axis equal;

subplot(2, 3, 6);
plot(real(S_jade(2, 1:100:end)), imag(S_jade(2, 1:100:end)), '.');
title('JADE分离信号2星座图'); axis equal;

% 重要参数调优示例
function optimal_params = parameter_tuning(X)
% 自动参数调优函数
    [m, T] = size(X);
    
    % 源数估计
    n_sources = estimate_source_number(X);
    
    % 信噪比估计
    snr_est = estimate_snr(X);
    
    % 根据信号特性选择算法
    if is_communication_signal(X)
        optimal_params.algorithm = 'jade';
        optimal_params.n_sources = n_sources;
    elseif is_nonstationary_signal(X)
        optimal_params.algorithm = 'fajd_tvar';
        optimal_params.order = 4; % TVAR模型阶数
        optimal_params.n_basis = 3; % 基函数数量
    else
        optimal_params.algorithm = 'fastica';
        optimal_params.n_sources = n_sources;
        optimal_params.max_iter = 500;
    end
end