中南大学数字信号处理课程设计任务书分解

目录

一、课程设计目的 (2)

二、课程设计内容 (2)

三、课程设计要求 (3)

四、程序设计原理及思路

4.1 设计原理 (3)

4.2 设计思路

4.2.1、第一大题 (6)

4.2.2、第二大题 (7)

4.2.3、第三大题

3.1 (7)

3.2 (7)

3.3 (8)

4.3、GUI图形界面 (8)

五、程序测试输出结果 (9)

六、总结 (15)

七、主要参考资料 (15)

八、程序源代码清单 (15)

中南大学

课程设计任务书

一、课程设计目的：

1．全面复习课程所学理论知识，巩固所学知识重点和难点，将理论与实践很好地结合起来。

2．提高综合运用所学知识独立分析和解决问题的能力；

3．熟练使用一种高级语言进行编程实现。

二、课程设计内容

课程设计选题组一：

一、一个连续信号含两个频率分量，经采样得

ππ?=-

x n n N

()=sin(2*0.125*n)+cos(2*(0.125+f)*n),0,1,,1

当N=16，Δf分别为1/16和1/64时，观察其频谱；当N=128时，Δf不变，其结果有何不同，为什么？绘出相应的时域与频域特性曲线，分析说明如何选择DFT参数才能在频谱分析中分辨出两个不同的频率分量。

二、对周期方波信号进行滤波

1）生成一个基频为10Hz的周期方波信号。

2）选择适当的DFT参数，对其进行DFT，分析其频谱特性，并绘出相应曲线。 3）设计一个滤波器，滤除该周期信号中40Hz以后的频率分量，观察滤波前后信号的时域和频域波形变化

4）如果该信号淹没在噪声中，试滤除噪声信号。

三、音乐信号处理：

1）获取一段音乐或语音信号，设计单回声滤波器，实现信号的单回声产生。给出单回声滤波器的单位脉冲响应及幅频特性，给出加入单回声前后的信号频谱。

2）设计多重回声滤波器，实现多重回声效果。给出多回声滤波器的单位脉冲响应及幅频特性，给出加入多重回声后的信号频谱。

3）设计全通混响器，实现自然声音混响效果。给出混响器的单位脉冲响应及幅频特性，给出混响后的信号频谱。

4）设计均衡器，使得不同频率的混合音频信号，通过一个均衡器后，增强或

削减某些频率分量**。（**可选做）

三、课程设计要求

1、使用MATLAB（或其它开发工具）编程实现上述内容，写出课程设计报告。滤波器设计题目应尽量避免使用现成的工具箱函数。为便于分析与观察，设计中所有频谱显示中模拟频率应以实际频率显示，数字频率应对 归一化。

2、课程设计报告的内容包括：

(1)课程设计题目和题目设计要求；

(2)设计思想和系统功能分析；

(3)设计中关键部分的理论分析与计算，参数设置，关键模块的设计思路；

(4)测试数据、测试输出结果，及必要的理论分析和比较

(5)总结，包括设计过程中遇到的问题和解决方法，设计心得与体会等；

(6)参考文献；

(7)程序源代码清单。

3、演示系统使用GUI界面或混合编程实现集成打包发布，

四、程序设计思路及分析

在开始动手之前，首先回顾了数字信号处理的理论知识，特别是关于采样和相关方面的知识。然后熟悉Matlab软件，最后是查找相关声音处理方面的资料，熟悉它的作用功能及研究的意义。

4.1设计原理

1.用DFT对连续信号进行谱分析

（1）原理

（2）频率分辨率与DFT参数的选择

频率分辨率是指所用的算法能将信号中两个靠得很近的谱峰分开的能力。

设是一个带限的连续时间信号，最高频率为f c，根据时域采样定理，采样频率f s>2f c，一般取。

对在时间长度为T p的一段上抽取N点，得到一个长度为N的有限长序

列x(n)，则有

由于f s对应于数字频率，对x(n)作N点DFT，则数字域的频率分辨率

此时，相应的模拟域的频率分辨率为

上式说明：如果保持采样点数N不变，要提高谱的分辨率(F减小)，必须降低采样速率，采样速率的降低会引起谱分析范围减少；如维持f s不变，为提高分辨率可以增加采样点数N。

2．用DFT进行谱分析的误差问题

（1）混叠现象

利用DFT逼近连续时间信号的傅里叶变换，为避免混叠失真，按照抽样定理的要求，采样频率至少是信号最高频率的两倍。

解决混叠问题的唯一方法是保证采样频率足够高。

（2）截断效应

利用DFT处理非时限序列时，须将该序列截断。设序列的频谱为，矩形窗函数的频谱为，则截断后序列的频谱为

由于矩形窗函数频谱的引入，使卷积后的频谱被展宽了，称为频谱泄露（截

断效应）。

减少方法：选择适当形状的窗函数，如汉宁窗或汉明窗等。

（3）栅栏效应

DFT是有限长序列的频谱等间隔采样，相当于透过一个栅栏去观察原来信号的频谱，这种现象称为栅栏效应。

减小栅栏效应的方法：末尾补零。

补零没有对原信号增加任何新的信息，因此不能提高频率分辨率。补零的目的：使数据N为2的整数次幂，以便于用快速傅里叶变换算法（FFT），而且补零还可对原X(k)做插值。

3. 巴特沃斯滤波器

巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。

一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。

巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。

4.FIR数字滤波器

FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，又称为非递归型滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。

在进入FIR滤波器前，首先要将信号通过A/D器件进行模数转换，把模拟信号转化为数字信号；为了使信号处理能够不发生失真，信号的采样速度必须满足奈奎斯特定理，一般取信号频率上限的4-5倍做为采样频率

5.IIR数字滤波器

IIR(Infinite Impulse Response)数字滤波器，又名“无限脉冲响应数字滤波器”，或“递归滤波器”。递归滤波器，也就是IIR数字滤波器，顾名思义，

具有反馈，一般认为具有无限的脉冲响应。

(1)、封闭函数

IIR数字滤波器的系统函数可以写成封闭函数的形式。

(2)、IIR数字滤波器采用递归型结构

IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。

(3)、借助成熟的模拟滤波器的成果

IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。

(4)、需加相位校准网络

IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。

4.2设计思路

1、第一大题

首先要在matlab中实现信号x(n)，然后分别按照Δf=1/16，Δf=1/64进行等间隔采样；再调用函数FFT对它进行离散傅里叶变换；画图调用subplot()和plot()来完成，并且要注意频率的归一化。

由以下得到的波形可以看出:数字时域中每格1/N,频域中不失真时应该在0.125和0.125+Δf的地方有幅值，如当N=16，Δf=1/16时，两频率分量应该是0.125和0.125+1/16，对应到谱线上是2*（1/N）、3*（1/N）即第2和3条，分辨率足够到把两谱线分开，但当Δf=1/64时由于两谱线靠的太近，1/N的分辨率不能分开他们，出现频谱混叠现象。根据奈奎斯特定理，采样频率应该保证大于两倍的截止频率，就能以免发生频谱混叠。所以，Δf=1/16时频谱具有两个清晰的双峰。

相对于N=16时，当N=128时，分辨率等于fs/N，Δf>=1/N时，主瓣宽度变窄，频率分辨率增大，不会发生频谱泄露，图中两谱线都能分清，且时域窗增大了，泄露现象有所改善，可见提高分辨率可以增加信号的截取长度，但旁瓣的相对幅度并不减少。

2、第二大题

首先对待滤波的信号进行频谱分析，观察信号频率分布的规律，从而确定数字滤波器的类型（FIR 滤波器、IIR 滤波器、自适应滤波器、小波滤波器等）。 在加性噪声的情况下，若信号的频谱与噪声的频谱基本不重叠，可以采用频率选择滤波器（FIR 滤波器、IIR 滤波器）。

若信号的频谱与噪声的频谱重叠较多，可以采用自适应滤波、小波滤波等。 若为乘性噪声，可以根据同态滤波的原理对信号进行预处理，然后再按照加性噪声的情况处理。

在确定了数字滤波器类型后，还需要根据信号时域特性、频域特性、或时频特性确定滤波器的设计参数，设计出相应的数字滤波器。

最后，利用该数字滤波器对信号进行滤波，在时域和频域观察信号滤波的主观及客观效果。若主观及客观效果满足要求，说明分析过程和滤波方法正确有效，若不满足要求，需要重新分析和设计。

针对题目要求的10Hz 周期方波频率，选择低通滤波器进行滤波，我优先选择了巴特沃斯低通滤波器，因为其具有平坦的通带与阻带。根据方波时域进行傅里叶变换后得到的频域图，分析正半峰容易得到通带边界频率与阻带截至频率，再由|20*log （输入/输出）|求出衰减，分别得到最大通带衰减与最小阻带衰减，由这四个巴特沃斯参数可以求出N 最大阶数与Wc 截至频率。最后通过butter 求出系统函数多项式的分子与分母。

3、第三大题

3.1 单回声

回音可以由简单的延时单元产生。直达声和在R 抽样周期后出现的一种单个回音，可以用FIR 滤波器产生，微分方程为：y[n]=x[n]+αx[n-R] |α|<1，传输函数为：H(z)=1+αz -R ，传递函数的幅频响应形状象梳子，这种滤波器又叫梳状滤波器。

梳妆滤波器可用于消除频谱中的谐波干扰和其他频谱等间隔分布的干扰。

3.2 多重回声

为了产生以间隔R 个抽样周期分开的具有指数衰减振幅的多重回声，可用一个以下形式传输函数的FIR 滤波器：

R

NR

N R N N R R z z z z z z H ---------=++++=ααααα111)()1(122 无限个振幅以指数衰减间隔为R 个抽样周期的多重回声可用以下形式传输

函数的IIR 滤波器生成：1||11

)(<-=-ααR

z z H IIR 多重回声滤波器的基本频率F R =F s /R ，通常锁定在伴音设备的基频上，比如基鼓拍子。

3.3 全通混响

声音在一个封闭的空间中到达听众时，包含几个部分：直达声、早期的反射和混响。早期的反射由几个空间上邻近的直达声的基本延迟和衰减组成，而混响则由密集的回声组成。上述多重回声滤波器不能提供自然声音混响。由其幅频特性可知，其幅度响应对于所有频率不是常数，收听效果不能令人满意。其次每秒回声数量太少会引起合成声的颤动，需要每秒约1000个回音才能生成没有颤动的反射声。

3.3.1全通滤波器

为了生成更真实的混响，提出一种有全通结构的混响器，传输函数为：

1||1)(<++=--αααR R

z z z H

3.3.2自然声音混响

4、GUI 图形界面

基于matlab 可视化图形用户界面（GUI ），使得界面设计在可视化状态下进行，设计过程简单直观，实现了“所见即所得”。

主要运用菜单编辑器工具，创造下拉式菜单和快捷菜单。在创建某个菜单项中，可以在右边显示某菜单的有关属性。利用菜单编辑器创建了‘第一大题’，‘第二大题’和‘第三大题’三个一级菜单项，并且在各个一级菜单下，创建多个子菜单。

五、程序测试输出结果

图一课程设计主界面

图二第一大题菜单及其子菜单

图三第二大题菜单及其子菜单

图四第三大题菜单及其子菜单

图五第一大题第一小问

问题分析：

观察N=16与128的频谱变化，当N=16时，会对信号产生截断，当N=128后，使主瓣变窄，减小泄露，提高了频率分辨率，但旁瓣的相对幅度并不减小。所以在矩形窗下，适当提高N, 可以降低频谱泄露。

观察f=1/64与1/16的频谱变化，当f=1/16容易发现，两峰之间没有发生频谱混叠。根据采样定理，w=л（对应模拟频率f=Fs/2）附近产生频谱混叠现象。因此，理论上必须Fs>=2fc，以免发生频谱混叠现象。

图六第一大题第二小问

图七第二大题第一小问

图八第二大题第二小问

图九第二大题第三小问

问题分析：

观察方波的单峰频谱，由于该图是低通，有横纵坐标可以读出通带边界频率和阻带截止频率，分别是40与50；幅值通带最大衰减与阻带最小衰减分别是0.7和0.1 。由于巴特沃斯滤波器函数调用格式是buttor（αp，αs，ωp，ωs），其中要对40与50归一化，所以是40/500，50/500；再根据20*|log(输入/输出)|，可以求出3，15 。最后由butter求出系统函数的多项式分子与分母。

图十第三大题第一小问

问题分析：

通过实验得到的频谱图分析发现，多重回声效果较单回声效果更加明显，如果延时调试得当，会得到很好的回声效果，产生类似在音乐厅的效果，自然混响的效果则较前两者更加明显。

图十一第三大题第二小问

图十二第三大题第三小问

六、总结

通过这次数字信号处理的课程设计，对信号分析、matlab、音频处理、GUI 图形界面都有更深层次的理解与应用。

在数字信号中，一开始画幅频图像总是得不到理想的效果，后来经过查找资料，知道了要对原函数进行了fft转换后，应该再归一化并且求出幅值。Matlab 的各种滤波器函数运用与subplot作图也得到了更好的灵活使用，例如梳妆滤波器，加性噪声和做出1：1：2型的图表；对多组程序的布局也有良好的条理分行。在音频处理上，经过网上搜索，我学会了成品音频都是经过多声道的归一融合，原本是存在多方位多乐器与回声的存在；还懂了音频比特，音频频率，回声利用延时单元创造，音频通道通过设置滤波器函数的参数得到。第一次利用GUI图形界面，也遇到了不少瓶颈，例如界面的简洁优化与合理布局，此外对于产生新选项窗口、按钮的触发、音频的播放和减少空白运行时间的问题都尝试了多种手段，最终也得到了较理想的效果。

总之，这次数字信号课程设计让我获益匪浅，不止知识上的丰富完善，还有遇到各种各样的瓶颈时，也可以坚持下来找到解决的方法。这些所学会的一切对我日后信号分析打下良好坚实的基础，对于其他问题也可以淡定分析解决。

七、主要参考资料

[1] Digital Signal Processing: A Computer-Based Approach，Sanjit K. K.

Mitra，2000

[2] Understanding Digital Signal Processing，Richard G. Lyons，科学出

版社

[3] Dual-tone multi-frequency [EB/OL] .https://www.360docs.net/doc/fb13690129.html,/wiki/

Dual-tone_multi-frequency

[4] 《数字信号处理》.丁玉美等西安电子科技大学出版社

[5]数字信号处理教程（第二版），程佩青，清华大学出版社，2001

[6] 数字信号处理，赵树杰等，西电出版社，1997.10

[7] 数字信号处理—时域离散随机信号处理，丁玉美等，西电出版社，2002.11

[8] MATLAB及在电子信息课程中的应用，陈怀琛等，电子工业出版社出版，2002.4

八、程序源代码清单

8.1主程序代码

function varargout = test6(varargin)

% TEST6 M-file for test6.fig

% TEST6, by itself, creates a new TEST6 or raises the existing

% singleton*.

% H = TEST6 returns the handle to a new TEST6 or the handle to

% the existing singleton*.

%

% TEST6('CALLBACK',hObject,eventData,handles,...) calls the local

% function named CALLBACK in TEST6.M with the given input arguments.

%

% TEST6('Property','Value',...) creates a new TEST6 or raises the

% existing singleton*. Starting from the left, property value pairs are % applied to the GUI before test6_OpeningFunction gets called. An

% unrecognized property name or invalid value makes property application % stop. All inputs are passed to test6_OpeningFcn via varargin.

%

% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one % instance to run (singleton)".

%

% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help test6

% Last Modified by GUIDE v2.5 25-Sep-2014 17:07:41

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;

gui_State = struct('gui_Name', mfilename, ...

'gui_Singleton', gui_Singleton, ...

'gui_OpeningFcn', @test6_OpeningFcn, ...

'gui_OutputFcn', @test6_OutputFcn, ...

'gui_LayoutFcn', [] , ...

'gui_Callback', []);

if nargin && ischar(varargin{1})

gui_State.gui_Callback = str2func(varargin{1});

end

if nargout

[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});

else

gui_mainfcn(gui_State, varargin{:});

end

% End initialization code - DO NOT EDIT

% --- Executes just before test6 is made visible.

function test6_OpeningFcn(hObject, eventdata, handles, varargin)

% This function has no output args, see OutputFcn.

% hObject handle to figure

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% varargin command line arguments to test6 (see VARARGIN)

% Choose default command line output for test6

handles.output = hObject;

% Update handles structure

guidata(hObject, handles);

% UIWAIT makes test6 wait for user response (see UIRESUME)

% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line. function varargout = test6_OutputFcn(hObject, eventdata, handles)

% varargout cell array for returning output args (see VARARGOUT);

% hObject handle to figure

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure

varargout{1} = handles.output;

% -------------------------------------------------------------------- function Untitled_1_Callback(hObject, eventdata, handles)

% hObject handle to Untitled_1 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% -------------------------------------------------------------------- function Untitled_2_Callback(hObject, eventdata, handles)

figure(1)

N=16;

f=1/16;

n=0:(N-1);

xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+f)*n);

X16k=fft(xn,16);

subplot(2,2,1);

stem(n,xn,'r');

xlabel('n');

ylabel('x(n)');

title('N=16,f=1/16的x(n)时域曲线');

subplot(2,2,2);

plot(2/N*n(1:N/2),M,'b');

xlabel('\omega/\pi');

ylabel('|X(e^j^\omega)|');

title('N=16,f=1/16的x(n)频域曲线');

f=1/64;

n=0:N-1;

xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+f)*n);

X16k=fft(xn,16);

M=abs(X16k(1:N/2));

subplot(2,2,3);

stem(n,xn,'r');

xlabel('n');

ylabel('x(n)');

title('N=16,f=1/64的x(n)时域曲线');

subplot(2,2,4);

plot(2/N*n(1:N/2),M,'b');

xlabel('\omega/\pi');

ylabel('|X(e^j^\omega)|');

title('N=16,f=1/64的x(n)频域曲线');

% hObject handle to Untitled_2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% -------------------------------------------------------------------- function Untitled_3_Callback(hObject, eventdata, handles)

figure(2)

N=128;

f=1/16;

n=0:(N-1);

xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+f)*n);

X128k=fft(xn,128);

subplot(2,2,1);

stem(n,xn,'r');

xlabel('n');

ylabel('x(n)');

title('N=128,f=1/16的x(n)时域曲线');

subplot(2,2,2);

plot(2/N*n(1:N/2),M,'b');

xlabel('\omega/\pi');

ylabel('|X(e^j^\omega)|');

title('N=128,f=1/16的x(n)频域曲线');

f=1/64;

n=0:N-1;

xn=sin(2*pi*0.125*n)+cos(2*pi*(0.125+f)*n);

X128k=fft(xn,128);

M=abs(X128k(1:N/2));

subplot(2,2,3);

stem(n,xn,'r');

xlabel('n');

ylabel('x(n)');

title('N=128,f=1/64的x(n)时域曲线');

subplot(2,2,4);

plot(2/N*n(1:N/2),M,'b');

xlabel('\omega/\pi');

ylabel('|X(e^j^\omega)|');

title('N=128,f=1/64的x(n)频域曲线');

% hObject handle to Untitled_3 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% --------------------------------------------------------------------

function Untitled_7_Callback(hObject, eventdata, handles)

figure(3)

Fs=1000;

T=1/10;

t=0:1/Fs:1;

x=square(2*pi*10*t);

subplot(211);

plot(t,x,'b');

title('基频为10Hz的周期方波信号时域');

xlabel('t/s');

ylabel('x的幅值');

axis([0,0.4,-1.5,1.5]);

N=T*Fs;

t1=0:1/Fs:1/Fs*(N-1);

x1=square(2*pi*10*t1);

f=(-N/2)*Fs/N:Fs/N:((N/2-1))*Fs/N;

fx=fftshift(fft(x1,N))/N;

subplot(212);

plot(f,abs(fx));

title('10Hz方波信号的频谱');

xlabel('\omega/\pi');

ylabel('|X(e^j^\omega)|');

% hObject handle to Untitled_7 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA)

% -------------------------------------------------------------------- function Untitled_10_Callback(hObject, eventdata, handles)

figure(4)

Fs=1000;

T=1/10;

t=0:1/Fs:1;

x=square(2*pi*10*t);

N=T*Fs;

t1=0:1/Fs:1/Fs*(N-1);

x1=square(2*pi*10*t1);

f=(-N/2)*Fs/N:Fs/N:((N/2-1))*Fs/N;

fx=fftshift(fft(x1,N))/N;

[n,Wn]=buttord(40/500,50/500,3,15);

[b,a]=butter(n,Wn);

y=filter(b,a,x);

subplot(211);

plot(t,y);

title('10Hz的周期方波通过低通滤波器的时域');

xlabel('t/s');

ylabel('x的幅值');