【腦電訊號】基於matlab小波睡眠監測【含Matlab原始碼 595期】
阿新 • • 發佈:2021-06-25
一、簡介
基於matlab的腦電波睡眠監測
二、原始碼
data0=rand(1,9999); %腦電訊號原始資料 tm=0.02; %取樣時間間隔 td=1:tm:30; %取時間1-30秒 data=data0(1:(30-1)/tm+1); %1-30秒的資料 figure(1) subplot(211); plot(td,data); xlabel('時間(秒)'),ylabel('腦波電壓'),title('1-30s腦電圖波'); %fft變換 Fs=1000; n=length(data); data1=fft(data,n);%傅立葉變換% df=Fs/length(data1); %頻域解析度 Fx=df*(0:length(data1)-1); %將橫軸變為頻率軸% figure(1); subplot(212); plot(Fx,abs(data1)); %繪製腦電波訊號的頻譜圖% axis([0 150 0 60]); title('頻譜圖');xlabel('頻率/Hz');ylabel('幅值'); %訊號預處理 %低通濾波——巴特沃斯濾波器 Fs=1000; fp=30; fs=40; Ap=1; As=30; Wp=fp/(Fs/2);%計算歸一化角頻率 Ws=fs/(Fs/2); [N,Wc]=buttord(Wp,Ws,Ap,As);%計算階數和截止頻率 [b,a]=butter(N,Wc,'low');%計算H(z)分子、分母多項式係數 [H,F]=freqz(b,a,500,Fs);%計算H(z)的幅頻響應,freqz(b,a,計算點數,取樣速率) figure(2) subplot(2,2,2) ; plot(F,20*log10(abs(H))) ; xlabel('頻率(Hz)'); ylabel('幅度(dB)') axis([0 100 -30 3]); grid on ; subplot(2,2,1) plot(F,abs(H)); xlabel('頻率(Hz)'); ylabel('幅度 ') ; title('低通濾波器'); axis([0 100 0 2]); grid on; subplot(2,2,3); pha=angle(H)*180/pi; plot(F,pha); xlabel('頻率(Hz)'); ylabel('相位(dB)') axis([0 100 -200 200]); grid on; % % %用低通巴特沃斯濾波器 Q=filter(b,a,data); figure(3) subplot(211); plot(td,Q); title('過巴斯後時域圖');xlabel('時間');ylabel('幅值'); %fft變換 n=length(Q); Q1=fft(Q,n);%傅立葉變換% df=Fs/length(Q1); %頻域解析度 Fx=df*(0:length(Q1)-1); %將橫軸變為頻率軸% figure(3); subplot(212); plot(Fx,abs(Q1)); %繪製腦電波訊號的頻譜圖% axis([0 50 0 60]); title('頻譜圖');xlabel('頻率/Hz');ylabel('幅值'); % %------小波閾值去噪 %------軟閾值小波去噪 [c,s]=wavedec2(Q,2,'db5'); [thr,sorh,keepapp] = ddencmp('den','wv',Q); [xc,cxc,lxc,perf0,perfl2]=wdencmp('gbl',Q,'sym4',2,thr,sorh,keepapp);% figure(4); subplot(211); plot(td,xc); title('小波去噪後時域圖');xlabel('時間');ylabel('幅值'); u=xc; %fft變換 n=length(xc); xc1=fft(xc,n);%傅立葉變換% df=Fs/length(xc1); %頻域解析度 Fx=df*(0:length(xc1)-1); %將橫軸變為頻率軸% figure(4); subplot(212); plot(Fx,abs(xc1)); %繪製聲音訊號的頻譜圖% axis([0 60 0 100]); title('小波去噪後頻譜圖');xlabel('頻率/Hz');ylabel('幅值'); %-- 第一段波過帶通濾波器------------ n=length(xc); fs=[0.1,5];fp=[0.5,3]; fo=1200; %取樣頻率 wp=2.*fp./fo;ws=2.*fs./fo; rp=1;as=40; f = design(fdesign.bandpass(fs(1),fp(1),fp(2),fs(2),as,rp,as,fo),'butter'); y1=filter(f,xc); yt1=fft(y1,n); df=Fs/length(yt1); %頻域解析度 Fx=df*(0:length(yt1)-1); %將橫軸變為頻率軸% figure(5); subplot(3,1,1),plot(td,y1); xlabel('時間(秒)'),ylabel('腦波電壓'),title('1-30s第一階段腦電圖波'); subplot(3,1,2); plot(Fx,abs(yt1)); %繪訊號的頻譜圖% axis([0 6 0 100]); title('第一階段頻譜圖');xlabel('頻率/Hz');ylabel('幅值'); %功率譜 fs=800;ts=1/fs; t=0:ts:2; nfft=64; power1=(norm(y1)^2/length(y1+1)); spow1=abs(fft(y1,nfft).^2); f=(0:nfft-1)/ts/nfft; f=f-fs/2; figure(5); subplot(3,1,3); plot(f,fftshift(spow1),'k'); title('第一階段功率譜圖');xlabel('頻率');ylabel('功率譜'); disp(['power1=',num2str(power1),'.']); %-- 第二段波過帶通濾波器------------ n=length(u); fs=[3,9];fp=[4,7]; fo=1200; %取樣頻率 wp=2.*fp./fo;ws=2.*fs./fo; rp=1;as=40; f = design(fdesign.bandpass(fs(1),fp(1),fp(2),fs(2),as,rp,as,fo),'butter'); y2=filter(f,u); yt2=fft(y2,n); df=Fs/length(yt2); %頻域解析度 Fx=df*(0:length(yt2)-1); %將橫軸變為頻率軸% figure(6); subplot(3,1,1); plot(td,y2); xlabel('時間(秒)');ylabel('腦波電壓');title('1-30s第二階段腦電圖波'); subplot(3,1,2); plot(Fx,abs(yt2)); %繪訊號的頻譜圖% axis([2 15 0 50]); title('第二階段頻譜圖');xlabel('頻率/Hz');ylabel('幅值'); %功率譜 fs=800;ts=1/fs; t=0:ts:2; nfft=64; power2=(norm(y2)^2/length(y2+1)); spow2=abs(fft(y2,nfft).^2); f=(0:nfft-1)/ts/nfft; f=f-fs/2; figure(6); subplot(3,1,3); plot(f,fftshift(spow2),'k'); title('第一階段功率譜圖');xlabel('頻率');ylabel('功率譜'); disp(['power2=',num2str(power2),'.']); %-- 第三段波過帶通濾波器------------ n=length(u); fs=[7,14];fp=[8,13]; fo=1200; %取樣頻率 wp=2.*fp./fo;ws=2.*fs./fo; f = design(fdesign.bandpass(fs(1),fp(1),fp(2),fs(2),as,rp,as,fo),'butter'); y3=filter(f,u); df=Fs/length(yt3); %頻域解析度 Fx=df*(0:length(yt3)-1); %將橫軸變為頻率軸% figure(7); subplot(3,1,1); plot(td,y3); xlabel('時間(秒)');ylabel('腦波電壓');title('1-30s第三階段腦電圖波'); subplot(3,1,2); plot(Fx,abs(yt3)); %繪訊號的頻譜圖% axis([4 16 0 30]); title('第三階段頻譜圖');xlabel('頻率/Hz');ylabel('幅值'); %功率譜 fs=800;ts=1/fs; t=0:ts:2; nfft=64; power3=(norm(y3)^2/length(y3+1)); spow3=abs(fft(y3,nfft).^2); f=(0:nfft-1)/ts/nfft; f=f-fs/2; figure(7); subplot(3,1,3); plot(f,fftshift(spow3),'k'); title('第三階段功率譜圖');xlabel('頻率');ylabel('功率譜'); disp(['power3=',num2str(power3),'.']);
三、執行結果
四、備註
版本:2014a
完整程式碼或代寫加1564658423