fid=fopen（‘109a.txt‘）;
C=textscan（fid‘%8c?%f?%*f‘‘headerlines‘2）;
fclose（fid）;
a=C{1};
M=C{2};
?k=length（a）
for?i=1:k
???c（i）=strread（a（i:）‘%*s?%f‘‘delimiter‘‘:‘）;
end
TIME=c‘;
plot（TIMEM）
clc;
%------------------------------低通濾波器濾除肌電信號------------------------------
Fs=1500;????????????????????????%采樣頻率
fp=80;fs=100;????????????????????%通帶截止頻率，阻帶截止頻率
rp=1.4;rs=1.6;????????????????????%通帶、阻帶衰減
wp=2*pi*fp;ws=2*pi*fs;???
[nwn]=buttord（wpwsrprs‘s‘）;?????%‘s‘是確定巴特沃斯模擬濾波器階次和3dB
[zPk]=buttap（n）;???%設計歸一化巴特沃斯模擬低通濾波器，z為極點，p為零點和k為增益
[bpap]=zp2tf（zPk）??%轉換為Ha（p）bp為分子系數，ap為分母系數
[bsas]=lp2lp（bpapwp）?%Ha（p）轉換為低通Ha（s）并去歸一化，bs為分子系數，as為分母系數

[hsws]=freqs（bsas）;?????????%模擬濾波器的幅頻響應
[bzaz]=bilinear（bsasFs）;?????%對模擬濾波器雙線性變換
[h1w1]=freqz（bzaz）;?????????%數字濾波器的幅頻響應
m=filter（bzazM（:1））;

figure
freqz（bzaz）;title（‘巴特沃斯低通濾波器幅頻曲線‘）;
??????
figure
subplot（211）;
plot（TIMEM（:1））;
xlabel（‘t（s）‘）;ylabel（‘mv‘）;title（‘原始心電信號波形‘）;grid;

subplot（212）;
plot（TIMEm）;
xlabel（‘t（s）‘）;ylabel（‘mv‘）;title（‘低通濾波后的時域圖形‘）;grid;
???
N=512
n=0:N-1;
mf=fft（M（:1）N）;???????????????%進行頻譜變換（傅里葉變換）
mag=abs（mf）;
f=（0:length（mf）-1）*Fs/length（mf）;??%進行頻率變換

figure
subplot（211）
plot（fmag）;axis（[01500150]）;grid;??????%畫出頻譜圖
xlabel（‘頻率（HZ）‘）;ylabel（‘幅值‘）;title（‘心電信號頻譜圖‘）;

mfa=fft（mN）;????????????????????%進行頻譜變換（傅里葉變換）
maga=abs（mfa）;
fa=（0:length（mfa）-1）*Fs/length（mfa）;??%進行頻率變換
subplot（212）
plot（famaga）;axis（[01500150]）;grid;??%畫出頻譜圖
xlabel（‘頻率（HZ）‘）;ylabel（‘幅值‘）;title（‘低通濾波后心電信號頻譜圖‘）;
????
wn=M（:1）;
P=10*log10（abs（fft（wn）.^2）/N）;
f=（0:length（P）-1）/length（P）;
figure
plot（fP）;grid
xlabel（‘歸一化頻率‘）;ylabel（‘功率（dB）‘）;title（‘心電信號的功率譜‘）;
%-----------------帶陷濾波器抑制工頻干擾-------------------??
%50Hz陷波器：由一個低通濾波器加上一個高通濾波器組成??
%而高通濾波器由一個全通濾波器減去一個低通濾波器構成??
Me=100;???????????????%濾波器階數??
L=100;????????????????%窗口長度??
beta=100;?????????????%衰減系數??
Fs=1500;??
wc1=49/Fs*pi;?????%wc1為高通濾波器截止頻率，對應51Hz??
wc2=51/Fs*pi?????;%wc2為低通濾波器截止頻率，對應49Hz??
h=ideal_lp（0.132*piMe）-ideal_lp（wc1Me）+ideal_lp（wc2Me）;?%h為陷波器沖擊響應??
w=kaiser（Lbeta）;??
y=h.*rot90（w）;?????????%y為50Hz陷波器沖擊響應序列??
m2=filter（y1m）;??
??
figure??
subplot（211）;plot（abs（h））;axis（[0?100?0?0.2]）;??
xlabel（‘頻率（Hz）‘）;ylabel（‘幅度（mv）‘）;title（‘陷波器幅度譜‘）;grid;??
N=512;??
P=10*log10（abs（fft（y）.^2）/N）;??
f=（0:length（P）-1）;??
subplot（212）;plot（fP）;??
xlabel（‘頻率（Hz）‘）;ylabel（‘功率（dB）‘）;title（‘陷波器功率譜‘）;grid;??
?????
figure??
subplot?（211）;?plot（TIMEm）;??
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘原始信號‘）;grid;??
subplot（212）;plot（TIMEm2）;??
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘帶阻濾波后信號‘）;grid;??
????
figure??
N=512??
subplot（211）;plot（abs（fft（m））*2/N）;axis（[0?100?0?1]）;??
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘原始信號頻譜‘）;grid;??
subplot（212）;plot（abs（fft（m2））*2/N）;axis（[0?100?0?1]）;??
xlabel（‘t（s）‘）;ylabel（‘幅值‘）;title（‘帶阻濾波后信號頻譜‘）;grid;????
%------------------IIR零相移數字濾波器糾正基線漂移-------------------??
Wp=1.4*2/Fs;?????%通帶截止頻率???
Ws=0.

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件??????118864??2017-12-26?19:55??109a.txt
?????文件????????4359??2018-01-03?16:19??Untitled3.m