%本程序采用巴特沃什帶通濾波器實現去心電信號工頻干擾%
%林宛華  中南大學生物醫學工程醫學圖像與信號處理實驗室   %
%E-mail: wanhualin01@163.com%
%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