驾驶员清醒和疲劳状态时 如何提取脑电信号特征？

编写一程序。。驾驶员清醒和疲劳状态时脑电信号特征提取，只需要把脑电信号转换成功率谱

脑电信号已采集

(资料图片)

t=0.001:0.001:1;

x=load("D:\EEGdata_data");     %读取文件

y=load("D:EEGdata_newdata\");

xx={};           %每个导联的数据存储

yy={};

n=1000;           %数据数目

sc=7;             %小波包的分解尺度

for i=1:1:8       %导联的数据分离

xx{i}=x(:,i);

yy{i}=y(:,i);

end

for i=1:1:8

%画出原始信号图像

figure

subplot(2,2,1)

plot(t,xx{i})

axis([min(t) max(t) 1.1*floor(min(xx{i})) 1.1*ceil(max(xx{i}))])

title("rest close 原始信号")

ylabel("幅值")

subplot(2,2,2)

plot(t,yy{i})

axis([min(t) max(t) 1.1*floor(min(yy{i})) 1.1*ceil(max(yy{i}))])

title("audio close 原始信号")

ylabel("幅值")

%fft_原始信号的频谱分析

xx1=fft(xx{i},n);

pxx1=xx1.*conj(xx1)/n;

yy1=fft(yy{i},n);

pyy1=yy1.*conj(yy1)/n;

%画出0-30hz内的功率谱图像

n=60;

f=1:n/2;

subplot(2,2,3);

plot(f,pxx1(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("rest_close功率谱")

subplot(2,2,4);

plot(f,pyy1(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("audio_close的功率谱")

%分解信号选择分解尺度为7，同时重构信号

wpt=wpdec(xx{i},sc,"db7","shannon");  %小波包分解信号

xx80=wprcoef(wpt,[sc,0]);              %重构信号

xx81=wprcoef(wpt,[sc,1]);

xx82=wprcoef(wpt,[sc,2]);

xx83=wprcoef(wpt,[sc,3]);

wpt=wpdec(yy{i},8,"db7","shannon");

yy80=wprcoef(wpt,[sc,0]);

yy81=wprcoef(wpt,[sc,1]);

yy82=wprcoef(wpt,[sc,2]);

yy83=wprcoef(wpt,[sc,3]);

%画出重构信号

figure

subplot(2,1,2);plot(yy80);

title("audio close delta");

ylabel("幅值");

subplot(2,1,1);plot(xx80);

title("rest close delta");

ylabel("幅值");

figure

subplot(2,1,1);plot(xx81);

title("rest close theta");

ylabel("幅值");

subplot(2,1,2);plot(yy81);

title("audio close theta");

ylabel("幅值");

figure

subplot(2,1,1);plot(xx82);

title("rest close alpha");

ylabel("幅值");

subplot(2,1,2);plot(yy82);

title("audio close alphta");

ylabel("幅值");

figure

subplot(2,1,1);plot(xx83);

title("rest close beta");

ylabel("幅值")

subplot(2,1,2);plot(yy83);

title("audio close beta");

ylabel("幅值");

n=1000;

%fft_重构信号的频谱分析

xx180=fft(xx80,n);

pxx180=xx180.*conj(xx180)/n;

xx181=fft(xx81,n);

pxx181=xx181.*conj(xx181)/n;

xx182=fft(xx82,n);

pxx182=xx182.*conj(xx182)/n;

xx183=fft(xx83,n);

pxx183=xx183.*conj(xx183)/n;

yy180=fft(yy80,n);

pyy180=yy180.*conj(yy180)/n;

yy181=fft(yy81,n);

pyy181=yy181.*conj(yy181)/n;

yy182=fft(yy82,n);

pyy182=yy182.*conj(yy182)/n;

yy183=fft(yy83,n);

pyy183=yy183.*conj(yy183)/n;

%画出重构信号的功率谱图

n=60;

f=1:n/2;

figure

subplot(2,1,1);

plot(f,pyy180(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("rest close delta的功率谱");

subplot(2,1,2);

plot(f,pyy180(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("audio close detta的功率谱");

figure

subplot(2,1,1);

plot(f,pxx181(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("rest close theta的功率谱");

subplot(2,1,2);

plot(f,pyy181(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("audio close theta的功率谱");

figure

subplot(2,1,1);

plot(f,pxx182(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("rest close alpha的功率谱");

subplot(2,1,2);

plot(f,pyy182(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("audio close alpha的功率谱");

figure

subplot(2,1,1);

plot(f,pxx183(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("rest close beta的功率谱");

subplot(2,1,2);

plot(f,pyy183(1:n/2));

ylabel("功率谱幅值(mv^2)");

title("audio close beta的功率谱");

end;

运行不粗来，得不到功率谱图.........TAT,求帮助！