matlab計算相對功率

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