• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

matlab计算相对功率

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

1、对脑电数据进行db4四层分解,因为脑电频率是在0-64HZ,分层后如图所示,

细节分量[d1 d2 d3 d4]

近似分量[a4]

重建细节分量和近似分量,然后计算对应频段得相对功率谱,重建出来得四个频段(αβθδ)都有14个通道,所以要计算4频段14通道共56个相对功率

 

 2、代码

function wavelet(signal)
A4Array = zeros(14,5000);
D4Array = zeros(14,5000);
D3Array = zeros(14,5000);
D2Array = zeros(14,5000);
  for i=1:14
   [C,L] = wavedec(signal(i,1:5000),4,'db4');%函数返回 3 层分解的各组分系数C(连接在一个向量里) ,向量 L 里返回的是各组分的长度。
   % [cD1,cD2,cD3,cD4] = detcoef(C,L,[1,2,3,4]);%抽取1234层细节系数
   % cA4 = appcoef(C,L,'d4',4);%抽取近似系数
   A4 = wrcoef('a',C,L,'db4',4);%重建4层近似,deta波
   A4Array(i,:) = A4;
   D4 = wrcoef('d',C,L,'db4',4);%重建4层细节,sita波
   D4Array(i,:) = D4;
   D3 = wrcoef('d',C,L,'db4',3);%重建3层细节,alpha波
   D3Array(i,:) = D3;
   D2 = wrcoef('d',C,L,'db4',2);%重建2层细节,beta波
   D2Array(i,:) = D2;
  end
     detaspectral(signal,A4Array);
     thetaspectral(signal,D4Array);
     alphaspectral(signal,D3Array);
     betaspectral(signal,D2Array);
end

  

detaspectral thetaspectral alphaspectral betaspectra的代码都是一样的
function alphaspectral(signal,dtest8theta)
  Fs=128;
  N=1024;Nfft=256;n=0:N-1;t=n/Fs;
  window=hanning(256);
  noverlap=128;
  dflag='none';
  for i=1:14
    x=signal(i,1:5000);
    powd(i,:)=psd(x,Nfft,Fs,window,noverlap,dflag);%计算未分频段,总数据的功率谱
    x1=dtest8theta(i,:);%某一频段的脑电数据
    powd1(i,:)=psd(x1,Nfft,Fs,window,noverlap,dflag);%计算该频段的功率谱
  end
  xdpowthetad = zeros(14,1);
  xdpowthetad=mean(abs(powd1),2)./mean(abs(powd),2);%计算相对功率,用分频段功率谱比上不分频段的。
  %save('G:\研三\音乐反馈数据\新算相对功率\xdpowthetad.mat','xdpowthetad');
  save('C:\Users\25626\Desktop\滤波后数据\14\相对功率谱\5 3\alphaspectra.mat','xdpowthetad');
end

  

function detaspectral(signal,dtest8theta)
  Fs=128;
  N=1024;Nfft=256;n=0:N-1;t=n/Fs;
  window=hanning(256);
  noverlap=128;
  dflag='none';
  for i=1:14
    x=signal(i,1:5000);
    powd(i,:)=psd(x,Nfft,Fs,window,noverlap,dflag);%计算未分频段,总数据的功率谱
    x1=dtest8theta(i,:);%某一频段的脑电数据
    powd1(i,:)=psd(x1,Nfft,Fs,window,noverlap,dflag);%计算该频段的功率谱
  end
  xdpowthetad = zeros(14,1);
  xdpowthetad=mean(abs(powd1),2)./mean(abs(powd),2);%计算相对功率,用分频段功率谱比上不分频段的。
  %save('G:\研三\音乐反馈数据\新算相对功率\xdpowthetad.mat','xdpowthetad');
  save('C:\Users\25626\Desktop\滤波后数据\14\相对功率谱\5 3\detaspectral.mat','xdpowthetad');
end

  

function betaspectral(signal,dtest8theta)
  Fs=128;
  N=1024;Nfft=256;n=0:N-1;t=n/Fs;
  window=hanning(256);
  noverlap=128;
  dflag='none';
  for i=1:14
    x=signal(i,1:5000);
    powd(i,:)=psd(x,Nfft,Fs,window,noverlap,dflag);%计算未分频段,总数据的功率谱
    x1=dtest8theta(i,:);%某一频段的脑电数据
    powd1(i,:)=psd(x1,Nfft,Fs,window,noverlap,dflag);%计算该频段的功率谱
  end
  xdpowthetad = zeros(14,1);
  xdpowthetad=mean(abs(powd1),2)./mean(abs(powd),2);%计算相对功率,用分频段功率谱比上不分频段的。
  %save('G:\研三\音乐反馈数据\新算相对功率\xdpowthetad.mat','xdpowthetad');
  save('C:\Users\25626\Desktop\滤波后数据\14\相对功率谱\5 3\betaspectral.mat','xdpowthetad');
end

  

function thetaspectral(signal,dtest8theta)
  Fs=128;
  N=1024;Nfft=256;n=0:N-1;t=n/Fs;
  window=hanning(256);
  noverlap=128;
  dflag='none';
  for i=1:14
    x=signal(i,1:5000);
    powd(i,:)=psd(x,Nfft,Fs,window,noverlap,dflag);%计算未分频段,总数据的功率谱
    x1=dtest8theta(i,:);%某一频段的脑电数据
    powd1(i,:)=psd(x1,Nfft,Fs,window,noverlap,dflag);%计算该频段的功率谱
  end
  xdpowthetad = zeros(14,1);
  xdpowthetad=mean(abs(powd1),2)./mean(abs(powd),2);%计算相对功率,用分频段功率谱比上不分频段的。
  %save('G:\研三\音乐反馈数据\新算相对功率\xdpowthetad.mat','xdpowthetad');
  save('C:\Users\25626\Desktop\滤波后数据\14\相对功率谱\5 3\thetaspectral.mat','xdpowthetad');
end

  


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
Delphi通过查找字符定位TADOQuery数据的位置发布时间:2022-07-18
下一篇:
数据类型表(DELPHI、C++)发布时间:2022-07-18
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap