功率谱密度代码

matlab脑电信号处理
t=0.001:0.001:1;
x=load('C:\Users\yxzhang\Desktop\rest_close.txt'); %读取文件
y=load('C:\Users\yxzhang\Desktop\audio_close.txt');

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