matlab计算相对功率
1、对脑电数据进行db4四层分解,因为脑电频率是在0-64HZ,分层后如图所示,
细节分量[d1 d2 d3 d4]
近似分量[a4]
重建细节分量和近似分量,然后计算对应频段得相对功率谱,重建出来得四个频段(αβθδ)都有14个通道,所以要计算4频段14通道共56个相对功率
2、代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | 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 |
1 | detaspectral thetaspectral alphaspectral betaspectra的代码都是一样的 |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | 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 |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | 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 |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | 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 |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | 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 |
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步