《DSP课程设计》word版.docx

DSP课程设计作者：9月9日-10月9日中南大学通信班目录1、绪论11.1简述21.2设计任务32、GUI图形界面设计42.1设计原理52.2设计效果63、语音信号去噪处理43.1设计流程图53.2设计步骤与过程64、语音信号的延时与混响44.1设计流程图54.2设计步骤与过程65、数字滤波器的设计与实现45.1设计流程图55.2设计步骤与过程66、心电信号的处理46.1设计流程图56.2设计步骤与过程67、设计小结47.1结果分析57.2心得体会98、附录48.1仿真程序58.2参考文献1、绪论 1.1简述数字信号处理课程是通信专业的一门重要专业基础课，是信息的数字化处理、存储和应用的基础。通过该课程的课程设计实践，使学生对信号与信息的采集、处理、传输、显示、存储、分析和应用等有一个系统的掌握和理解；巩固和运用在数字信号处理课程中所学的理论知识和实验技能，掌握数字信号处理的基础理论和处理方法，提高分析和解决信号与信息处理相关问题的能力，为以后的工作和学习打下基础。 本次设计是用MATLAB语言对语音信号进行采样分析，并设计数字滤波器对信号进行滤波，比较滤波前后信号特性的变化。用MATLAB提供的GUI界面开发环境设计用户图形界面使布局编程简化。 1.2设计任务1、 语音信号去噪处理主要要求：1) 在windows系统下的录音机录制一段1s左右的语音信号作为原声信号，在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数；2) 画出语音信号的时域波形，对采样后的语音进行fft变换，得到信号的频谱特性；对语音信号分别加入正弦噪声和白噪声，画出加噪信号的时域波形和频谱图；3) 根据对加噪语音信号谱分析结果，确定滤除噪声滤波器的技术指标，设计合适的数字滤波器，并画出滤波器的频域响应；4) 用所设计的滤波器对加噪的信号进行滤波，在同一个窗口画出滤波前后信号的时域图和频谱图，对滤波前后的信号进行对比，分析信号变化；5) 利用sound(x)回放语音信号，验证设计效果。2、 语音信号的延时和混响主要要求：1) 利用Windows下的录音机或其他软件，录制一段自己的语音信号，时间控制在1s左右，并对录制的信号进行采样；2) 语音信号的频谱分析，画出采样后语音信号的时域波形和频谱图；3) 将信号加入延时和混响，再分析其频谱，并与原始信号频谱进行比较；4) 设计几种特殊类型的滤波器：单回声滤波器，多重回声滤波器，全通结构的混响器，并画出滤波器的频域响应；5) 用自己设计的滤波器对采集的语音信号进行滤波；6) 分析得到信号的频谱，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；7) 回放语音信号。3、 数字滤波器的设计及实现主要要求：1) 调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，观察st的时域波形和幅频特性曲线；2) 要求将st中的三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率，要求滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB； 3) 编程调用MATLAB滤波器设计函数分别设计这三个数字滤波器，并绘图显示其幅频特性曲线； 4) 调用滤波函数filter，用所设计的三个滤波器分别对复合信号st进行滤波，分离出st中的三路不同载波频率的调幅信号，并绘图显示滤波后信号的时域波形和频谱，观察分离效果。4、 心电信号的处理主要要求：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、图形界面的设计2.1设计原理 图形用户界面GUI (Graphics User Interface) 是由各种图形对象,如图形窗口、图轴、菜单、按钮、文本框等构建的用户界面,是人机交流信息的工具和方法。GUI 设计即可以基本的MATLAB 程序设计为主,也可以鼠标为主利用GUIDE 工具进行设计。GUIDE ( Graphics User Interface Design Environ2ment) 是一个专用于GUI 程序设计的快速开发环境,使用者通过鼠标就能迅速地产生各种GUI 控件,并随心所欲地改变它们的外形、大小及颜色等,从而帮助用户方便地设计出各种符合要求的图形用户界面。调用GUIDE 的方法有2 种,在MATLAB 命令窗口中输入guide 命令,或在MATLAB 主菜单中点击File New GUI 即可打开一个可编辑的新窗口。在GUIDE 设计环境中,需要用到的工具有属性编辑器、控件布置编辑器、菜单编辑器、对象浏览器、网格标尺设置编辑器以及GUI 应用属性设置编辑器等。GUI 设计面板是GUI 设计工具应用的平台,面板上部提供了菜单和常用工具按钮,左边提供了多种GUI 控件,如按钮、单选按钮、复选框、文本框等。进行 GUI 设计时,首先单击GUI 面板左边所需的控件,然后在右边的图形界面编辑区中再次单击某一恰当的位置,这时将在该位置上为图形界面添加一相应的控件,接下来,通过属性编辑器和对齐编辑器对各控件设置相关属性和进行界面布置,以完善界面功能。2.2设计效果3、语音信号去噪处理3.1设计流程图开始语音信号的采集（wavread函数），画时域图快速傅里叶变换，并且画频谱图设定滤波器性能指标，通带截止频率fb=，阻带截止频率fc=，通带波纹Ap=，阻带波纹As=脉冲响应不变法法设计巴特沃斯滤波器窗函数法设计FIR滤波器验证并进行频谱分析脉冲响应不变法巴特沃斯滤波器对语音信号去噪流程图开始语音信号的采集（wavread函数），画时域图快速傅里叶变换，并且画频谱图设定滤波器性能指标，通带截止频率fb=1100，阻带截止频率fc=1200，通带波纹Ap=1，阻带波纹As=20脉冲响应不变法法设计巴特沃斯滤波器验证并进行频谱分析设计好的滤波器进行滤波处理比较滤波前后语音信号的波形及频谱回放语音信号结束3.2设计步骤与过程(1)语音信号的获取选择一个wav文件作为分析的对象，可以利用Windows下的录音机或其他软件，录制一段自己的话音，在MATLAB中，y,fs,bits=wavread(Blip,N1 N2);用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。N1 N2表示读取的值从N1点到N2点的值。sound(y); 用于对声音的回放。向量y则就代表了一个信号，也即一个复杂的“函数表达式”，也可以说像处理一个信号的表达式一样处理这个声音信号。下面是语音信号在MATLAB中的语言程序，它实现了语音的读入与打开，并绘出了语音信号时域波形，然后对语音信号进行频谱分析。在MATLAB中，可以利用函数fft对信号进行快速傅里叶变化，得到信号的频谱特性。在频谱特性中分析最大值的位置（可能有几个），它代表的频率和时域的采样时间有关，相邻的两点之间的距离为。其中，N是离散傅里叶变换用的点数，是采样的时间，前面在读取 wav文件时得到了采样频率。既然知道了该声波的频谱，按频率就可以反演它的时域值，利用以上分析的主要峰值来重构声波。由于没有考虑相位和其他的频谱分量，所以波形和原来的波形相差甚大，但大体的频率是没有错的。 x,fs,bits=wavread(H:bat.wav); N=length(x) fssound(x,fs,bits); %播放语音信号y1=fft(x,n); %对信号做2048点FFT变换f=fs*(0:2047)/n;figure(1)magy1=abs(y1);angy1=angle(y1);subplot(3,1,1),plot(x);title(原始信号波形)subplot(3,1,2),plot(magy1);title(原始信号幅值)subplot(3,1,3),plot(angy1);title(原始信号相位)figure(2)freqz(x) %绘制原始语音信号的频率响应图title(频率响应图)figure(3)plot(f,abs(y1(1:2048);title(原始语音信号频谱)xlabel(Hz);ylabel(fudu); axis(0 4500 0 400)（2）wav语音信号加噪声在MATLAB软件平台下，给原始的语音信号叠加上噪声，噪声类型分为如下几种：（1）单频噪色（正弦干扰）；（2）高斯随机噪声。绘出加噪声后的语音信号时域和频谱图，在视觉上与原始语音信号图形对比，也可通过Windows播放软件从听觉上进行对比，分析并体会含噪语音信号频谱和时域波形的改变。由于原理相似这里以正弦干扰为例加以介绍。x,fs,bits=wavread( H:bat.wav); %读取语音信号的数据，赋给变量xN=length(x);y1=fft(x,n); %对信号做n点FFT变换f=fs*(0:511)/n;t=(0:length(x)-1)/22050;x1=0.05*sin(2*pi*10000*t);x2=x+x1;sound(x2,fs,bits);figure(1)subplot(2,1,1)plot(x) %做原始语音信号的时域图形title(原语音信号时域图)subplot(2,1,2)plot(x2) %做原始语音信号的时域图形title(加正弦噪声后语音信号时域图)xlabel(time n);ylabel(fudu);y2=fft(x2,n);figure(2)subplot(2,1,1)plot(abs(y1)title(原始语音信号频谱);xlabel(Hz);ylabel(fudu);subplot(2,1,2)plot(abs(y2)title(加噪语音信号频谱);xlabel(Hz);ylabel(fudu);axis(0 4500 0 300);wavwrite(x2,fs, H:bat加噪.wav);（3）巴特沃斯低通滤波对加入高斯随机噪声和正弦噪声的语音信号进行滤波。用双线性变换法设计了巴特沃斯数字低通IIR滤波器对两加噪语音信号进行滤波，并绘制了巴特沃斯低通滤波器的幅度图和两加噪语音信号滤波前后的时域图和频谱图。fb = 1000;fc = 1200;fs = 22050;wp=0.1*pi;ws=0.4*pi;Rp=1;Rs=15;Fs=22050;Ts=1/Fs;wp1=2/Ts*tan(wp/2); %将模拟指标转换成数字指标ws1=2/Ts*tan(ws/2); N,Wn=buttord(wp1,ws1,Rp,Rs,s); %选择滤波器的最小阶数Z,P,K=buttap(N); %创建butterworth模拟滤波器Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,Wn); bd,ad=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换h,w=freqz(bd,ad); figure(1)subplot(111);plot(w*fs/(2*pi),abs(h)grid;title(滤波器的性能分析); pause;figure(2)x,fs,bits=wavread( H:bat加噪.wav );n=length(x);f=fs*(0:(n/2-1)/n;X=fft(x);z=filter(bd,ad,x);subplot(211);plot(x);title(原始信号的波形);subplot(212);plot(z);title(滤波后信号的波形);pause;figure(3)sound(z,fs,bits);subplot(211);plot(f,abs(X(1:n/2);title(原始信号的频谱);xlabel(Hz);Z=fft(z);subplot(212);plot(f,abs(Z(1:n/2); title(滤波后的信号频谱);xlabel(Hz);wavwrite(z,fs, H:bat巴滤.wav);（4）汉明窗的FIR低通滤波使用窗函数法,选用海明窗设计了数字FIR低通滤波器对加了正弦噪声的语音信号进行滤波，并绘制了滤波器滤波后的语音信号时域图和频谱图。%FIR滤波 x,fs,bits=wavread(H:bat加噪.wav);n=length(x);wp=0.25*pi;ws=0.3*pi;wdelta=ws-wp;N=ceil(6.6*pi/wdelta); %取整t=0:(size(x)-1);wn=(0.2+0.3)*pi/2;b=fir1(N,wn/pi,hamming(N+1); %选择窗函数，并归一化截止频率f1=fftfilt(b,x);figure(1)freqz(b,1,512)h1,w1=freqz(b,1);plot(w1*fs/(2*pi),20*log10(abs(h1);figure(2)subplot(2,1,1)plot(t,x)title(滤波前的时域波形);subplot(2,1,2)plot(t,f1);title(滤波后的时域波形);sound(f1); %播放滤波后的语音信号F0=fft(f1,1024);f=fs*(0:511)/1024;figure(3)y2=fft(x,1024);subplot(2,1,1);plot(f,abs(y2(1:512); %画出滤波前的频谱图title(滤波前的频谱)xlabel(Hz);ylabel(fuzhi);subplot(2,1,2)F1=plot(f,abs(F0(1:512); %画出滤波后的频谱图title(滤波后的频谱)xlabel(Hz);ylabel(fuzhi);wavwrite(f1,fs,h:batF滤.wav); 4、语音信号的延时与混响4.1设计流程图4.2设计步骤与过程一、信号采样（1）采样频率采样频率是指计算机每秒钟采集多少个声音样本，采样频率越高，即采样的间隔时间越短，则在单位时间内计算机得到的声音样本数据就越多，对声音波形的表示也越精确。只有采样频率高于声音信号最高频率的两倍时，才能把数字信号表示的声音还原成为原来的声音。（2）采样位数即采样值或取样值，用来衡量声音波动变化的参数，是指声卡在采集和播放声音文件时所使用数字声音信号的二进制位数。声卡的位客观地反映了数字声音信号对输入声音信号描述的准确程度。声卡的主要的作用之一是对声音信息进行录制与回放，在这个过程中采样的位数和采样的频率决定了声音采集的质量。二、混响与延迟(1)混响效果主要是用于增加音源的融合感。自然音源的延时声阵列非常密集、复杂，所以模拟混响效果的程序也复杂多变。常见参数有以下几种：混响时间：能逼真的模拟自然混响的数码混响器上都有一套复杂的程序，其上虽然有很多技术参数可调，然而对这些技术参数的调整都不会比原有的效果更为自然，尤其是混响时间。高频滚降：此项参数用于模拟自然混响当中，空气对高频的吸收效应，以产生较为自然的混响效果。一般高频混降的可调范围为0.11.0。此值较高时，混响效果也较接近自然混响；此值较低时，混响效果则较清澈。扩散度：此项参数可调整混响声阵密度的增长速度，其可调范围为010，其值较高时，混响效果比较丰厚、温暖；其值较低时，混响效果则较空旷、冷僻。预延时：自然混响声阵的建立都会延迟一段时间，预延时即为模拟次效应而设置。声阵密度：此项参数可调整声阵的密度，其值较高时，混响效果较为温暖，但有明显的声染色；其值较低时，混响效果较深邃，切声染色也较弱。频率调制：这是一项技术性的参数，因为电子混响的声阵密度比自然混响稀疏，为了使混响的声音比较平滑、连贯，需要对混响声阵列的延时时间进行调制。此项技术可以有效的消除延时声阵列的段裂声，可以增加混响声的柔和感。调治深度：指上述调频电路的调治深度。(2)延时就是将音源延迟一段时间后，再欲播放的效果处理。依其延迟时间的不同，可分别产生合唱、镶边、回音等效果。当延迟时间在335ms之间时人耳感觉不到滞后音的存在，并且他与原音源叠加后，会因其相位干涉而产生梳状滤波效应，这就是镶边效果。如果延迟时间在50ms以上时，其延迟音就清晰可辨，此时的处理效果才是回音。回音处理一般都是用于产生简单的混响效果。延时、合唱、镶边、回音等效果的可调参数都差不多，具体有以下几项：*延时时间（Dly），即主延时电路的延时时间调整。*反馈增益（FB Gain），即延时反馈的增益控制。*反馈高频比（Hi Ratio），即反馈回路上的高频衰减控制。*调制频率（Freq），指主延时的调频周期。*调制深度（Depth），指上述调频电路的调制深度。*高频增益（HF），指高频均衡控制。*预延时（Ini Dly），指主延时电路预延时时间调整。*均衡频率（EQ F），这里的频率均衡用于音色调整，此为均衡的中点频率选择。由于延时产生的效果都比较复杂多变，如果不是效果处理专家，建议使用设备提供的预置参数，因为这些预置参数给出的处理效果一般都比较好。三、离散傅立叶变换在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅立叶变换和逆变换。下面介绍这些函数。函数FFT用于序列快速傅立叶变换。函数的一种调用格式为 y=fft(x)其中，x是序列，y是序列的FFT，x可以为一向量或矩阵，若x为一向量，y是x的FFT。且和x相同长度。若x为一矩阵，则y是对矩阵的每一列向量进行FFT。如果x长度是2的幂次方，函数fft执行高速基2FFT算法；否则fft执行一种混合基的离散傅立叶变换算法，计算速度较慢。函数FFT的另一种调用格式为 y=fft(x,N)式中，x，y意义同前，N为正整数。函数执行N点的FFT。若x为向量且长度小于N，则函数将x补零至长度N。若向量x的长度大于N，则函数截短x使之长度为N。若x 为矩阵，按相同方法对x进行处理。经函数fft求得的序列y一般是复序列，通常要求其幅值和相位。MATLAB提供求复数的幅值和相位函数：abs，angle，这些函数一般和FFT同时使用。函数abs(x)用于计算复向量x的幅值，函数angle(x)用于计算复向量的相角，介于 和 之间，以弧度表示。函数unwrap(p)用于展开弧度相位角p ,当相位角绝对变化超过 时，函数把它扩展至 。用MATLAB工具箱函数fft进行频谱分析时需注意：（1） 函数fft返回值y的数据结构对称性。（2） 频率计算。（3） 作FFT分析时，幅值大小与FFT选择点数有关，但不影响分析结果。四、滤波器设计本次设计用IIR滤波器对信号进行滤波，函数名为filter函数filter的调用格式为 y=filter(b,a,x)该格式采用数字滤波器对数据进行滤波，既可以用于IIR滤波器，也可以用于FIR滤波器。其中向量b和a分别表示系统函数的分子、分母多项式的系数，若a1，此时表示FIR滤波器，否则就是IIR滤波器。该函数是利用给出的向量b和a，对x中的数据进行滤波，结果放入向量y。（1）单回声滤波器回音可以由简单的延时单元产生。直达声和在R抽样周期后出现的一种单个回音，可以用FIR滤波器产生，微分方程为：yn=xn+xn-R |1传输函数为：H(z)=1+z-R传递函数的幅频响应形状象梳子，这种滤波器又叫梳状滤波器。（2）多重回声滤波器为了产生以间隔R个抽样周期分开的具有指数衰减振幅的多重回声，可用一个以下形式传输函数的FIR滤波器： 无限个振幅以指数衰减间隔为R个抽样周期的多重回声可用以下形式传输函数的IIR滤波器生成： IIR多重回声滤波器的基本频率FR=Fs/R，通常锁定在伴音设备的基频上，比如基鼓拍子。（3）无限个回声滤波器的系统函数: , 1 （4）全通结构的混响器的系统函数:全通滤波器的传递函数公式为H(Z)=y(Z)/X(Z)=(-K+Z(-m)/(1-K*Z(-m)其中m为回声延时取样，k为反馈系数。用直接1型表示这个传递函数则为：y(n)=k*x(n)+x(n-m)+k*y(n-m)可见其实际上是一个简单的IIR滤波器，时间n的输出有时间n的输入和m点之前的输入与输出计算而得。由于这个IIR滤波器的频率响应为水平直线，所以被称为全通滤波器。原始语音信号时域波形和频谱y,fs,bits=wavread();%用于读取语音，采样值放在向量y中，fs表示采样频率(Hz)，bits表示采样位数。sound(y,fs,bits);%播放原始信号Y=fft(y);%对原始信号取傅里叶变换subplot(2,2,1:2);plot(y);title(原始信号时域波形);subplot(2,2,3);plot(abs(Y);title(原始信号幅频);subplot(2,2,4);plot(angle(Y);title(原始信号相频);采样后语音信号的时域波形和频谱y,fs,bits=wavread();Y=fft(y,6001);%抽取6001点subplot(2,2,1:2);plot(y);title(采样信号时域波形);subplot(2,2,3);plot(abs(Y);title(采样信号幅频);subplot(2,2,4);plot(angle(Y);title(采样信号相频);sound(y,fs,bits);%回放采样信号对采样后的信号延时y,fs,bits=wavread();z=zeros(3000,1);y;%延迟3000Z=fft(z,6001);subplot(2,2,1:2);plot(z);title(延时后时域波形);subplot(2,2,3);plot(abs(Z);title(延时后幅频);subplot(2,2,4);plot(angle(Z);title(延时后相频);sound(y,fs,bits);%回放延迟信号对采样后的信号混响y,fs,bits=wavread();z=zeros(3000,1);y;y=y;zeros(3000,1);%补3000y1=y+z;%将原始信号与延迟信号相加，产生混响信号Y1=fft(y1,6001);subplot(2,2,1:2);plot(y1);title(混响的时域波形);subplot(2,2,3);plot(abs(Y1);title(混响的幅频);subplot(2,2,4);plot(angle(Y1);title(混响的相频);sound(y,fs,bits);%回放混响信号单回声滤波器程序设计y,fs,bits=wavread();z=zeros(3000,1);y;y=y;zeros(3000,1);a=0.5;R=5000;%滤波器阶数设置，其值越高，回声越明显Bz1=1,zeros(1,R-1),a;%单回声滤波器系统函数分子Az1=(1);%单回声滤波器系统函数分母y2=filter(Bz1,Az1,y);%单回声滤波器系统函数 Y2=fft(y2,6001);h,w=freqz(Bz1,Az1);%求设计的滤波器频谱subplot(3,2,1);plot(abs(h);title(单回声滤波器幅频响应);subplot(3,2,2); plot(angle(h);title(单回声滤波器相频响应);subplot(3,2,3:4);plot(y2);title(单回声滤波器时域图 );subplot(3,2,5);plot(abs(Y2);title(单回声滤波器幅频 );subplot(3,2,6);plot(angle(Y2);title(单回声滤波器相频 );sound(y2,fs,bits); 多重回声滤波器程序设计y,fs,bits=wavread();z=zeros(3000,1);y;y=y;zeros(3000,1);a=0.5;N=5;R=9000;%滤波器阶数设置，其值越高，回声越明显Bz1=1,zeros(1,N*R-1),-aN;%多重回声滤波器系统函数分子Az1=1,zeros(1,R-1),-a;%多重回声滤波器系统函数分母y2=filter(Bz1,Az1,y);Y2=fft(y2,6001);h,w=freqz(Bz1,Az1);subplot(3,2,1);plot(abs(h);title(多重回声滤波器幅频响应);subplot(3,2,2); plot(angle(h);title(多重回声滤波器相频响应);subplot(3,2,3:4);plot(y2);title(多重回声滤波器时域图 );subplot(3,2,5);plot(abs(Y2);title(多重回声滤波器幅频 );subplot(3,2,6);plot(angle(Y2);title(多重回声滤波器相频 );sound(y2,fs,bits); 无限个回声滤波器程序设计y,fs,bits=wavread();z=zeros(3000,1);y;y=y;zeros(3000,1);a=0.5;R=5000;%滤波器阶数设置，其值越高，回声越明显Bz=0,zeros(1,R-1),1;%无限个回声滤波器系统函数分子Az=1,zeros(1,R-1),-a;%无限个回声滤波器系统函数分母y2=filter(Bz,Az,y);Y2=fft(y2,6001);h,w=freqz(Bz,Az);subplot(3,2,1);plot(abs(h);title(无限个回声滤波器幅频响应);subplot(3,2,2); plot(angle(h);title(无限个回声滤波器相频响应);subplot(3,2,3:4);plot(y2);title(无限个回声滤波器时域图 );subplot(3,2,5);plot(abs(Y2);title(无限个回声滤波器信号幅频 );subplot(3,2,6);plot(angle(Y2);title(无限个回声滤波器信号相频 );sound(y2,fs,bits); 全通结构滤波器程序设计y,fs,bits=wavread(.wav);y=y(:,1);z=zeros(3000,1);y;y=y;zeros(3000,1);a=0.5;R=5000;%滤波器阶数设置，其值越高，回声越明显Bz1=a,zeros(1,R-1),1;%全通滤波器系统函数分子Az1=1,zeros(1,R-1),a;%全通滤波器系统函数分母h,w=freqz(Bz1,Az1);yy2=filter(Bz1,Az1,y);YY2=fft(yy2,6001);subplot(3,2,1);plot(abs(h);title(全通滤波器幅频响应);subplot(3,2,2); plot(angle(h);title(全通滤波器相频响应);subplot(3,2,3:4);plot(yy2);title(全通结构时域图 );subplot(3,2,5);plot(abs(YY2);title(全通结构幅频 );subplot(3,2,6);plot(angle(YY2);title(全通结构相频 );sound(yy2,fs,bits);部分结果展示：5、数字滤波器的设计与实现调用函数mstg产生st，自动绘图显示st的时域波形和幅频特性曲线调用ellipord和ellip分别设计三个椭圆滤波器，并绘图显示其幅频响应特性曲线。调用filter，用三个滤波器分别对信号st进行滤波，分离出三路不同载波频率的调幅信号y1(n)、y2(n)和y3(n)绘图显示y1(n)、y2(n)和y3(n)的时域波形和幅频特性曲线结束5.1设计流程图5.2设计步骤与过程通过观察，三路信号的载波频率分别为250，500，1000Hz,所以低通滤波器的截至频率250Wc500，带通滤波器的截至频率250W1500,500W21000;高通滤波器的截至频率500Wc1000; 调用ellipord计算椭圆DF阶数N和通带截止频率wp. 调用ellip计算椭圆带通DF系统函数系数向量B和A.产生信号N=1600;Fs=10000;T=1/Fs;Tp=N*T; %采样频率Fs=10kHz，Tp为采样时间t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;T=1/Fs;t=0:T:(N-1)*T;k=0:N-1;fc1=Fs/10;fm1=fc1/10;fc2=Fs/20;fm2=fc2/10;fc3=Fs/40;fm3=fc3/10;xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);st=xt1+xt2+xt3;fxt=fft(st,N); %计算信号st的频谱subplot(2,1,1)plot(t,st);grid;xlabel(t/s);ylabel(s(t);axis(0,Tp/8,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(幅度);%三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。低通滤波fp=280;fs=450;%fp必须大于250，fs必须小于500wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;%DF指标（低通滤波器的通、阻带边界频）N,wp0=ellipord(wp,ws,rp,rs);%调用ellipord计算椭圆DF阶数N和通带截止频率wpB,A=ellip(N,rp,rs,wp0);%调用ellip计算椭圆带通DF系统函数系数向量B和Ay1t=filter(B,A,st);%滤波器软件实现figure;subplot(2,1,1);H1,w=freqz(B,A,1000);%低通的频率响应m=abs(H1);plot(w/pi,20*log(m/max(m);grid on;title(低通滤波损耗函数曲线);xlabel(w);ylabel(H1);subplot(2,1,2);plot(t,y1t);title(低通滤波后的波形);xlabel(t);ylabel(y1t);带通滤波Fs=10000;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,wp0=ellipord(wp,ws,rp,rs);B,A=ellip(N,rp,rs,wp0);y2t=filter(B,A,st);figure;subplot(2,1,1);H2,w=freqz(B,A,1000);m=abs(H2);plot(w/pi,20*log(m/max(m);grid on;title(带通滤波损耗函数曲线);xlabel(w);ylabel(H2);subplot(2,1,2);plot(t,y2t);title(带通滤波后的波形);xlabel(t);ylabel(y2t);高通滤波Fs=10000;fp=550;fs=900;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;N,wp0=ellipord(wp,ws,rp,rs);B,A=ellip(N,rp,rs,wp0,high);y3t=filter(B,A,st);figure;subplot(2,1,1);H3,w=freqz(B,A,1000);m=abs(H3);plot(w/pi,20*log(m/max(m);grid on;title(高通滤波损耗函数曲线);xlabel(w);ylabel(H3);subplot(2,1,2);plot(t,y3t);title(高通滤波后的波形);xlabel(t);ylabel(y3t);部分结果展示：6、心电信号的处理6.1设计流程图开始向心电信号中加入白噪声，50HZ工频干扰，带限chirp噪声使用Butterworth型低通数字滤波器对信号过滤使用切比雪夫I型数字低通滤波器对信号过滤使用带阻滤波器对信号过滤进行效果比较结束6.2设计步骤与过程一、设计原理（1）、巴特沃斯低通滤波器巴特沃斯低通滤波器的平方幅度响应为其中，n为滤波器的阶数，为低通滤波器的截止频率。 该滤波器具有 一些特殊的性质： 对所有的n,都有当时, ； 对所有的n,都有当时, ；是的单调递减函数,即不会出现幅度响应的起伏； 当时,巴特沃斯滤波器趋向于理想的低通滤波器； 在处平方幅度响应的各级导数均存在且等于0，因此在该点上取得最大值，且具有最大平坦特性。 图1展示了2阶、4阶、8阶巴特沃斯低通滤波器的幅频特性。可见阶数n越高，其幅频特性越好，低频检测信号保真度越高，过渡带变窄，即衰减加剧，但半功率点不变。图1 巴特沃斯低通滤波器的幅频特性（2） 切比雪夫I型数字低通滤波器（1）确定数字低通滤波器的技术指标：通带截止频率p、通带衰减ap、阻带截止频率s、阻带衰减as切比雪夫滤波器的振幅平方特性如图2所示：图2 切比雪夫滤波器的振幅平方特性(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。如果采用脉冲响不变法，边界频率的转换关系为:如果采用双线性变换法，边界频率的转换关系为(3) 按照模拟低通滤波器的技术指标设计模拟低通滤波器。(4) 利用双线性变换法将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。（5）数字低通技术指标为：p=0.4rad, a p=1dB; s=0.5rad, a s=40Db（6）模拟低通的技术指标为：归一化截止角频率wp=2pi*Fs/Ft; ws=2pi*Fs/Ft(7)利用模拟切比雪夫滤波器设计数字滤波器。通带截止频率为：wp=0.4*pi; 阻带截止频率为：ws=0.5*pi;通带最大衰减为：Rp=1;阻带最大衰减为：As=15;设定周期为1s；模拟低通滤波器的生成：b,a=cheby1(n,1,Wn,low,s);满足设计指标的最小阶数和截止频率：Wnn,Wn=cheb1ord(OmegaP,OmegaS,1,40,s)。最后实现输入输出、幅频特性、相频特性的图形。（3） IIR数字滤波器的性质无限长冲激响应( IIR) 数字滤波器是数字滤器的一种，数字滤波器根据其冲激响应函数的时域特性,还包括有限长冲激响应( FIR)数字滤波器。IIR数字滤波器的特征是：具有无限持续时间冲激响应；需要用递归模型来实现，这可以从其差分方程得出，也可以从其系数函数为：得出。 数字巴特沃思滤波器属于IIR滤波器，该类滤波器具有特定的性质和设计方法。目前比较成熟的IIR数字滤波器设计方法有两种：1）直接法目前所用的方法主要是：零极点累试法、频域幅度平方误差最小法和时域单位脉冲响应逼近法。直接法的最大优点在于可以设计任意幅频特征的滤波器。2)间接法，目前所用的方法主要是：冲激响应不变法、阶跃响应不变法和双线性法。它们都是借助于 已经成熟的现有低通滤波器原型进行设计，即对数字低通数字滤波器，先将数字低通滤波器的技术指标按希望的设计方法转换为模拟低通滤波器的技术指标，再按指定的模拟低通滤波器的类型设计模拟滤波器H(s)，然后，将模拟滤波器的系统函数H(s)从s平面转换到z平面，得到数字低通滤波器的系统函数H(z)；如所设计的数字滤波器为高通、带通或带阻滤波器，则可借助模拟滤波器的频带变换转换为低通模拟滤波器。 由于直接法设计巴特沃思滤波器相对复杂，在不需要任意幅频特征的情况下，一般采用问接法，同时由于冲激响应不变法和脉冲响应不变法，从s平面转换到z平面的映射为多值映射，容易造成频谱混叠，故而本文采用不会产生频谱混叠的双线性变换法。 二、 课程设计步骤及结果分析（1）心电数据的导入将老师给的心电信号原始数据存于桌面，然后我选用桌面上的第一组数据a=load(C:Documents and Settingshnie桌面心电信号数据a01.txt); （2）绘出心电信号的时域图和频谱图将导入的两行数据分别用t,b来替换，然后通过调用plot函数来画出时域图，然后通过对1000个心电数据的幅值进行FFT运算，再次调用plot函数来绘出频域图，具体设计如下：figure(1);subplot(2,1,1);t=a(1:1000,1);b=a(1:1000,2);plot(t,b);title(原始波形图);xlabel(时间（s）);ylabel(幅值（A）);y1=fft(a(:,2),1000);f1=100*(0:999)/1000;subplot(2,1,2);plot(f1,abs(y1);title(原始频谱图);xlabel(频率（Hz）);ylabel(幅度(dB);图通过导入的心电信号数据发现在其频谱图上的020Hz和80100Hz之间的幅值比较大，而在3070Hz之间的幅值相对较小。（3） 加入噪声干扰这里我加入的噪声是：白噪声，50Hz的工频噪声，带限chirp噪声。（4）白噪声通过用s来代表加入白噪声后的信号，并进行数字滤波器的