q_ica_filter.m

用于脑电分析的各类源码

function data=q_ica_filter(data,trail_t,fs,frq,kept_c_num)
% 命令格式为：data=q_ica_filter(data,trail_t,fs,frq,kept_c_num)
% 使用ica对多导数据分解后只选取frq频率成份最高的kept_c_num个分量重建信号
% 输入参数： data，输入数据，每行为一导信号
%           trail_t,单次trail的持续时间
%           fs,采样频率
%           frq，保留频段，格式为一维数组
%           kept_c_num，保留分量个数
% 输出参数： data，重构信号
%
% 注：      输入数据为n个trail顺序连接而成。

l=length(data);
chan_num=size(data,1);
trail_n=floor(l/trail_t/fs);

data_new=[];
for i=1:trail_n
    i% one trail
    temp_index=[1:trail_t*fs]+(i-1)*trail_t*fs;
    data_temp=data(:,temp_index);
    % ica process
    [ww,ss]=runica(data_temp,'extended',-1);
    ica_temp=ww*ss*data_temp;
    
    % get rid of components with low 10hz psd
    for j=1:chan_num
        psd_temp(j,:)=q_psfd_cal(ica_temp(j,:),frq,fs);
    end
    psd_amp=mean(abs(psd_temp),2);  % psd mean value of each component in frq
    [a order_ind]=sort(psd_amp);    % sort from lowest to highest
    min_order_ind=order_ind(1:chan_num-kept_c_num);   % low components no
    ica_temp(min_order_ind,:)=0*ica_temp(min_order_ind,:); % kept high psd components
    
    % data reconstuction 
    data_temp=inv(ww*ss)*ica_temp;
    data_new=[data_new,data_temp];
end

data=data_new;