1.软件版本
matlab2021a
2.本算法理论知识
FAST方法，步骤如下所示：
3.核心代码
clc;
clear;
close all;
warning off;
addpath 'func\'
rng(3);
%读取数据选择
%标准信号
load edf_data.mat
%含伪迹信号
load edf_data1.mat
%%
%fast ICA算法
y01 = sensor3C;
y02 = sensor4C;
%信号归一化
data= data/max(abs(data));
y01 = sensor3C/max(abs(sensor3C));
y02 = sensor4C/max(abs(sensor4C));
figure;
subplot(211);
plot(y01);
title('观测信号1');
subplot(212);
plot(y02);
title('观测信号2');
y0 = [y01;y02];
%%
%步骤1：去掉均值
data_avg= mean(data);
y01_avg = mean(y01);
y02_avg = mean(y02);
data = data - data_avg;
y01 = y01-y01_avg;
y02 = y02-y02_avg;
data = func_smooth2(data,1);
y01 = func_smooth2(y01,1);
y02 = func_smooth2(y02,1);
y1 = [y01;y02];
%%
%步骤2：白化
%协方差
y1_cov = cov(y1');
%特征值分解
[E,D] = eig(y1_cov);
%白化矩阵Q
Q = inv(sqrt(D))*(E)';
%白化后信号
y1_white = Q*y1;
y2 = y1_white;
figure;
subplot(211);
plot(y2(1,:));
title('白化后信号1');
subplot(212);
plot(y2(2,:));
title('白化后信号2');
%FASTICA算法
data_fica = y2;
[Rnum,Cnum] = size(data_fica);
Register_Matrix = zeros(Rnum,Rnum);
Error2 = cell(1,Rnum);
W_sum2 = cell(1,Rnum);
Err2 = cell(1,Rnum);
for ker=1:Rnum%观测信号个数
index =1;
%设置最大迭代次数
%
global_error = 1e-10;
maxIteration = 500;
%初始值
Wt = randn(Rnum,1)/4;
Wt = Wt/norm(Wt);
Err = [];
Error = [];
W_sum = [];
while index <= maxIteration | abs(abs(Wt'*Wt)-1) < global_error
W_n_1 = Wt;
Ka = 1;
u = 0.4;
t = data_fica'*Wt;
yexp = exp(-Ka*t.^2/2);
g = t.*yexp;
s1 =(1-Ka*t.^2);
s2 = exp(-Ka*t.^2/2);
dg = s1.*s2;
%加权求和判决
if sum(abs(t'*g*Wt)) > sum(abs(data_fica*g))
Wt_sum=(1-u)*t'*g*Wt + u*data_fica*g;
else
Wt_sum=(u)*t'*g*Wt + (1-u)*data_fica*g;
end
tmps =(Wt_sum)/Cnum;
Wt = tmps - mean(dg)*Wt;
Wt = Wt/norm(Wt);
%FAST ICA计算公式
%正交化检测
if issparse(Register_Matrix) == 0
Rs = abs(qr(Register_Matrix));
else
Rs = Register_Matrix*Register_Matrix';
end
Rs(1,2) = 0;
Rs(2,1) = 0;
RR = Rs;
%正交化
Wt = Wt - RR*Wt;
Wt = Wt/norm(Wt);
if abs(abs(Wt'*W_n_1)-1) < global_error
Register_Matrix(:,ker) = Wt;
break;
else
Register_Matrix(:,ker) = Register_Matrix(:,ker);
end
if norm(Wt - W_n_1) < global_error
Register_Matrix(:,ker) = Wt;
break;
else
Register_Matrix(:,ker) = Register_Matrix(:,ker);
end
index = index+1;
Error = [Error,abs(abs(Wt'*W_n_1)-1)];
W_sum = [W_sum,sum(sum(Wt))];
%计算ICA后的矩阵
Y_final = Register_Matrix'*Q*y0;
Y_final = Y_final/max(max(abs(Y_final)));
Err = [Err,mean(Y_final(2,:)-data)];
end
Err2{ker} = Err;
Error2{ker} = Error;
W_sum2{ker} = W_sum;
end
figure;
plot(Error2{1},'b-o','linewidth',2);
xlabel('迭代次数');
ylabel('收敛误差');
grid on
%计算ICA后的矩阵
Y_final = Register_Matrix'*Q*y0;
Register_Matrix
Y_final = Y_final/max(max(abs(Y_final)));
%将混合矩阵重新排列并输出
figure;
subplot(311);
plot(Y_final(1,:));
title('伪迹');
subplot(312);
plot(Y_final(2,:));
title('EEG去伪迹信号');
subplot(313);
plot(data,'b');
hold on
plot(Y_final(2,:),'r');
legend('原始信号','去伪迹后信号');
N=1024;
wn = Y_final(1,:);
Pxx=10*log10(abs(fft(wn).^2)/N);
f=(0:length(Pxx)-1)/length(Pxx)*1e2;
figure;
plot(f(1:floor(length(f)/2)),Pxx(1:floor(length(f)/2)));
xlabel('频率');
ylabel('功率（dB）');
title('周期图法N=256')
D1 = Pxx(1:floor(length(f)/2));
dif1= std(D1)
wn = Y_final(2,:);
Pxx=10*log10(abs(fft(wn).^2)/N);
f=(0:length(Pxx)-1)/length(Pxx)*1e2;
figure;
plot(f(1:floor(length(f)/2)),Pxx(1:floor(length(f)/2)));
xlabel('频率');
ylabel('功率（dB）');
title('周期图法N=256')
D2 = Pxx(1:floor(length(f)/2));
dif2= std(D2)
4.操作步骤与仿真结论
5.参考文献
[1] D Maino, Fa Rusi A , Baccigalupi C , et al. All-sky astrophysical component separation with Fast Independent Component Analysis (FastICA)[J]. Monthly Notices of the Royal Astronomical Society, 2010(1):53-68.
A28-50
