MATLAB-脑电信号处理程序

基于matlab的脑电信号处理 程序设计，包括GUI界面设计

S1=‘E:\医学信号处理课程设计\normal1.txt‘;
Y1=load（S1）;
Y=Y1‘;
x=Y（1:）;y=Y（2:）;
n=length（x）;%求取样点数
T=0.004;fs=1/T;
fk=（[0:n-1]/n/T）;%横坐标变换
Y1=fft（y）;Y1r=abs（Y1）;A1r=angle（Y1）;%进行fft变换并求出频谱图
figure （1）%显示时频域
subplot 311;plot（xy）;title（‘脑电信号时域图‘）;xlabel（‘时间/s‘）;ylabel（‘幅值/uV‘）;
subplot 312;plot（fk20*log10（Y1r））;    title（‘脑电信号频域图‘）;xlabel（‘频率/Hz‘）;ylabel（‘20*log10（H）/dB‘）;axis（[0 125 0 100]）; 
subplot 313; plot（fkA1r）; axis（[0 125 -4 4]）;title（‘脑电信号相位图‘）;xlabel（‘频率/Hz‘）;
%矩形窗
wp=0.232*pi;ws=0.248*pi;
deltaw=ws-wp;                  
N=ceil（1.8*pi/deltaw）%定义滤波器的阶数 
wc=30/（fs/2）;
H1=fir1（Nwcstepseq（N+1））;%使用fir1函数设计滤波器
[h1w]=freqz（H11512250）;
 sf=filter（H11y）;
 Fsf=fft（sf）;
figure （2）%显示FIR滤波器
subplot 211;plot（w20*log10（abs（h1）））;title（‘低通滤波器模型‘）;xlabel（‘频率/Hz‘）; ylabel（‘20*log10（H）/dB‘）;axis（[0 125 -100 10]）;%参数分别是频率与幅值
subplot 212;plot（wangle（h1））;
figure（3）%低通滤波后的信号图像
subplot（211）;plot（xsf）;title（‘低通滤波后的信号时域图‘）;xlabel（‘时间/s‘）;ylabel（‘幅值/uV‘）;
subplot（212）;plot（fk（abs（Fsf）））;axis（[0 125 0 70000]）; title（‘低通滤波后的信号频域图‘）;xlabel（‘频率/Hz‘）; %ylabel（‘20*log10（H）/dB‘）;
%巴特沃斯
wp0=[4 8];ws0=[2 10];Rp=1;  Rs=18;   %设置通带波纹系数
wp1=wp0*pi*2/fs;               %设置通带频率
ws1=ws0*pi*2/fs;                %设置阻带频率
wp=2*fs*tan（wp1./2）;
ws=2*fs*tan（ws1./2）;
[NWn]=buttord（wpwsRpRs‘s‘）       %求巴特沃斯滤波器阶数，输出参数N代表满足设计要求的滤波器的最小阶数，Wn是等效低通滤波器的截止频率
[zpk]=buttap（N）;
[bpap]=zp2tf（zpk）;
Bp=wp（2）-wp（1）;
Wp0=sqrt（wp（1）*wp（2））;
[bsas]=lp2bp（bpapWp0Bp）;
[bzaz]=bilinear（bsasfs）;
[Hbwb]=freqz（bzaz125fs）;%125个数据点
h=（filter（bzazy））;%原信号时域滤波
H=fft（h）;
figure （4）%显示IIR滤波器
subplot 211;plot（wb20*log10（abs（Hb）））;title（‘带通滤波器模型‘）;xlabel（‘频率/Hz‘）; ylabel（‘20*log10（H）/dB‘）;
subplot 212; plot（wbangle（Hb））; title（‘脑电信号相位图‘）;xlabel（‘频率/Hz‘）;%axis（[0 125 -4 4]）;
figure （5）%带通滤波前后的信号频域
subplot 211;plot（xy）;title（‘带通滤波前时域‘）;xlabel（‘时间/s‘）;ylabel（‘幅值/uV‘）;axis（[0 60 -400 400]）;
subplot 212;plot（xh）;title（‘带通滤波后时域‘）;xlabel（‘时间/s‘）;ylabel（‘幅值/uV‘）;axis（[0 60 -400 400]）;
figure （6）%带通滤波前后的信号时域
subplot 211;plot（fkY1r）;title（‘带通滤波前频域‘）;xlabel（‘频率/Hz‘）;axis（[0 125 0 70000]）;
subplot 212;plot（fkabs（H））;title（‘带通滤波后频域‘）;xlabel（‘频率/Hz‘）;axis（[0 125 0 40000]）;
%陷波器
[mn]=findpeaks（abs（Fsf）‘minpeakdistance‘270）;
for i=2:4;
    a1=n（i）/61.44;
    rp=1;rs=60;Fs=250;
    wp=[a1-0.4 a1+0.4]*2/Fs;ws=[a1 a1+0.2]*2/Fs;
    [NWn]=ellipord（wpwsrprs）;
    [ba]=ellip（NrprsWn‘stop‘）;
    [Hf]=freqz（ba512Fs）;
    Y=filter（basf）;
    h2=fft（Y）;
    sf=Y;
end
    figure （7）
    subplot 211;plot（xY）;title（‘陷波滤波后时域‘）;xlabel（‘时间/s‘）;ylabel（‘幅值/uV‘）;axis（[0 60 -400 400]）;
    subplot 212;plot（fk（abs（h2）））;title（‘陷波滤波后频域‘）;xlabel（‘频率/Hz‘）;axis（[0 125 0 70000]）;
%频谱图
    Hs=spectrum.yulear;  %周期图法
    figure （8）;
    [Pxxf]= pwelch（Y1r33[]15360Fs）;%33为信号截取长度
    plot（fPxx）;axis（[012504000000]）;xlabel（‘频率/Hz‘）;title（‘频谱图‘）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3271  2019-01-16 09:28  all_figure.m

     文件        370  2019-01-01 19:54  box.m

     文件       1655  2019-01-16 09:06  butterfilter.m

     文件        460  2019-01-14 10:44  eeg_fft.m

     文件      17479  2019-01-16 08:44  eeg_filter.m

     文件        914  2019-01-07 11:10  eeg_fir.m

     文件        198  2019-01-12 14:14  eeg_load.m

     文件        747  2019-01-14 10:33  eeg_power.m

     文件        133  2019-01-01 19:47  ideallp.m

     文件         77  2019-01-12 13:21  stepseq.m

----------- ---------  ---------- -----  ----

                25304                    10