ecg_process_buttord_bandstop_filter.asv

计算脑电信号非线性参数

%本程序采用巴特沃什带通滤波器%
%林宛华  中南大学生物医学工程医学图像与信号处理实验室   %
%E-mail: [e-mail]%
%2008年3月10日%
function ECG_process_buttord_bandstop_filter()
clear
clc
load ECG1; 
data2 = data2';
t=0:0.001:0.999;
%data2 = zeros(1,1000);
%t=0:0.001:1.3;
%x=sin(2*pi*50*t)+sin(2*pi*300*t);
%f=x+3.5*randn(1,length(t));
noise = sin(2*pi*50*t);
data_noise = data2 + noise;
subplot(2,1,1),plot(t,data2);Ylabel = ('幅值');Xlabel = ('时间');
y = fft(data2,1000);
p = y.*conj(y)/1000;  %计算功率谱密度
ff = 1000*(0:499)/1000;    %计算变换后不同点所对应的频率值
subplot(2,1,2),plot(ff,p(1:500));
figure,subplot(2,1,1),plot(t,data_noise);Ylabel = ('幅值');Xlabel = ('时间');
y1 = fft(data_noise,1000);
p1 = y1.*conj(y1)/1000;  %计算功率谱密度
ff1 = 1000*(0:499)/1000;    %计算变换后不同点所对应的频率值
subplot(2,1,2),plot(ff1,p1(1:500));

w1=45/500;
w2=55/500;
[B,A]=butter(1,[w1,w2],'stop');

s=filter(B,A,data_noise);
figure,subplot(2,1,1),plot(t,s);Ylabel = ('幅值');Xlabel = ('时间');
y11 = fft(s,1000);
p11 = y11.*conj(y11)/1000;  %计算功率谱密度
ff11 = 1000*(0:499)/1000;    %计算变换后不同点所对应的频率值
subplot(2,1,2),plot(ff11,p11(1:500));
for m=0:4;
   % close;
    s=filter(B,A,s);
    figure,subplot(2,1,1),plot(t,s);Ylabel = ('幅值');Xlabel = ('时间');
    y11 = fft(s,1000);
    p11 = y11.*conj(y11)/1000;  %计算功率谱密度
    ff11 = 1000*(0:499)/1000;    %计算变换后不同点所对应的频率值
    subplot(2,1,2),plot(ff11,p11(1:500));
end