数字信号处理课程设计报告

目录
一、课程设计要求
二、设计过程
（1）设计题目
（2）设计源代码
（3）设计结果
（4）结果分析
三、设计总结与心得体会
四、课程设计指导书
一、课程设计要求
1、课程设计指导书
①《数字信号处理（第二版）》，丁玉美等，西安电子科技大学出版社；
②《MATLAB 及在电子信息课程中的应用》，陈怀琛等，电子工业出版社。

2、课程设计内容：
⑴语音信号去噪处理
主要要求：
1)在windows系统下的录音机录制一段1s左右的语音信号作为原声信号，在
MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数；
2)画出语音信号的时域波形，对采样后的语音进行fft变换，得到信号的频谱
特性；对语音信号分别加入正弦噪声和白噪声，画出加噪信号的时域波形和频谱图；
3)根据对加噪语音信号谱分析结果，确定滤除噪声滤波器的技术指标，设计合
适的数字滤波器，并画出滤波器的频域响应；
4)用所设计的滤波器对加噪的信号进行滤波，在同一个窗口画出滤波前后信号
的时域图和频谱图，对滤波前后的信号进行对比，分析信号变化；
5)利用sound(x)回放语音信号，验证设计效果。

⑵语音信号的延时和混响
主要要求：
1)利用Windows下的录音机或其他软件，录制一段自己的语音信号，时间控制
在1s左右，并对录制的信号进行采样；
2)语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；
3)将信号加入延时和混响，再分析其频谱，并与原始信号频谱进行比较；
4)设计几种特殊类型的滤波器：单回声滤波器，多重回声滤波器，全通结构的
混响器，并画出滤波器的频域响应；
5)用自己设计的滤波器对采集的语音信号进行滤波；
6)分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的
信号进行对比，分析信号的变化；
7)回放语音信号。

⑶数字滤波器的设计及实现
主要要求：
1)调用信号产生函数mstg产生三路抑制载波调幅信号相加构成的复合信号st，
观察st的时域波形和幅频特性曲线；
2)由要求将st中的三路调幅信号分离，通过观察st的幅频特性曲线，分别确定
可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率，要求滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB；
3)编程调用MATLAB滤波器设计函数分别设计这三个数字滤波器，并绘图显
示其幅频特性曲线；
4)调用滤波函数filter，用所设计的三个滤波器分别对复合信号st进行滤波，
分离出st中的三路不同载波频率的调幅信号，并绘图显示滤波后信号的时域波形和频谱，观察分离效果。

⑷心电信号的处理
主要要求：
1)在MATLAB软件平台下，给原始的心电信号叠加上噪声或干扰，干扰类型
分为如下几种：白噪声、工频干扰（50Hz）、谐波干扰（二次、三次谐波为主，分别为100Hz、150Hz）绘出叠加噪声后的心电信号时域和频谱图，在视觉上与原始心电信号图形对比，分析频域基本特征变化。

2)给定滤波器的规一化性能指标（参考指标，实际中依据每个同学所叠加噪声
情况而定），例如：通带截止频率wp=0.25*pi, 阻通带截止频率ws=0.3*pi; 通带最大衰减Rp=1 dB; 阻带最小衰减Rs=15 dB；
3)采用窗函数法设计各型FIR滤波器（低通、高通、带通、带阻中的至少2种
类型），来对叠加干扰前后的心电信号进行滤波处理，绘出滤波器的频域响应及滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析或解释；
4)采用双线性变换法利用不同的原型低通滤波器（Butterworth型与切比雪夫I
型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少2种类型）绘出滤波器的频域响应，并用这些数字滤波器对含噪心电信号分别进行滤波处理；比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析或解释；
5)心电信号波形观察、频谱观察，对滤波后的心电信号观察其时域、频域特征
变化。

3、具体要求
⑴、使用 MATLAB（或其它开发工具）编程实现上述内容，写出课程设计报告。

⑵、课程设计报告的内容包括：
①课程设计题目和题目设计要求；
②设计思想和系统功能结构及功能说明；
③设计中关键部分的详细描述和介绍，采用流程图描述关键模块的设计思路；
④总结，包括设计过程中遇到的问题和解决方法，心得体会等；
⑤参考文献；
⑥程序源代码清单。

4、考核方式
课程考核分三部分，一部分是上机率，占 20%；第二部分是检查成绩，最后两次上机为
检查时间，占 50%；第三部分为课程设计报告，占 30%。

注意：
⑴、使用 GUI 界面或混合编程实现仿真程序，酌情加分；
⑵、若发现程序或课程设计报告雷同，一律不及格。

⑶、主要参考资料
[1] S. K. Mitra. Digital Signal Processing: A Computer Based Approach, 3rd Edition [M], New York,
USA: McGraw-Hill, 2000
[2] R. G . Lyons. Understanding Digital Signal Processing, 2nd Edition [M]. New Jersey, USA:
Prentice Hall, 2005
[3] 程佩青. 数字信号处理教程, 第二版[M]. 北京: 清华大学出版社, 2001
[4] 赵树杰等. 数字信号处理[M]. 西安: 西安电子科技大学出版社, 1997
[5] 丁玉美等. 数字信号处理—时域离散随机信号处理[M]. 西安: 西安电子科技大学出版社,
2002
[6] 陈怀琛等. MATLAB 及在电子信息课程中的应用[M], 北京: 电子工业出版社出版, 2002
⑷、课程设计进度安排
序号阶段内容合计（天）
一设计准备 1
二方案选择及初步设计 2
三目标项目设计实现及调试 3
四撰写课程设计报告 2
五上机检查成绩 2
总计（2 周） 10
二.设计过程
第一题、语音信号去噪处理
1.设计要求：
（1）在windows系统下的录音机录制一段1s左右的语音信号作为原声信号，在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数；
（2）画出语音信号的时域波形，对采样后的语音进行fft变换，得到信号的频谱特性；对语音信号分别加入正弦噪声和白噪声，画出加噪信号的时域波形和频谱图；
（3）根据对加噪语音信号谱分析结果，确定滤除噪声滤波器的技术指标，设计合适的数字滤波器，并画出滤波器的频域响应；
（4）用所设计的滤波器对加噪的信号进行滤波，在同一个窗口画出滤波前后信号的时域图和频谱图，对滤波前后的信号进行对比，分析信号变化；
（5）利用sound(x)回放语音信号，验证设计效果。

2.设计步骤：
（1）找到7s的语音信号，利用函数wavread对语音信号进行信号读取；（2）计算样本时刻和频谱图的频率，并进行N+1点FFT变换；
（3）加噪声为5000Hz的正弦信号正弦噪声，采用awgn函数加信噪比为10的高斯白噪声；
（4）设计滤波器；
（5）绘出相应的时域、频域图；
（6）利用sound函数进行原始信号的语音播放，加噪声音播放，以及滤波之后的语言播放。

3.设计实现：
（1）时域图与频谱图（加正弦）
录入原始信号的时域图：
加入正弦信号后的时域图：
滤波后的时域图：
录入原始信号的频域图：
加入正弦信号后的频率图：
滤波后的频域图：
采用巴斯低通滤波器滤除正弦波：
（2）具体代码实现：
[x,fs,bits]=wavread('E:\mcpass.wav');%原信号
n=size(x,1); %提取采样信号的长度
t=(0:length(x)-1)/fs; %计算样本时刻
f=fs*(0:(n+1)/2-1)/n+1; %计算频域图的频率
X=fft(x,n+1); %进行N+1点FFT变换
ts=0:1/fs:(size(x)-1)/fs; %将所加噪声信号的点数调整到与原始信号相同
s=x+0.05*sin(2*pi*5000*ts)'; %加噪声为5000Hz的正弦信号正弦噪声
S=fft(s,n+1); %加正弦噪声后的频域
%正弦滤波
wp=2000/fs*2*pi; %2000为通带截止频率
ws=3000/fs*2*pi; %3000为阻带下限截止频率
Rp=4; %通带波纹
Rs=25; %阻带波纹
T=1/fs;Fs=1/T; %定义采样间隔
Wp=2/T*tan(wp/2); %计算对应的数字频率
Ws=2/T*tan(ws/2);
[N,wn]=buttord(Wp,Ws,Rp,Rs,'s'); %计算滤波器介数和截止频率
[c,d]=butter(N,wn,'s'); %计算滤波器系统函数分子分母系数
[B,A]=bilinear(c,d,Fs); %双线性变换得到数字滤波器系统函数分子分母系数[Hb,Wc]=freqz(B,A);
sf=filter(B,A,s); %对加噪信号进行滤波
Sf=fft(sf,n+1); %对滤波后进行N+1点FFT变换
%绘图部分
figure(3);
plot(fs*Wc/(2*pi),20*log10(abs(Hb)));title('巴斯低通滤波器频域响应图');
xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(3,1,1);
plot(t,x);title('原信号时域')
xlabel('时间(s)');
ylabel('幅度');
figure(2);
subplot(3,1,1);
plot(f,abs(X(1:(n+1)/2)));title('原信号频域')
xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(3,1,2);
plot(t,s);title('加正弦信号后的时域')
xlabel('时间(s)');
ylabel('幅度');
figure(2);
subplot(3,1,2);
plot(f,abs(S(1:(n+1)/2)));title('加正弦信号后的频域图') xlabel('频率(Hz)');
ylabel('幅度');
figure(1)
subplot(3,1,3);
plot(t,sf);title('滤波后的时域图');
xlabel('时间(s)');
ylabel('幅度');
figure(2)
subplot(3,1,3);
plot(f,abs(Sf(1:(n+1)/2))); title('滤波后的频域图'); xlabel('频率(Hz)');
ylabel('幅度');
sound(x);
sound(s);
sound(sf);
（3）时域图与频域图（加白噪声）
加白噪声后的时域图和滤除之后的时域图：
加白噪声和滤除之后的频域图：
采用blackman函数滤波：
具体代码实现：
[x,fs,bits]=wavread('E:\hbsong.wav');
N=size(x,1);
t=(0:length(x)-1)/fs;
f=fs*(0:(N+1)/2-1)/N+1;
X=fft(x,N+1);
%加高斯白噪声
z=awgn(x,20); %对信号加信噪比为10的高斯白噪声N1=size(z,1); %提取采样信号的长度
t=(0:length(z)-1)/fs; %计算样本时刻
f=fs*(0:(N1+1)/2-1)/N1+1;
Z=fft(z,N1+1);
Wp=2500/fs*2*pi;
Ws=3000/fs*2*pi; %计算对应的数字频率
B=Ws-Wp;
n=ceil(1*pi/B);
wc=(Wp+Ws)/2;
b=fir1(n-1,wc/pi,'stop',blackman(n));%blackman窗函数滤波
[H,w]=freqz(b,1);
y=fftfilt(b,z);
t1=(0:length(y)-1)/fs;
Y=fft(y,N1+1);
subplot(2,2,1);
plot(t,z);title('加高斯白噪声后时域图');
subplot(2,2,2);
plot(f,abs(Z(1:(N1+1)/2)));title('滤波前信号频谱图')
figure
plot(fs*w/(2*pi),20*log10);title('blackman函数频域响应图');
title('频率响应')
xlabel('频率（Hz）');
ylabel('幅度');
subplot(2,2,3);
plot(t1,y)
title('滤波后信号时域图');
xlabel('时间(s)');
ylabel('幅度');
subplot(2,2,4);
plot(f,abs(Y(1:(N1+1)/2)));
title('滤波后信号频谱图')
xlabel('频率（Hz）');
ylabel('幅度');
sound(y,fs)
第二题、语音信号的延时和混响
1.设计要求：
（1）利用Windows下的录音机或其他软件，录制一段自己的语音信号，时间控制在1s左右，并对录制的信号进行采样；
（2）语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；
将信号加入延时和混响，再分析其频谱，并与原始信号频谱进行比较；
（3）设计几种特殊类型的滤波器：单回声滤波器，多重回声滤波器，全通结构的混响器，并画出滤波器的频域响应；
（4）用自己设计的滤波器对采集的语音信号进行滤波；
（5）分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；
（6）回放语音信号。

2. 设计步骤：
（1）录入原始声音信号；
（2）计算样本时刻和频谱图的频率，并进行N+1点FFT变换；
（3）加入单回声；
（4）设计单回声滤波器幅频响应函数；
（5）绘出相应时域与频域图；
（6）利用sound函数进行原始信号的语音播放，加单回声后语言播放，以及滤
除之后的语言播放；
3.设计实现：
（1）时域图和频域图（加单回声）原始信号时域图：
滤波后的时域图：
原始信号频域图：
滤波后频域图：
单回声滤波器幅频响应：
（2）具体代码实现（单回声）：
[x,fs,bits]=wavread('E:\mcpass.wav');%原声音信号
n=size(x,1);
t=(0:length(x)-1)/fs;
f=fs*(0:(n+1)/2-1)/n+1;
X=fft(x,n+1);
a=0.6;%单回声滤波
R=fs*a;
B=[1,zeros(1,R-2),a];
A=[1,zeros(1,R-1)];
y = filter(B,A,x);
Y=fft(y,n+1);
[H,W]=freqz(B,A); %求单回声滤波器幅频响应函数%绘图部分
figure(2);
plot(fs*W/(2*pi),20*log10(abs(H)));title('频率响应');
xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,1);
plot(t,x);title('原信号时域');
xlabel('时间(s)');
ylabel('幅度');
subplot(4,1,2);
plot(f,abs(X(1:(n+1)/2)));title('原信号频域');
xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,3);
plot(t,y);title('滤波后时域图');
xlabel('时间(s)');
ylabel('幅度');
figure(1);
subplot(4,1,4);
plot(f,abs(Y(1:(n+1)/2)));title('滤波后频域图');
xlabel('频率(Hz)');
ylabel('幅度');
sound(x,fs,bits);
sound(y,fs,bits);
（3）时域图与频域图（多重回声）
原始信号时域图与多重信号时域图：
原始信号时域图与多重信号频域图：
多重回声滤波器幅频响应：
（4）具体代码实现（多回声）
%原声音信号
[x,fs,bits]=wavread('E:\mcpass.wav');
N=size(x,1);
t=(0:length(x)-1)/fs;
f=fs*(0:(N+1)/2-1)/N+1;
%多重回声
a=0.5;
R=fs*a;
B=[1,zeros(1,R-2)];
A=[1,zeros(1,R-1),a];
yd = filter(B,A,x); %滤波器函数
%频率响应
[Hb,Wc]=freqz(B,A); %求多重回声滤波器幅频响应函数
%绘图部分
figure(2);
plot(fs*Wc/(2*pi),20*log10(abs(Hb)));title('多重滤波器频域响应图') xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,1);
plot(t,x);title('原信号时域');
xlabel('时间(s)');
ylabel('幅度');
X=fft(x,N+1);
subplot(4,1,3);
plot(f,abs(X(1:(N+1)/2)));title('原信号频域');
xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,2);
plot(t,yd);title('多重信号时域图')
xlabel('时间(s)');
ylabel('幅度');
Yd=fft(yd,N+1); %多重回声滤波后信号FFT变换subplot(4,1,4);
plot(f,abs(Yd(1:(N+1)/2)));title('多重信号频域图');
xlabel('频率(Hz)');
ylabel('幅度');
sound(x);
sound(yd);
（5）时域图与频域图（全通滤波器）
原始信号时域图与全通滤波器时域图：
原始信号时域图与全通滤波器频域图：
全通滤波器的频域响应图：
（6）具体代码实现（全通滤波器）：
%原声音信号
[x,fs,bits]=wavread('E:\mcpass.wav');
N=size(x,1);
t=(0:length(x)-1)/fs;
f=fs*(0:(N+1)/2-1)/N+1;
%全通结构混响
a=0.5;
R=fs*a;
B=[a,zeros(1,R-2),1];
A=[1,zeros(1,R-2),a];
yd = filter(B,A,x);
%频率响应
[Hb,Wc]=freqz(B,A); %求全通滤波器幅频响应函数figure(2);
plot(fs*Wc/(2*pi),20*log10(abs(Hb)));title('全通结构频域响应图') xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,1);
plot(t,x);title('原信号时域');
xlabel('时间(s)');
ylabel('幅度');
X=fft(x,N+1);
subplot(4,1,3);
plot(f,abs(X(1:(N+1)/2)));title('原信号频域');
xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,2);
plot(t,yd);title('全通滤波器时域图');
xlabel('时间(s)');
ylabel('幅度');
Yd=fft(yd,N+1);
subplot(4,1,4);
plot(f,abs(Yd(1:(N+1)/2)));title('全通滤波器频域图');
xlabel('频率(Hz)');
ylabel('幅度');
sound(x);
sound(yd);
（7）时域图与频域图（延时与混响）
原始信号时域图：
延时时域图：
混响时域图：
原始信号频域：
延时信号频域图：
混响信号频域图：
（8）具体代码实现（延时与混响）：
[x,fs,bits]=wavread('E:\mcpass.wav');
N=size(x,1);
t=(0:length(x)-1)/fs;
f=fs*(0:(N+1)/2-1)/N+1;
xd=[zeros(300,1);x]; %信号延时
td=(0:length(xd)-1)/fs;
Xd=fft(xd);
Xd1=fftshift(Xd);
dFs =fs/length(xd);
xd1=[x;zeros(300,1)];%信号混响
x1=xd1+xd;
Xhun=fft(x1);%混响信号fft变换
Xh2=fftshift(x1);%平移，中心为0频率
dFs = fs/length(x1);
%绘图部分
figure(1);
subplot(4,1,1);
plot(t,x);title('原信号时域');
xlabel('时间(s)');
ylabel('幅度');
X=fft(x,N+1);
subplot(4,1,2);
plot(f,abs(X(1:(N+1)/2)));title('原信号频域');
xlabel('频率(Hz)');
ylabel('幅度');
subplot(4,1,3);
plot(td,xd);title('延时信号时域图')
xlabel('时间(s)');
ylabel('幅度');
subplot(4,1,4);
plot([-fs/2:dFs: fs/2-dFs],abs(Xd)); title('延时信号频域图') xlabel('频率(Hz)');
ylabel('幅度');
figure(2);
subplot(2,1,1);
plot(td,x1);title('混响信号时域图');
xlabel('时间(s)');
ylabel('幅度');
figure(2);
subplot(2,1,2);
plot([-fs/2:dFs: fs/2-dFs],abs(Xh2)); title('混响信号频域图'); xlabel('频率(Hz)');
ylabel('幅度');
sound(x);
sound(xd);
sound(x1);
第三题、数字滤波器的设计及实现
1.设计要求：
（1）调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，观察st的时域波形和幅频特性曲线；
（2）要求将st中的三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率，要求滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB；
（3）编程调用MATLAB滤波器设计函数分别设计这三个数字滤波器，并绘图显示其幅频特性曲线；
（4）调用滤波函数filter，用所设计的三个滤波器分别对复合信号st进行滤波，分离出st中的三路不同载波频率的调幅信号，并绘图显示滤波后信号的时域波形和频谱，观察分离效果。

2.设计步骤：
（1）产生三路调幅信号；
（2）三路信号相加组成复合信号；
（3）分别设计低通、带通、高通滤波器，对三路信号滤波；
（4）绘出相应的时域、频域图；
（5）观察分离效果；
3.设计实现：
（1）时域图与频域图（第一路）：
第二路：
第三路：
4.具体代码实现：
N=1600; %N为信号st的长度。

Fs=10000;T=1/Fs;Tp=N*T; %采样频率Fs=10kHz，Tp为采样时间
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10; %第1路调幅信号的载波频率fc1=1000Hz, fm1=fc1/10; %第1路调幅信号的调制信号频率fm1=100Hz fc2=Fs/20; %第2路调幅信号的载波频率fc2=500Hz fm2=fc2/10; %第2路调幅信号的调制信号频率fm2=50Hz fc3=Fs/40; %第3路调幅信号的载波频率fc3=250Hz, fm3=fc3/10; %第3路调幅信号的调制信号频率fm3=25Hz xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号
st=xt1+xt2+xt3; %三路调幅信号相加
fxt=fft(st,N); %计算信号st的频谱
figure(1);
subplot(2,1,1);
plot(t,st);title('三路调幅信号时域图')
grid;
xlabel('t/s');ylabel('s(t)');
axis([0,Tp/2,min(st),max(st)]);
title('(a) s(t)的波形');
subplot(2,1,2);
stem(f,abs(fxt)/max(abs(fxt)),'.');
grid;
title('(b) s(t)的频谱');
axis([0,Fs/5,0,1.2]);
xlabel('f/Hz');ylabel('幅度');
Fs=10000;T=1/Fs; %采样频率
%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st
%低通滤波器设计与实现
fp=280;fs=450;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
[N,wp]=ellipord(wp,ws,rp,rs);
[B,A]=ellip(N,rp,rs,wp);
y1t=filter(B,A,st);
[Hb,Wc]=freqz(B,A);
figure(2);
subplot(2,1,1);
plot(Fs*Wc/(2*pi),20*log10(abs(Hb)));title('低通滤波器频域响应图')
subplot(2,1,2);
plot(t,y1t);title('第一路低频信号时域图')
%带通滤波器设计与实现
fpl=440;fpu=560;fsl=275;fsu=900;
wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;
[N,wp]=ellipord(wp,ws,rp,rs);
[B,A]=ellip(N,rp,rs,wp);
y2t=filter(B,A,st);
[Hb1,Wc1]=freqz(B,A);
figure(3);
subplot(2,1,1);
plot(Fs*Wc1/(2*pi),20*log10(abs(Hb1)));title('带通滤波器频域响应图')
subplot(2,1,2);plot(t,y2t);title('第二路中频信号时域图')
%高通滤波器设计与实现
fp=890;fs=600;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp [B,A]=ellip(N,rp,rs,wp,'high'); %调用ellip计算椭圆带通DF系统函数系数向量B和A
y3t=filter(B,A,st); %滤波器软件实现
[Hb2,Wc2]=freqz(B,A);
figure(4);
subplot(2,1,1);
plot(Fs*Wc2/(2*pi),20*log10(abs(Hb2)));title('高通滤波器频域响应图')
subplot(2,1,2);plot(t,y3t);title('第三路高频信号时域图')
第四题、心电信号的处理
1.设计要求：
（1）在MATLAB软件平台下，给原始的心电信号叠加上噪声或干扰，干扰类型分为如下几种：白噪声、工频干扰（50Hz）、谐波干扰（二次、三次谐波为主，分别为100Hz、150Hz）绘出叠加噪声后的心电信号时域和频谱图，在视觉上与原始心电信号图形对比，分析频域基本特征变化。

（2）给定滤波器的规一化性能指标（参考指标，实际中依据每个同学所叠加噪声情况而定），例如：通带截止频率wp=0.25*pi, 阻通带截止频率ws=0.3*pi; 通带最大衰减Rp=1 dB; 阻带最小衰减Rs=15 dB；
（3）采用窗函数法设计各型FIR滤波器（低通、高通、带通、带阻中的至少2种类型），来对叠加干扰前后的心电信号进行滤波处理，绘出滤波器的频域响应及滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析或解释；（4）采用双线性变换法利用不同的原型低通滤波器（Butterworth型与切比雪夫I型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少2种类型）绘出滤波器的频域响应，并用这些数字滤波器对含噪心电信号分别进行滤波处理；比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析或解释；
（5）心电信号波形观察、频谱观察，对滤波后的心电信号观察其时域、频域特征变化。

2.设计步骤：
（1）读取原始心电信号；
（2）心电信号加50HZ工频信号；
（3）心电信号加入谐波；
（4）心电信号加入白噪声；
（5）设计相应的滤波器进行不同噪声的滤除；
（6）绘出相应的时域频域图；
3.设计实现：
（1）各时域频域图：
原始心电信号时域图以及加入个噪声后的时域图：
原始心电信号频域图以及加入个噪声后的时频域图：
（2）具体代码实现：
%原心电信号
xl=load('E:\心电信号-新.txt');
x=xl(:,2);
N=size(x,1);
fs=1000;
t=(0:length(x)-1)/fs;
f=fs/N*(0:(N+1)/2-1)+1;
X=fft(x,N);
%加50hz工频干扰
y=x+sin(2*50*pi*t)';
Y=fft(y,N);
%加100hz、150hz谐波干扰
y1=x+sin(2*100*pi*t)'+cos(2*150*pi*t)';
Y1=fft(y1,N);
%加白噪声
yb=awgn(x,20);
Yd=fft(yb,N);
figure(1);
subplot(4,1,1);
plot(t,x);title('原信号时域图')
xlabel('时间(s)');
ylabel('幅度');
figure(2);
subplot(4,1,1)
plot(f,abs(X(1:(N+1)/2)));title('原信号频域图')
xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,2);
plot(t,y);title('加50hz工频信号时域图')
xlabel('时间(s)');
ylabel('幅度');
figure(2);
subplot(4,1,2);
plot(f,abs(Y(1:(N+1)/2)));title('加50hz工频信号频域图') xlabel('频率(Hz)');
ylabel('幅度');
figure(1)
subplot(4,1,3);
plot(t,y1);title('加谐波后的时域图')
xlabel('时间(s)');
ylabel('幅度');
figure(2)
subplot(4,1,3);
plot(f,abs(Y1(1:(N+1)/2)));title('加谐波后的频域图') xlabel('频率(Hz)');
ylabel('幅度');
figure(1);
subplot(4,1,4);
plot(t,y1);title('加白噪后的时域图');
figure(2);
subplot(4,1,4);
plot(f,abs(Y1(1:(N+1)/2)));title('加白噪后的频域图');
（3）滤除部分：
滤除50HZ信号后的时域图（blackman低通滤波器）：
滤除50HZ信号后的频域图：
blackman低通滤波器滤波：
具体代码：
%原心电信号
xl=load('E:\心电信号-新.txt');
x=xl(:,2);
N=size(x,1);
fs=1000;
t=(0:length(x)-1)/fs;
f=fs/N*(0:(N+1)/2-1)+1;%f=fs/N*(0:N-1);
X=fft(x,N);
%加50hz工频干扰
y=x+sin(2*50*pi*t)';
Y=fft(y,N);
%blackman低通滤波器
wp=30/fs*pi;ws=40/fs*pi;
B=ws-wp;
M=ceil(12*pi/B)-1;
bl=fir1(M,(ws+wp)/2/pi,'low',blackman(M+1));
yl=fftfilt(bl,y);
Yl=fft(yl,N);
figure(1);
subplot(3,1,1);
plot(t,x);title('原信号时域图')
figure(2);
subplot(3,1,1)
plot(f,abs(X(1:(N+1)/2)));%plot(f,abs(X));
title('原信号频域图')
figure(1);
subplot(3,1,2);
plot(t,y);title('加50hz工频信号时域图')
figure(2);
subplot(3,1,2);
plot(f,abs(Y(1:(N+1)/2)));%plot(f,abs(Y));
title('加50hz工频信号频域图')
figure(1);
subplot(3,1,3);
plot(t,yl);title('滤波后信号时域图')
figure(2)
subplot(3,1,3);
plot(f,abs(Yl(1:(N+1)/2)));%plot(f,abs(Yl));
title('滤波后信号频域图')
figure(3)
[h2,w2]=freqz(bl,1);
plot(w2/(2*pi)*fs,20*log10(abs(h2)));title('blackman低通滤波器频率响应图')；滤除50HZ信号后的时域图（切比雪夫滤波器）：
滤除50HZ信号后的频域图（切比雪夫滤波器）：
切比雪夫滤波器：
具体代码实现：
%原心电信号
xl=load('E:\心电信号-新.txt');
x=xl(:,2);
N=size(x,1);
fs=1000;
t=(0:length(x)-1)/fs;
f=fs/N*(0:(N+1)/2-1)+1;
X=fft(x,N);
%加50hz工频干扰
y=x+sin(2*50*pi*t)';
Y=fft(y,N);
%切比雪夫滤波器
wp=20/fs*2*pi;ws=30/fs*2*pi;
wp1=2*fs*tan(wp/2);
ws1=2*fs*tan(ws/2);
Rp=1;
Rs=15;
[N1,w]=cheb1ord(wp1,ws1,Rp,Rs,'s'); [B1,A1]=cheby1(N1,1,w,'low','s'); [B,A]=bilinear(B1,A1,fs);
[Hk,wk]=freqz(B,A);
yf=filter(B,A,y);
Yf=fft(yf,N);
figure(1);
subplot(3,1,1);
plot(t,x);title('原信号时域图')
figure(2);
subplot(3,1,1)
plot(f,abs(X(1:(N+1)/2)));title('原信号频域图')
figure(1);
subplot(3,1,2);
plot(t,y);title('加50hz工频信号时域图')
figure(2);
subplot(3,1,2);
plot(f,abs(Y(1:(N+1)/2)));title('加50hz工频信号频域图')
figure(1);
subplot(3,1,3)
plot(t,yf);title('滤波后信号时域图')
figure(2);
subplot(3,1,3)
plot(f,abs(Yf(1:(N+1)/2)));title('滤波后信号频域图')
figure(3)
plot(fs*wk/(2*pi),20*log10(abs(Hk)));title('切比雪夫滤波器频域响应图')；滤除白噪声后的时域图（blackman低通滤波器）：
滤除白噪声后的频域图（blackman低通滤波器）：
blackman低通滤波器：
具体代码实现：
%原心电信号
xl=load('E:\心电信号-新.txt');
x=xl(:,2);
N=size(x,1);
fs=1000;
t=(0:length(x)-1)/fs;
f=fs/N*(0:(N+1)/2-1)+1;
X=fft(x,N);
%加白噪声
yb=awgn(x,20);
Yd=fft(yb,N);
%blackman低通滤波器
wp=70/fs*pi;ws=90/fs*pi;
Bl=ws-wp;
Ml=ceil(12*pi/Bl)-1;
bl=fir1(Ml,(ws+wp)/2/pi,'low',blackman(Ml+1));
yl=fftfilt(bl,y1);
Yl=fft(yl,N);
%绘图部分
figure(1);
subplot(3,1,1);
plot(t,x);title('原信号时域图')
figure(2);
subplot(3,1,1)
plot(f,abs(X(1:(N+1)/2)));title('原信号频域图')
figure(1);
subplot(3,1,2);
plot(t,y1);title('加白噪后的时域图');
figure(2);
subplot(3,1,2);
plot(f,abs(Y1(1:(N+1)/2)));title('加白噪后的频域图');
figure(1);
subplot(3,1,3);
plot(t,yl);title('滤波后信号时域图')
figure(2)
subplot(3,1,3);
plot(f,abs(Yl(1:(N+1)/2)));%plot(f,abs(Yl));
title('滤波后信号频域图')
figure(3)
[h2,w2]=freqz(bl,1);
plot(w2/(2*pi)*fs,20*log10(abs(h2)));title('blackman低通滤波器频率响应图');
三、设计总结与心得体会
在课程设计的这段时间，我获益匪浅。

不但进一步掌握了数字信号处理的基础知识及MATLAB的基本操作，还使我了解了信号的产生、采样及频谱分析的方法。

我进一步了解到凡事都需要耐心，细心仔细是成功的保证。

虽然在做的过程中遇到了一些问题，但是我都通过自己的努力解决了它们，证明了自己的能力。

这次课程设计对我各方面的综合能力有了很大的提高，对我以后的工作，实践都有很大的帮助。

在此次课程设计当中，我经常把C语言的语法知识照搬到MA TALAB设计中，从而导致调试失败，所以下次用此类语言做课程设计时，应事先学习下这类语言的基本语法，以免与其他语言相混淆。

还有就是有些不定参数存在时，可先取定值，用于调试，这样可以节约调试时间，从而提高效率。

本次课程设计不但让我又学到了一些知识，而且也提高了我的综合能力。

使我在各个方面都得到了锻炼，以后有这样的机会一定会更加的很好利用，它不仅可以提高学习的针对性而且可以很好的锻炼动手能力以及自己的逻辑设计能力和处理问题的能力，希望在以后这方面的能力会很好的加强。

四、课程设计指导书
[1] 《数字信号处理(第二版)》.丁玉美等西安电子科技大学出版社
[2] 《数字信号处理及其MATLAB实现》，陈怀琛等译，电子工业出版社；
[3] 《MATLAB及在电子信息课程中的应用》，陈怀琛等，电子工业出版社忽略此处..。