用matlab分析脑电,在matlab中如何做功率谱分析

%function zhao A1=load('E:\biology\brain\BCI\file\toyekun\zhaochunfeng20080319.dat'); C3=A1(:,1); Cz=A1(:,2); C4=A1(:,3); B1=load('E:\biology\brain\BCI\file\toyekun\zhaochunfeng20080319.txt'); dist=B1(:,2); index=B1(:,6); %clear B; depend=length(C3); %B1=zeros(160,9000); %采样率1000Hz，trial长度9s %D1=B1; j=0; k=0; for i=1:160 if dist(i)<1.5 j=j+1; LC3(j,:)=C3(index(i):index(i)+8999); % 1599是根据采样率来决定的，200*8 LCz(j,:)=Cz(index(i):index(i)+8999); LC4(j,:)=C4(index(i):index(i)+8999); else k=k+1; RC3(k,:)=C3(index(i):index(i)+8999); RCz(k,:)=Cz(index(i):index(i)+8999); RC4(k,:)=C4(index(i):index(i)+8999); end end ALC3=mean(LC3(1:80,:)); ARC3=mean(RC3(1:80,:)); ALCz=mean(LCz(1:80,:)); ARCz=mean(RCz(1:80,:)); ALC4=mean(LC4(1:80,:)); ARC4=mean(RC4(1:80,:)); %while thinking nothing A2=load('E:\biology\brain\BCI\file\toyekun\zhaochunfeng20080319_2.dat'); c3=A2(:,1); cz=A2(:,2); c4=A2(:,3); B2=load('E:\biology\brain\BCI\file\toyekun\zhaochunfeng20080319_2.txt'); dist2=B2(:,2); index2=B2(:,6); %clear B; depend=length(c3); %B2=zeros(160,9000); %采样率1000Hz，trial长度9s %D2=B2; for i=1:40 nc3(i,:)=c3(index2(i):index2(i)+8999); % 1599是根据采样率来决定的，200*8 ncz(i,:)=cz(index2(i):index2(i)+8999); nc4(i,:)=c4(index2(i):index2(i)+8999); end c3=mean(nc3(1:40,:)); cz=mean(ncz(1:40,:)); c4=mean(nc4(1:40,:)); % All frequency values are in Hz. Fs = 1000; % Sampling Frequency Fstop1 = 7; % First Stopband Frequency Fpass1 = 8; % First Passband Frequency Fpass2 = 65; % Second Passband Frequency Fstop2 = 67; % Second Stopband Frequency Astop1 = 60; % First Stopband Attenuation (dB) Apass = 1; % Passband Ripple (dB) Astop2 = 80; % Second Stopband Attenuation (dB) match = 'stopband'; % Band to match exactly % Construct an FDESIGN object and call its BUTTER method. h = fdesign.bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, ... Astop2, Fs); Hd = design(h, 'butter', 'MatchExactly', match); ALC3=filter(Hd,ALC3); ALCz=filter(Hd,ALCz); ALC4=filter(Hd,ALC4); ARC3=filter(Hd,ARC3); ARCz=filter(Hd,ARCz); ARC4=filter(Hd,ARC4); c3=filter(Hd,c3); cz=filter(Hd,cz); c4=filter(Hd,c4); figure subplot(2,2,1); spectrogram(ALC3,500,250,512,1000); title('left hand C3'); %figure subplot(2,2,2); spectrogram(ALC4,500,250,512,1000); title('left hand C4'); %figure subplot(2,2,3); spectrogram(ARC3,500,250,512,1000); title('right hand C3'); %figure subplot(2,2,4); spectrogram(ARC4,500,250,512,1000); title('right hand C4'); figure subplot(3,1,1); spectrogram(c3,500,250,512,1000); title('C3'); %figure subplot(3,1,2); spectrogram(cz,500,250,512,1000); title('Cz'); %figure subplot(3,1,3); spectrogram(c4,500,250,512,1000); title('C4');