mu.zip_能量 脑电_能量矩_脑电 能量_脑电能量计算_运动想象mu

close all clear all %=====================================读取脑电数据,并显示波形 load('F:\2003\1\dataset3\graz_data\dataset_BCIcomp1.mat'); count=135; %组数 y_train(count) c3=x_train(:,1,count); c4=x_train(:,3,count); figure(1); subplot(2,1,1) plot((1:1152)/128,c3); xlabel('时间'); ylabel('C3通道脑电信号幅度'); subplot(2,1,2) plot((1:1152)/128,c4); xlabel('时间'); ylabel('C4通道脑电信号幅度'); l=1152; fs=128; %==================滤波器设计 wn=[0.125 0.1875]; %8~12Hz带通滤波 wn1=[0.28 0.375]; %18~24Hz带通滤波 b = fir1(50,wn,'DC-0'); figure(2) freqz(b,1,64,fs); %==================信号滤波 l1=filter(b,1,c3); l2=filter(b,1,c4); nf=2^nextpow2(1152); k=fs*linspace(0,1,nf)-fs/2; %归一化 L1=fftshift(fft(l1,nf)); L2=fftshift(fft(l2,nf)); figure(3); plot(k,abs(L1)); hold on; plot(k,abs(L2),'r'); figure(4); subplot(2,1,1) plot((1:1152)/128,l1); xlabel('时间'); ylabel('滤波后C3通道信号幅度'); subplot(2,1,2) plot((1:1152)/128,l2); xlabel('时间'); ylabel('滤波后C4通道信号幅度'); %======================================计算基于mu节律二阶矩能量，固定窗长512，每次移动1 %==========数据不到512，变窗口长度，每次移动1 m1(1)=l1(1).^2; m2(1)=l2(1).^2; for i=1:511 m1(i+1)=i*m1(i)/(i+1)+(l1(i+1).^2)/(i+1); m2(i+1)=i*m2(i)/(i+1)+(l2(i+1).^2)/(i+1); end %==========数据超过512，固定窗口长度512，每次移动1 for j=512:1151 m1(j+1)=m1(j)-(l1(j-511).^2-l1(j+1).^2)/512; m2(j+1)=m2(j)-(l2(j-511).^2-l2(j+1).^2)/512; end figure(5) plot((1:1152)/128,m1(1:1152)); hold on; plot((1:1152)/128,m2(1:1152),'r'); % axis([0 9 0 0.03]) xlabel('时间'); ylabel('基于mu节律的能量（二阶矩）'); title('EEG(C3,C4)信号的在线二阶矩时序图') % %============================从第3秒开始显示 % figure(6) % plot((384:1152)/128,m1(384:1152)); % hold on; % plot((384:1152)/128,m2(384:1152),'r'); % % axis([3 9 0 0.03]) % xlabel('时间'); % ylabel('基于mu节律的能量（二阶矩）'); % title('EEG(C3,C4)信号的在线二阶矩时序图') % % %================结果 C3-C4>0 左手 C3-C4<0 右手 % for n=384:1152 % sub(n)=m1(n)-m2(n); % if(sub(n)>0) % result(n)=1; % else % result(n)=-1; % end % end % figure(7) % plot((384:1152)/128,result(384:1152)); % axis([3 9 -2 2]) % %