Phân tích không gian pha bao gồm việc xác định thời gian độ trễ \( t \) và số chiều nhúng \( m \). Phương pháp thông tin tương hỗ trợ để xác định thời gian độ trễ. Chúng ta sẽ tính toán số chiều liên kết của dữ liệu một chiều để xác minh đặc trưng hỗn loạn.

Khi phân tích tín hiệu điện não bộ, tôi đã phát hiện ra một hiện tượng thú vị: những biến động dường như ngẫu nhiên có thể chứa các luật định tính. Kỹ thuật phân tích không gian pha giúp nâng dữ liệu chuỗi thời gian một chiều lên không gian nhiều chiều, tiết lộ các đặc trưng động lực học của hệ thống gốc. Bây giờ chúng ta sẽ thực hành phân tích không gian pha bằng MATLAB, từ việc lựa chọn thời gian độ trễ, xác định số chiều nhúng cho đến xác minh đặc trưng hỗn loạn.

Bắt đầu bằng việc chuẩn bị dữ liệu thử nghiệm, ở đây tôi sử dụng hệ thống Lorenz tiêu chuẩn để tạo chuỗi hỗn loạn:

% Tạo dữ liệu hệ thống Lorenz
sigma = 10; beta = 8/3; rho = 28;
f = @(t,x) [sigma*(x(2)-x(1)); x(1)*(rho-x(3))-x(2); x(1)*x(2)-beta*x(3)];
[t,x] = ode45(f,[0:0.01:100],[1;1;1]);
dulieu = x(:,1); % Lấy thành phần x làm dữ liệu thử nghiệm

Bước 1: Sử dụng thông tin tương hỗ trợ tìm thời gian độ trễ \( \tau \)

Hàm tương quan tự do chỉ phản ánh mối quan hệ tuyến tính, trong khi thông tin tương hỗ trợ có thể nắm bắt mối quan hệ phi tuyến tính. Khi thông tin tương hỗ trợ đạt cực tiểu đầu tiên, điều này cho thấy dữ liệu tại độ trễ đó mang lại sự độc lập lớn nhất.

function tgian = thongtin_tuonghoptu(dulieu, gioihan_tau)
    do_dai = length(dulieu);
    thongtin = zeros(1,gioihan_tau);
    
    for t = 1:gioihan_tau
        dichchuyen = dulieu(t+1:end);
        goc = dulieu(1:end-t);
        
        % Đếm phân phối xác suất liên hợp bằng biểu đồ cột hai chiều
        [P_xy,canh] = histcounts2(goc, dichchuyen, 'BinMethod','sqrt');
        P_xy = P_xy / sum(P_xy(:));
        
        % Tính xác suất biên
        P_x = sum(P_xy,2);
        P_y = sum(P_xy,1);
        
        % Tính thông tin tương hỗ trợ
        hoply = P_xy > 0;
        thongtin(t) = sum(P_xy(hoply) .* log2(P_xy(hoply)./(P_x(hoply(:,1)) .* P_y(hoply(:,2))')));
    end
    
    % Tìm cực tiểu đầu tiên
    [~,tgian] = findpeaks(-thongtin, 'NPeaks',1);
end

Gọi hàm `tgian = thongtin_tuonghoptu(dulieu, 50)` để nhận được thời gian độ trễ tốt nhất. Phần chính của mã là ước lượng phân phối xác suất liên hợp thông qua biểu đồ cột, khi độ lệch giữa chuỗi dịch chuyển và chuỗi gốc giảm tối đa, \( \tau \) lúc đó là độ trễ tối ưu.

Bước 2: Sử dụng phương pháp hàng xóm giả mạo để xác định số chiều nhúng \( m \)

Khi tăng số chiều không còn giảm đáng kể số điểm hàng xóm giả mạo thì không gian pha đã được mở rộng đầy đủ.

function sochiendoi = hangxom_giamao(dulieu, tgian, gioihan_m)
    do_dai = length(dulieu);
    ty_le_giamao = zeros(1,gioihan_m);
    
    for chieudai=1:gioihan_m
        % Xây dựng không gian pha
        du_lieu_phang = nhung(dulieu, chieudai, tgian);
        
        % Tìm hàng xóm gần nhất của mỗi điểm
        [~,khoangcach1] = knnsearch(du_lieu_phang(1:end-1,:), du_lieu_phang(2:end,:));
        
        % Thay đổi khoảng cách khi tăng thêm một chiều
        du_lieu_phang_tieptheo = nhung(dulieu, chieudai+1, tgian);
        [~,khoangcach2] = knnsearch(du_lieu_phang_tieptheo(1:end-1,:), du_lieu_phang_tieptheo(2:end,:));
        
        % Tính tỷ lệ hàng xóm giả mạo
        giamao = abs(khoangcach2 - khoangcach1) ./ khoangcach1 > 0.15;
        ty_le_giamao(chieudai) = sum(giamao)/length(giamao);
    end
    
    % Xác định m khi tỷ lệ ngừng giảm đáng kể
    sochiendoi = find(diff(ty_le_giamao) < 0.05, 1);
end

function du_lieu_phang = nhung(dulieu, chieudai, tgian)
    do_dai = length(dulieu);
    du_lieu_phang = zeros(do_dai-(chieudai-1)*tgian, chieudai);
    for i=1:chieudai
        du_lieu_phang(:,i) = dulieu((1:do_dai-(chieudai-1)*tgian) + (i-1)*tgian);
    end
end

Tại đây, khi thay đổi khoảng cách vượt quá 15% thì được coi là hàng xóm giả mạo. Trong thực tế, ngưỡng này có thể điều chỉnh dựa trên đặc tính của dữ liệu, thường nằm trong khoảng 10%-20%.

Bước 3: Xác minh hỗn loạn bằng số chiều liên kết

Hiện tượng bão hòa của số chiều liên kết là dấu hiệu của hệ thống hỗn loạn, đối lập với sự tăng vĩnh cửu của quá trình ngẫu nhiên.

function D2 = chieudai_lienket(dulieu, tgian, chieudai)
    du_lieu_phang = nhung(dulieu, chieudai, tgian);
    do_dai = size(du_lieu_phang,1);
    rs = logspace(log10(0.1*std(dulieu)), log10(0.5*std(dulieu)), 20);
    C = zeros(size(rs));
    
    for k=1:length(rs)
        r = rs(k);
        % Tính tích phân liên kết
        ma_tran_khoangcach = pdist2(du_lieu_phang, du_lieu_phang);
        C(k) = sum(ma_tran_khoangcach(:) < r) / (do_dai*(do_dai-1));
    end
    
    % Phù hợp với khu vực tuyến tính
    he_so = diff(log(C))./diff(log(rs));
    D2 = mean(he_so(5:15)); % Lấy khu vực ổn định giữa
end

Kết quả cho thấy số chiều liên kết tụ hội ở mức 2.05, rõ ràng thấp hơn số chiều nhúng 3, cho thấy hệ thống có một hấp dẫn tử hỗn loạn chiều thấp. Nếu là nhiễu ngẫu nhiên, số chiều sẽ tiếp tục tăng theo số chiều nhúng.

Ví dụ Áp Dụng

Tính toán số chiều liên kết của tín hiệu EEG đo được:

eeg = load('eeg_data.mat').signal; 
tgian = thongtin_tuonghoptu(eeg, 30);
chieudai = hangxom_giamao(eeg, tgian, 8);
D2 = chieudai_lienket(eeg, tgian, chieudai);

Khi \( D2 \) ổn định ở mức 3-5, cho thấy có đặc trưng phát ban nổ của loạn thần kinh. Chỉ số phi tuyến tính này tiết lộ các đặc tính phi ổn định của tín hiệu não bộ tốt hơn so với phân tích phổ truyền thống.

Một số lời khuyên thực hành:

1. Độ dài dữ liệu phải ít nhất là 100*(\( \tau \)*(\( m \)-1)), nếu không kết quả phân tích kém.
2. Khi tính thông tin tương hỗ trợ có thể sử dụng ước lượng mật độ hạt nhân thay vì biểu đồ cột.
3. Trong việc tính số chiều liên kết có thể thêm cửa sổ Theiler để loại bỏ các điểm liên quan theo thời gian.

Phân tích không gian pha giống như trang kính đa chiều giúp chúng ta nhìn thấy các biến động hỗn loạn dưới góc độ mới, biến những biến động hỗn loạn trở thành cấu trúc hấp dẫn tử tinh vi. Sự chuyển đổi từ bề mặt đến bản chất là sức mạnh của phân tích phi tuyến.