数字信号处理课程设计说明书.doc

数字信号处理课程设计说明书

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

个人收集整理 勿做商业用途
1引言
随着信息时代和数字世界的到来，数字信号处理已成为一门极其重要的学科和技术领域。数字信号处理在语音通信,图像，自动控制,雷达,军事，航空航天，医疗和家用电器等众多领域得到了广泛的应用。其中MATLAB是数字信号处理中极其重要的软件.
快速傅里叶变换（FFT)作为计算和分析工具，在众多学科领域（如信号处理、图像处理、生物信息学、计算物理、应用数学等)有着广泛的应用.在高速数字信号处理领域，如雷达信号处理，FFT的处理速度往往是整个系统设计性能的关键所在。FFT算法的基本思想就是利用权函数的周期性、对称性、特殊性及周期N的可互换性，将较长序列的DFT运算逐次分解为较短序列的DFT运算。
快速傅里叶变化的特点:
（1)复数运算：傅立叶变换是基于复数的，因此首先知道复数的运算规则,在FFT算法中，只涉及复数的加、减和乘法三种运 。 
（2)蝶形变换：普通的FFT算法称为基2的FFT算法，这种算法的核心是蝶形变换。  
(3)w数组。
（4）复数数组排序,在基2的蝶形变换中，复数数组需要重新排.
本次课设内容：录制一段个人的语音信号,并对其进行采样；画出采样的后的时域波形和频谱图；在Matlab环境下编写基2DIT-FFT算法；利用自己编写好的算法对已采集的语音信号进行频谱分析，并画出语音信号的时域和频谱图,并与Matlab数字信号处理工具箱中的fft函数进行 对比研究来验证自编算法的正确性。最后设计一个信号处理界面，可以实现对输入信号的选择和不同点的FFT运算的选择。
本次课程设计目的:初步掌握了MATLAB语言的主要特点和作用、MATLAB的GUI设计及信号的基本运算的实现。全面复习课程所学理论知识,巩固所学知识重点和难点,根据课堂讲授内容，将理论与实践很好地结合起来，消化课堂所讲解的内容;通过调试积累经验，逐渐培养学生的编程能力、用计算机解决实际问题的能力。提高综合运用所学知识独立分析和解决问题的能力；熟练使用一种高级语言进行编程实现;熟悉Matlab的工作环境及运行状况。

2 MATLAB
2.1 MATLAB简介
MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，代表了当今国际科学计算软件的先进水平.主要应用于工程计算、信号处理与通讯、图像处理、信号检测、金融建模设计与分析、控制系统设计以及计算生物学等众多应用领域。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多,并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件.MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连MATLAB开发工作界面接其他编程语言的程序等,主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 　　
MATLAB特点：
（1）高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；
（2)具有完备的图形处理功能，实现计算结果和编程的可视化;
（3）友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握;
（4)功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等),为用户提供了大量方便实用的处理工具。
MATLAB具有方便的数据可视化功能，以将向量和矩阵用图形表现出来。高层次的作图包括二维和三维的可视化、图象处理、动画和表达式作图。可用于科学计算和工程绘图。MATLAB不仅具有高层绘图能力，而且还具有底层绘图能力-句柄绘图方法。使用户可以用来开发各专业的专业图形。图形用户界面(GUI）是一种提供人机交互的工具和方法.
2.2 MATLAB基本语法
（1）语言简洁紧凑，使用方便灵活，库函数及其丰富.程序书写形式自由,利用丰富的库函数避开繁杂的子程序编写任务，压缩了一切不必要的编程工作。
（2)运算符丰富，由于MATLAB使用C语言编写的,它提供了和C语言几乎一样多的运算符,灵活使用MATLAB的运算符将使程序变得极为简短.
(3）程序限制不严格，程序设计自由度大。例如，在MATLAB里,用户无需对矩阵预定义就可使用。
（4）程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行.
（5）MATLAB的图形功能强大，在FORTRAN和C语言里,绘图都很不容易，但在MATLAB里，数据的可视化非常简单。它还具有较强的编辑图形界面的能力。
（6)MATLAB的缺点是：它和其他高级程序相比，程序的执行速度较慢。由于MATLAB的程序不用编辑等预处理，也不生成可执行文件，程序为解释执行，所以速度较慢。
（7）功能强大的工具箱是MATLAB的另一特色，它包含两个部分：核心部分和各种可选的工具箱。核心部分中有数百个核心内部函数。其工具箱又分为两类：功能性工具箱和学科性工具箱。

3快速傅里叶变换（FFT）算法
3。1 FFT的分类
快速傅里叶变换（FFT）是为提高DFT运算速度而采用的一种算法，快速算法的种类较多，最基本的算法是时域2FFT（DIT-FFT）和频域基2（DIT-FFT），本课程设计在MATALB环境下编写基2 DIT-FFT算法。
3.2 FFT的基本算法分析
3。2.1 FFT的基本算法
快速傅里叶变换(FFT）是为提高DFT运算速度而采用的一种算法。
对一个有限长度序列x（n)的N点的DFT为：X(k）=∑x（n）W^knN （k=0，1，……，N—1;n=0,1,……，N—1;W=e^-j2π/N)
当N=4时，X（k）可展开为：
X(0）= x(0）W^0＊4+ x（1）W^0＊4 +x（2）W^0＊4+ x(3）W^0＊4
X（1)= x(0）W^0*4+ x(1)W^1*4 +x（2)W^2＊4+ x（3)W^3＊4
X（2）= x（0)W^0＊4+ x(1）W^2*4 +x（2）W^4＊4+ x（3）W^6*4
X（3）= x（0）W^0*4+ x(1）W^3＊4 +x(2)W^6*4+ x（3）W^9＊4
从上式可以看出,要求4点的DFT,需要16次的复数乘法运算，12次复数乘法运算算。由此类推，要求出N点的DFT，需要N^2次复数乘法运算,N＊(N-1）次复数加法运算。当N值较大时,要完成的复数乘法运算和复数加法运算得次数都非常多，无论是用通用计算机还是用DSP芯片,都需要消耗大量的时间，不适合于对实时处理要求高的场合。为了能实时处理DFT，要想减少DFT的运算量可以有两个途径:第一是降N，N的值减小了,运算量就减少了;第二是利用旋转因子的周期性和对称性，可约性。利用这两个途径实现DFT的快速傅里叶变换（FFT）,FFT算法基本上可分为时域抽取法和频域抽取法.
W=e^—j2π/N的性质：
（1）周期性
（2）共轭对称性
(3）可约性
本程序是用基2的按时间抽取的FFT算法（DIT—FFT），设序列x（n)的长度为N，且N满足N=2^M,M为正整数。若N不能满足上述关系,可以将序列x（n）补零实现，则x（n）的N点DFT为：
X（k)=∑x（n）W^knN （k=0,1，……，N—1；n=0,1,……，N—1；W=e^—j2π/N）将n分为奇数与偶数两部分。
按时间抽取基2-FFT算法的基本思路是将N点序列按时间下标的奇偶分为两个N/2点序列，计算这两个N/2点序列的N/2点DFT，计算量可减小约一半；每一个N/2点序列按照同样的划分原则，可以划分为两个N/4点序列，最后,将原序列划分为多个2点序列,将计算量大大降低。
按时间下标的奇偶将N点x（n）分别抽取组成两个N/2点序列,分别记为x1(n)和x2(n），将x（n）的DFT转化为x1(n）和x2（n）的DFT的计算。

用蝶形运算可表式为如图3—1所示：
X(k)
X(k)
X(k+)=X(k)-WX(k)
X(k)=X(k)+WX(k)

W
图3—1 蝶形运算
以此类推，还可以把x1（n)和x2(n）按n值得奇偶分为两个序列，这样就达到了降N得目的，从而减少了运算量。

FFT对DFT的数学运算量改进：
直接采用DFT进行计算，运算量为N^2次复数乘法和N*（N—1）次复数乘法.
当采用M次FFT时，由N=2^M求得M=logN,运算流图有M级蝶形，每一级都由N/2个蝶形运算构成,这样每一级蝶形运算都需要N/2次复数乘法和N次复数加法.M级运算共需要复数乘法次数为C=N/2*M，复数加法次数为C=N＊M。
当N值较大时，FFT减少运算量的特点表现的越明显。
3。2.2基2 DIT-FFT算法程序流程图
第L级中,每个蝶形运算的两个输入相距同一旋转因子对应的蝶形运算相隔点.编程思想：（1)从输入端开始逐级进行，共进行M即运算；（2）在进行第L级运算时，依次要求出个旋转因子，然后计算每个旋转因子所对应的2个蝶形运算.外层控制第L级顺序运算；中间层控制不同种的旋转因子；内层控制同种旋转因子对应的蝶形运算.

图3-2 DIT—FFT运算程序框图
3.2。3 FFT的运算规律
（1）原位运算
①N=2^M的FFT共M级运算，每级有N/2蝶形原位计算，当数据输入到存储器以后，每一组蝶形运算后,结果仍然存放在这同一组存储器的同一位置,不需要另辟存储空间,直接最后输出.
②同一级的蝶形运算每个蝶形运算的输入数据对其他级输入没有影响.
（2）倒序运算的规律
输入序列先按自然顺序存入存储单元,然后经变址运算来实现倒位序排列,用J表示倒序的十进制数，对N=2^M，M位的二进制数从左到右各位数权值位N/2，N/4，N/8……2,1。因此，最高位加1相当于J+N/2。
①.如果最高位为0，则直接得到下一个倒序值，J+N/2;
②。如果最高位为1，则最高位为0（J-N/2），次高位加1（J+N/4）。
③。以此类推，直到最后一位二进制数字。
例如 ，N=8时如表3—1：
表3—1 码位倒序(N=8）
顺序
倒序
十进制
二进制
二进制
十进制
0
000
000
0
1
001
100
4
2
010
010
2
3
011
110
6
4
100
001
1
5
101
101
5
∑6
110
011
3
7
111
111
7
3.2。4 DIT-FFT算法与直接计算DFT运算量的比较
当N=2^M时,其运算流图有M级蝶形，每一级都有N/2个蝶形运算构成。每一个蝶形运算需要1次复数乘和2次复数加.所以每一级运算都需要N/2次复数乘和N次复数加。
M级运算总共需要的复数乘法次数为：

M级运算总共需要的复数加法次数为:

如直接计算DFT，复数乘法为次，复数加法N*（N-1）次。
当N=1024时用FFT需计算5120次复数乘法和10240次复数加法，与直接计算相比运算量大大降低.

4功能的实现
4。1程序所需实现的功能
首先对语音信号进行采样，得到原始语音信号，再对原始语音信号进行抽样，得到N点的序列，再对其进行N点FFT和编写的N点变换,在实验的相应步骤中输出波形，并将N点FFT和编写的N点变换得到的波形进行比较,以验证编写程序的正确性。
4.2具体功能的实现
4.2.1对语音信号进行采样
信号必须为wav格式,并运用wavread对其进行采样。
采样形式为：［y,fs］=wavread（’qing1。wav'）;
其中采样频率及相应音频信号的幅值分别记录在y和fs中。
用freqz函数对已采集的信号y求得频率响应。
形式为：［h,f]=freqz（y）；
其中采样点频率及相应频响值分别记录在f和h中。
然后用plot将上述图像画出。结果如图4-1所示:

图4-1 语音信号的采集
4.2。2输入N的值并进行校正
因为在进行蝶形变化时，N的值必须为2的整数次幂，所以必须对输入的N值进行校正，此处利用函数nextpow2,给N赋值为比其稍大的2的M次幂。即2^I〉=N，M=min（I）。
形式为：
N=input（’请输入N的值：'）；
M=nextpow2 （N); ％ 2M为大于等于N的最小的二的整数次幂的数字
N=2^M;
4.2。3对采样后的信号取前N点
对采样后的信号进行抽样，即取前N个点，并赋值给新的序列x,
形式为：
n=[1：N]；
x=y（n）；
用freqz函数对已抽样的信号x求得频率响应。
形式为：［hh，ff]=freqz（x）；
并用plot对其进行输出：例如输入的N为1024，结果如图4—2所示:

图4-2信号的抽样
4。2.4对抽样后的序列进行一维快速傅立叶变换（FFT）
对抽样后的序列x进行一维快速傅立叶变换（FFT）
形式为：fh=fft（x，N);结果如图4—4示。
4。2。5对抽样后的序列进行倒序
利用函数dec2bin将十进制转换为二进制；
利用函数fliplr将二进制进行倒序；
利用函数bin2dec将二进制转换为十进制。
由于MATLAB中的输出是以1开始的，所以实现倒序的程序为：
p=［0：N-1］;
bN=zeros(N—1，M）;
dbN=zeros(N—1,M）；
bN=dec2bin(p)；
dbN=fliplr（bN）;
i=bin2dec（dbN）+1；
P=［x(i)]；
其中:bN为p的二进制，dbN为bN的倒序，i为dbN的十进制，新序列P为n点序列x倒序后的序列。
4.2.6对倒序后的序列进行蝶形变换
先计算旋转因子：=
蝶形运算是实现FFT的重要方法，其原理图如图4-3所示：
X(k)
X(k)
X(k+)=X(k)-WX(k)
X(k)=X(k)+WX(k)

W
图4—3 蝶形运算表示
两个N/2点DFT 和在蝶形运算图左边，乘以的旋转因子用箭头及系数表示，蝶形运算图右边是它们的运算组合，上面表示求和运算,下面表示求差运算.
计算时从最低级逐级计算，每计算出一个旋转因子,就进行用到此旋转因子的所有蝶形运算。

程序为：
for L=1:M ％L为级数
WN=exp（（—j*2＊pi)/(2^L））;
for Q=1:(2^(L—1）) %Q为第L级的旋转因子个数
WN1=WN^(Q—1）；
for k=Q：2^L：(N—2^（L-1）)
kp=k+2^（L-1）;
t=P（kp）＊WN1; %碟形运算的乘积项
P（kp）=P（k)-t；
P（k）=P（k）+t；
end
end
end
最后把P赋值给一新序列X,输出X.
运行结果如图4-4示：

图4—4结果比较

5 总结
经过这次实习，使我对MATLAB这个软件有了新的更高的认识。平时的实验虽然也用这软件，不过老师已经给出大部分的参考程序,自己稍微改一下就可出来实验结果，与之前相比，这次课程设计更能锻炼学生的自主学习能力和思维能力,使我们的视界更开阔。还有非常重要的一点：MATLAB中Help中的信息,很全面，很具体，有使用技巧，有范例。
这次课程设计的主要内容是录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后的语音信号波形的时域波形和频谱图；在MATLAB环境下编写基于2DIT—FFT算法，利用自己编写的算法对已采集的语音信号进行频谱分析，画出语音信号的时域图和频谱图，并与MATLAB数字信号处理工具箱中的fft函数进行对比研究，验证自编算法的正确性。
通过这次的数字信号处理的课程设计，对METLAB编程语言的基本概念和语法等有了更深一步的理解，使课堂上的知识能应用到实践中。更进步掌握了蝶形运算的实现过程。虽然在一开始编程时总是出错，遇到的问题也很多,通过上网或请教老师同学，最终会把问题解决掉.课程设计不仅仅要的是最后的结果，最重要的是过程。在编程和做界面是好多不懂得,这就要靠自己的自学能力了.因此，在大学应多锻炼自学能力，遇到困难学会自己解决。在做课设时离不开他人的帮助，无论在以后的学习还是工作都要学会与他人和作.
总的来说这次课程设计收获还是挺大的,让我意识到仅仅学会课堂上的知识是不够的.激励我们课下要学会利用时间对自己感兴趣的内容多了解一些.对自己以后在找工作或社会上生存会有很大帮助的。
最后感谢实习以来老师的细心教诲，使我学到了很多的实用的知识。

参考文献
[1］ 范寿康 DSP技术与DSP芯片．北京:电子工业出版社
［2］ 万永格 数字信号处理的MATLAB实现．北京：科学出版社,2007
[3］ 程佩青 数字信号处理教程．北京:清华大学出版社,2001
［4］ 高西全 丁玉美等 数字信号处理．北京:电子工业出版社,2009
附录
％ 语音信号y
[y，fs]=wavread(’qing1。wav')； %fs是采样频率
t=（1:length(y））/fs;
figure（1） ％第一张图
subplot（2，1，1)；
plot(t,y)；
title（'语音信号y的时域波形’）；
xlabel(’时间');
ylabel（'幅度’）;
% 得到y的频率响应h
［h,f]=freqz(y); %f为采样点，h为频率响应
subplot（2，1，2）；
plot((2*pi*f/fs）,abs(h))；
title（’语音信号y的频谱图'）;
xlabel（'频率/Hz’）;
ylabel（'幅度'）；
％ 输入N的值
N=input(’请输入N的值：’)；
M=nextpow2 （N）；
N=2^M；
％取前N个点
n=［1:N];
if N>length（y)
y=[y’,zeros（1,(N—length（y)）)］;
end
x=y(n)；
figure(2) ％第二张图
subplot(2，1，1);
plot（n,x)；
title(’采样得到的信号y（n)的时域波形');
xlabel（’n’）；
ylabel('幅度’）；
[hh,ff]=freqz（x);
subplot(2,1，2)；
plot(（2＊pi*ff/fs),abs(hh））;
title('采样得到的信号y(n)的频谱图’）；
xlabel('频率/Hz')；
ylabel(’幅度’）；
％ 得到y的fft
fh=fft（x，N)；
figure（3) %第三张图
subplot（2，1，1）；
plot(n,abs（fh）)；
title('信号y(n)的N点fft'）；
xlabel（'n’)；
ylabel(’幅度'）；
％计算总级数M
M=log2（N);
% 实现倒序
p=［0:N-1］；
bN=zeros(N-1，M）;
dbN=zeros(N—1,M）；
bN=dec2bin（p)； %bN为p的二进制
dbN=fliplr（bN)； %dbN为bN的倒序
i=bin2dec(dbN）+1; %i为dbN的十进制
P=［x(i)］; ％P为最初的x的n点
％ 计算旋转因子
for L=1:M ％L为级数
WN=exp（(—j＊2*pi）/（2^L));
for Q=1:（2^（L—1)) %Q为第L级的旋转因子个数
WN1=WN^(Q-1）;
for k=Q：2^L:（N-2^(L-1）)
kp=k+2^（L—1)；
t=P(kp)*WN1; %碟形运算的乘积项
P（kp)=P（k)—t;
P(k）=P(k)+t;
end
end
end
X=P; ％X用来储存最后的结果
subplot(2，1，2);
plot(n,abs(X))；
title(’信号y（n）经编写程序的N点运算波形’）；
xlabel（’n'）；
ylabel（'幅度’）;

17