《数字信号处理原理及实现》课程小结

第一篇：《数字信号处理原理及实现》课程小结

时间过得好快，转眼半学期结束了。这半学期数字信号的学习让我受益匪浅。前两章和信号与线性系统相关，介绍了离散时间信号与系统的时域分析方法最深刻的是采样，时域采样定理与采样恢复和采样内插公式。第二章介绍了离散时间信号与系统的频域分析，DTFT与Z变换，系统函数的零极点分布。第三章主要讲了离散傅里叶变换DFT及其性质，和频域采样定理。第四章介绍了傅里叶快速变换FFT，熟悉了其原理特点及方法。五六两章分别介绍了IIR和FIR滤波器，知道了IIR的脉冲响应不变法与双线性变换法及其优缺点，并学会其MATLAB应用设计滤波器，FIR的窗函数法与频域采样法设计滤波器及其MATLAB实现。第七章主要介绍了IIR和FIR滤波器的基本网络结构，通过老师上课习题的练习基本掌握了其结构图的画法。

先说说对课程的建议吧，张晓光老师是个很负责讲课思路也很清晰的老师，知道从学生的角度来讲问题，根据学生的反应来调整课程进度。我们都很喜欢这样的老师，老师开新课之前总是列提纲复习上节课讲的知识，每章结束都根据章节的重要性开一节总结课，这种方式个人觉得很好，希望老师坚持。但是，感觉老师讲题讲的不是很多，或许是课时原因，但我觉得每章结束后开一节例题课，把知识点融进去，毕竟大学生现在做题比较少，这样强制一下效果会更好。这次考试的试题觉得有不少都见过，有的是课后题，但做起来还是有点吃力，应该就是习题练的少，计算跟不上去。至于教材，我觉得编的很好，每章都有相关的MATLAB编程方法，在原理讲清之后就来实践，免去了学生盲目做实验，提高了效率。还有就是老师也很重视实验，总是把相关的MATLAB语句语义讲解清楚，这样我们在编程序时也就相对容易点。但好像老师讲程序时都注重程序的意思了，希望老师以后再讲程序时把它先部分后整体，就是在讲完程序意思后把程序设计思路或框架结构，及各部分要实现什么再讲讲，这样有助于学生设计时设计思路更清晰。再说说考试，老师分卷面成绩和实验成绩及平时成绩，将实验单独考试，可见对实验的重视，也说明MATLAB的重要性，这样确实提高了学生的重视心理，虽然实验做完了，但做完50道题并看完相关讲解，我又收获了不少，理清了设计方法与思路，所以我觉的考试方式还是挺不错的，锻炼了我们各方面的知识。

数字信号课程结束了，真希望您还能教我们别的课。

小组成员：陈文斌、李亚伟、王猛、汪子雄、吴官宝

第二篇：MATLAB实现数字信号处理

数字信号处理

说 明 书

目录

一．摘要…………………………………3 二．课程设计目的………………………3 三．设计内容……………………………3 四．设计原理……………………………4 4.1．语音信号的采集…………………………….4 4.2．滤波器……………………………………….4 4．21．IIR滤波器原理…………………………………….4 4.22．FIR滤波器原理………………………………………5 五．设计步骤……………………………6 5.1录制女音………………………………………6 5.2采样语音信号并画出时域波形和频谱图……7 5.3采用双线性变换法设计IIR滤波器…………10 5.4窗函数法设计FFR滤波器………………......12 5.5用IIR滤波器对信号进行滤波………………14 5.6用FIR滤波器对信号进行滤波………………16 5.7男女声语音信号频谱特点分析………………19 5.8有背景噪声的信号分析………………………20 六．心得体会…………………………….22 七．参考文献…………………………….23

一．摘要：

这次课程设计的主要目的是综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB或者ＤＳＰ开发系统作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。通过对声音的采样，将声音采样后的频谱与滤波。

MATLAB全称是Matrix Laboratory，是一种功能强大、效率高、交互性好的数值和可视化计算机高级语言，它将数值分析、矩阵运算、信号处理和图形显示有机地融合为一体，形成了一个极其方便、用户界面友好的操作环境。经过多年的发展，已经发展成为一种功能全面的软件，几乎可以解决科学计算中所有问题。MATLAB软件还提供了非常广泛和灵活的用于处理数据集的数组运算功能。

在本次课程设计中，主要通过MATLAB来编程对语音信号处理与滤波，设计滤波器来处理数字信号并对其进行分析。

二．课程设计目的：

综合运用本课程的理论知识进行频谱分析以及滤波器设计，通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

三．设计内容：

内容：录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；给定滤波器的性能指标，采用窗函数 法和双线性变换法设计滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；回放语音信号；换一个与你性别相异的人录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点；再录制一段同样长时间的背景噪声叠加到你的语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除。

四．设计原理：

4.1．语音信号的采集

熟悉并掌握MATLAB中有关声音（wave）录制、播放、存储和读取的函数，在MATLAB环境中，有关声音的函数有：

a：y=wavrecord(N,fs,Dtype);利用系统音频输入设备录音，以fs为采样频率，默认值为11025，即以11025HZ进行采样。Dtype为采样数据的存储格式，用字符串指定，可以是：‘double’、‘single’、’int16’、‘int8’其中只有int8是采用8位精度进行采样，其它三种都是16位采样结果转换为指定的MATLAB数据；

b：wavplay(y,fs)；利用系统音频输出设备播放，以fs为播放频率，播放语音信号y；

c：wavwrite((y,fs,wavfile);创建音频文件； d：y=wavread(file);读取音频文件；

关于声音的函数还有sound()；soundsc()；等。4.2滤波器： 4．21.IIR滤波器原理

冲激响应不变法是使数字滤波器在时域上模拟滤波器，但是它们的缺点是产生频率响应的混叠失真，这是由于从s平面到z平面是多值的映射关系所造成的。

双线性变换法是使数字滤波器的频率响应与模拟滤波器的频率响应相似的一种变换方法。为了克服多值映射这一缺点，我们首先把整个s平面压缩变换到某一中介的s1平面的一条横带里，再通过变换关系将此横带变换到整个z平面上去，这样就使得s平面与z平面是一一对应的关系，消除了多值变换性，也 就消除了频谱混叠现象。

双线性法设计IIR数字滤波器的步骤：

1）将数字滤波器的频率指标{ k}由Wk＝（2/T）*tan(wk),转换为模拟滤波器的频率指标{k}.2)由模拟滤波器的指标设计H(s).3)由H(s)转换为H(z)21z1H(z)H(s)sT1z1

4.22.FIR滤波器原理

FIR滤波器与IIR滤波器特点不同，IIR滤波器的相位是非线性的，若需线性相位则要采用全通网络进行相位校正。而有限长单位冲激响应（FIR）数字滤波器就可以做成具有严格的线性相位，同时又可以具有任意的幅度特性。

由于FIR系统的冲激响应就是其系统函数各次项的系数，所以设计FIR滤波器的方法之一可以从时域出发，截取有限长的一段冲激响应作为H（z）的系数，冲激响应长度N就是系统函数H（z）的阶数。只要N足够长，截取的方法合理，总能满足频域的要求。这种时域设计、频域检验的方法一般要反复几个回合，不像IIR DF设计靠解析公式一次计算成功。给出的理想滤波器频率响应是，它是w的周期函数，周期

由傅立叶反变换导出，即

hd(n)1Hd(ejw)ejwndw2，再将hd(n)与窗函数，因此可展开成傅氏级数w(n)相乘就可以得到h(n)。、的计算可采用傅氏变换的现成公式和程序，窗函数也是现成的。但整个设计过程不能一次完成，因为窗口类型和大小的选择没有解析公式可一次算，整个设计可用计算机编程来做。

窗函数的傅式变换W(ejω)的主瓣决定了H(ejω)过渡带宽。W(ejω)的旁瓣大小和多少决定了H(ejω)在通带和阻带范围内波动幅度，常用的几种窗函数有：

矩形窗

w(n)=RN(n)；

Hanning窗

；

Hamming窗

；

Blackmen窗

；

Kaiser窗。

式中Io(x)为零阶贝塞尔函数。

五．设计步骤：

5.1录制女音：

利用MATLAB中的函数录制声音。function nvyin()fs=11025;

%采样频率

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

开始录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);y=wavrecord(3*fs,fs,'double');

%录制声音3秒

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放录音');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];5 disp(str);wavplay(y,fs);

%播放录音

str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);disp('

播放录音结束');str=['@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'];disp(str);wavwrite(y,fs,'原女音');

%声音的存储

5.2采样语音信号并画出时域波形和频谱图

读取语音信号，画出其时域波形和频谱图，与截取后的语音信号的时域波形和频谱图比较，观察其变化。程序如下：

[x,fs,bits]=wavread('女音.wav');

%读取声音

N=length(x);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半 figure(1);subplot(2,2,1);

%把画图区域划分为2行2列，指定第一个图 plot(t, x);

%画出声音采样后的时域波形 title('原女音信号的时域波形');

%给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid;

%添加网格

y=fft(x);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,2,3);

%把画图区域划分为2行2列，指定第三个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

%添加网格

subplot(2,2,4);

%把画图区域划分为2行2列，指定第四个图 if y~=0

%判断指数是否为0

plot(k,20*log10(abs(y(n))));

%画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid;

%添加网格

%实际发出声音落后录制动作半拍的现象的解决 siz=wavread('女音.wav','size');x1=wavread('女音.wav',[3500 32076]);

%截取语音信号 N=length(x1);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

figure(2);subplot(2,2,1);

%把画图区域划分为2行2列，指定第一个图 plot(t,x1);

%画出声音采样后的时域波形 title('截取后女音信号的时域波形');

%给图形加注标签说明 xlabel('时间/t');ylabel('振幅/A');grid;

%添加网格

y1=fft(x1);

%对信号做N点FFT变换

subplot(2,2,3);

%把画图区域划分为2行2列，指定第三个图 k=(n-1)*f0;

%频域采样点

plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('FFT变换后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

%添加网格

subplot(2,2,4);

%把画图区域划分为1行2列，指定第二个图 if y1~=0

%判断指数是否为0

plot(k,20*log10(abs(y1(n))));

%画信号频谱的分贝图 end xlabel('Hz');ylabel('振幅/分贝');title('FFT变换后声音的频谱特性');grid;

%添加网格

原女音信号的时域波形10.5A/幅0振-0.5-10123时间/tFFT变换后声音的频谱特性FFT变换后声音的频谱特性30050A200贝/值分/幅0幅100振00200040006000-***频率/HzHz 截取后女音信号的时域波形10.5振幅/A0-0.5-10123FFT变换后声音的频谱特性50时间/tFFT变换后声音的频谱特性300200振幅/分贝幅值/A01000020004000频率/Hz6000-5002000Hz40006000

结果分析：

由原女音信号的时域波形可知录取开始时实际发出声音大概落后3500个采样点，我们把前3500点去除即可解决实际发出声音落后录制动作半拍的现象。由原女音的的频谱图和截取后声音的频谱图可看出，对声音的截取并不会影响它们频谱分布。

5.3采用双线性变换法设计IIR滤波器：

人的声音频率一般在（1~~4）kHZ之间，则我们只需要设计一个带通滤波器即可滤去声音频带以外的无用噪声，得到比较清晰的声音。根据声音的频谱图分析，设计一个带通滤波器性能指标如下：

fp1＝1000 Hz，fp2＝3000 Hz，fsc1＝500 Hz，fsc2＝3500Hz，As＝100dB，Ap＝1dB，fs=10000 程序如下：

%iir带通的代码： %w=2*pi*f/fs Ap=1;

%通带波纹系数

Az=100;

%最小阻带衰减

wp=[0.2 0.6];

%归一化通带数字截止频率 wz=[0.1 0.7];

%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估计契比雪夫I型滤波器阶数 [b,a]=cheby1(N,Ap,wn);

%N指定滤波器阶数，wn归一化

截 %止频率，Ap通带波动

[h,w]=freqz(b,a);

%求数字滤波器的复频率响应 figure(1);subplot(2,1,1);plot(w/pi,abs(h));

%绘制数字滤波器的频谱图 grid;xlabel('omega/pi');ylabel('振幅（幅值）');title('契比雪夫Ⅰ型带通滤波器的幅频响应');subplot(2,1,2);if abs(h)~=0

%判断指数是否为0

plot(w/pi,20*log10(abs(h)));

%绘制数字滤波器频谱的分贝图 end grid;xlabel('omega/pi');ylabel('振幅(分贝)');title('契比雪夫Ⅰ型带通滤波器的幅频响应');契比雪夫Ⅰ型带通滤波器的幅频响应1振幅（幅值）0.5000.10.20.50.60.70.8/契比雪夫Ⅰ型带通滤波器的幅频响应0.30.40.910振幅(分贝)-200-400-60000.10.20.30.40.5/0.60.70.80.91

5.4窗函数法设计FFR滤波器

线性相位FIR滤波器通常采用窗函数法设计。窗函数法设计FIR滤波器的基本思想是：根据给定的滤波器技术指标，选择滤波器长度N和窗函数ω（n），使其具有最窄宽度的主瓣和最小的旁瓣。其核心是从给定的频率特性，通过加窗确定有限长单位脉冲响应序列h(n)。工程中常用的窗函数共有6种，即矩形窗、巴特利特（Bartlett）窗、汉宁（Hanning）窗、汉明（Hamming）窗、布莱克曼（Blackman）窗和凯泽（Kaiser）窗。

这次设计我采用的是布莱克曼来设计给定数字带通滤波器的参数如下： wp1=0.3pi, wp2=0.6pi, wz1=0.2pi, wz2=0.7pi, Ap=1dB, Az=70dB 程序如下：

Ap=1;

%通带波纹系数 Az=100;

%最小阻带衰减 fs=10000;

%采样频率 wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));

%---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW);

%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);

%---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);%理想带通滤波器 wdWindow=blackman(N);

%布拉克曼窗 hr=wdWindow.*hdWindow';

%点乘

n=0:N-1;

%阶数 subplot(2,2,1);stem(n,wdWindow);

%绘制布拉克曼窗时域波形 xlabel('时间');ylabel('振幅');title('布拉克曼窗');[H,W]=freqz(hr,1);

%求滤波器频率响应 subplot(2,2,3);plot(W/pi,abs(H))

%绘制滤波器频域波形 xlabel('omega/pi');ylabel('振幅');title('FIR带通滤波器幅频特性');subplot(2,2,4);

if abs(H)~=0

%判断指数是否为0

plot(W/pi,20*log10(abs(H)));

%画滤波器频谱的分贝图 end xlabel('omega/');ylabel('振幅/分贝');title('FIR带通滤波器幅频特性');grid;

%添加网格 %---ideallp（）函数（非系统自有函数）在系统安装目录的WORK子目录ideallp.m function hd = ideallp(wc,N);% 理想低通滤波器的脉冲响应子程序 % hd = 点0 到 N-1之间的理想脉冲响应 % wc = 截止频率(弧度)% N = 理想滤波器的长度

tao =(N-1)/2;

% 理想脉冲响应的对称中心位置 n = [0:(N-1)];

% 设定脉冲响应长度 m = n-tao + eps;

% 加一个小数以避免零作除数

hd = sin(wc*m)./(pi*m);

% 理想脉冲响应

布拉克曼窗1振幅0.500406080时间FIR带通滤波器幅频特性500振幅/分贝20FIR带通滤波器幅频特性1.51振幅-50-100-15000.5/10.5000.5/1

5.5用IIR滤波器对信号进行滤波

用自己设计的IIR滤波器分别对采集的信号进行滤波，在Matlab中，IIR滤波器利用函数filter对信号进行滤波。程序如下： [x,fs,bits]=wavread('女音.wav');N=length(x);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半 y=fft(x);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,1);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

%iir带通的代码：

Ap=1;

%通带波纹系数

Az=100;

%最小阻带衰减

wp=[0.2 0.6];

%归一化通带数字截止频率 wz=[0.1 0.7];

%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估计契比雪夫I型滤波器阶数

[b,a]=cheby1(N,Ap,wn);

%N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x);

%对声音滤波 wavplay(x1)wavwrite(x1,'IIR滤波后女音.wav');N=length(x1);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

y=fft(x1);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,2);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

滤波前女音的频谱特性300幅值/A***030004000频率/Hz滤波后女音的频谱特性500060006040幅值/A***0频率/Hz400050006000

结果分析：

由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3400Hz的频谱成分。回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。

5.6用FIR滤波器对信号进行滤波

用自己设计的FIR滤波器分别对采集的信号进行滤波，在Matlab中,FIR滤波器利用函数fftfilt对信号进行滤波 程序如下：

[x,fs,bits]=wavread('女音.wav');N=length(x);

%计数读取信号的点数

t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

y=fft(x);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,1);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('滤波前女音的频谱特性');

%给图形加注标签说明 xlabel('频率/omega');ylabel('幅值/A');grid;

%FIR带通滤波器代码 fs=10000;wp1=0.3*pi;wp2=0.6*pi;wz1=0.2*pi;wz2=0.7*pi;wc1=(wz1+wp1)/2;wc2=(wz2+wp2)/2;deltaW=min((wp1-wz1),(wz2-wp2));

%---取两个过渡带中的小者 N0=ceil(2*5.5*pi/deltaW);

%---查表7-3（P342）布拉克曼窗 N=N0+mod(N0+1,2);

%---确保N为奇数 hdWindow=ideallp(wc2,N)-ideallp(wc1,N);wdWindow=blackman(N);hr=wdWindow.*hdWindow';x1=fftfilt(hr,x);

%对声音滤波 wavplay(x1)wavwrite(x1,'FIR滤波后女音.wav');N=length(x1);

%计数读取信号的点数 t=(1:N)/fs;

%信号的时域采样点 f0=fs/N;

%采样间隔 n=1:N/2;

%取信号的一半

y=fft(x1);

%对信号做N点FFT变换 k=(n-1)*f0;

%频域采样点

subplot(2,1,2);

%把画图区域划分为2行1列，指定第一个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('l滤波后女音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;

滤波前女音的频谱特性300200幅值/A***004000频率/l滤波后女音的频谱特性500060006040幅值/A20005001000***03000频率/Hz***0

结果分析：

由上面滤波前后的频谱图可看出，滤波器滤除了小于1000Hz和大于3500Hz的频谱成分。和用IIR滤波器滤波一样，回放语音信号，由于低频和高频成分被滤除，声音变得较低沉。5.7男女声语音信号频谱特点分析

换一个男音录制同样一段语音内容，分析两段内容相同的语音信号频谱之间有什么特点。程序如下：

[x,fs,bits]=wavread('女音.wav');

%读取声音

N=length(x);

%计数读取信号的点数 t=(1:N)/fs;

f0=fs/N;

n=1:N/2;

y=fft(x);

k=(n-1)*f0;

subplot(2,1,1);

plot(k,abs(y(n)));

title('FFT变换后女音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid;

[x,fs,bits]=wavread('明明.wav');

N=length(x);

t=(1:N)/fs;

f0=fs/N;

n=1:N/2;

y=fft(x);

k=(n-1)*f0;

subplot(2,1,2);

plot(k,abs(y(n)));

title('FFT变换后男音的频谱特性');xlabel('频率/omega');ylabel('幅值/A');grid;

%信号的时域采样点

%采样间隔

%取信号的一半

%对信号做N点FFT变换

%频域采样点

%把画图区域划分为2行1列，指定第一个图%绘制原始语音信号的幅频响应图

%给图形加注标签说明

%添加网格

%读取声音

%计数读取信号的点数

%信号的时域采样点

%采样间隔

%取信号的一半

%对信号做N点FFT变换

%频域采样点

%把画图区域划分为2行1列，指定第二个图%绘制原始语音信号的幅频响应图

%给图形加注标签说明

%添加网格

axis([0 6000 0 300])；

%改变横纵坐标便于比较频谱图

FFT变换后女音的频谱特性300200幅值/A***00频率/FFT变换后男音的频谱特性***200幅值/A***00频率/400050006000

结果分析：

通过比较上面女音频谱图和男音频谱图可知，男音的频谱集中在低频部分，高频成分底，谱线较平滑，声音听起来低沉。5.8有背景噪声的信号分析

从硬盘中把一段噪声（频谱能量集中在某个小范围内）叠加到语音信号中，分析叠加前后信号频谱的变化，设计一个合适的滤波器，能够把该噪声滤除； 程序如下：

z=wavread('女音.wav',[1 24000]);

%读取声音在1-24000之间 f=wavread('noise.wav',[1 24000]);x=z+f;wavplay(x);fs=11025;N=length(x);f0=fs/N;

%采样间隔

n=1:N;

%取信号的一半 y=fft(x,N);%对信号做N点FFT变换

k=(n-1)*f0;

%频域采样点

subplot(2,1,1);

%把画图区域划分为1行2列，指定第二个图 plot(k,abs(y(n)));

%绘制原始语音信号的幅频响应图 title('加噪声后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格

%iir带通滤波器的代码：

Ap=1;

%通带波纹系数

Az=70;

%最小阻带衰减

wp=[0.2 0.7];

%归一化通带数字截止频率 wz=[0.1 0.8];

%归一化阻带数字截止频率 [N,wn]=cheb1ord(wp,wz,Ap,Az);

%估计契比雪夫I型滤波器阶数

[b,a]=cheby1(N,Ap,wn);

%N指定滤波器阶数，wn归一化截止频率，Ap通带波动 x1=filter(b,a,x);

%对声音滤波 wavplay(x1)；

wavwrite(x1,'滤除噪音后女音.wav');N=length(x1);f0=fs/N;

%采样间隔 n=1:N;

%取信号的一半

y1=fft(x1,N);

%对信号做fs点FFT变换

subplot(2,1,2);

%把画图区域划分为1行2列，指定第二个图 k=(n-1)*f0;

%频域采样点

plot(k,abs(y1(n)));

%绘制原始语音信号的幅频响应图 title('滤除噪声后声音的频谱特性');

%给图形加注标签说明 xlabel('频率/Hz');ylabel('幅值/A');grid;%添加网格

加噪声后声音的频谱特性3000幅值/A***0008000频率/Hz滤除噪声后声音的频谱特性***030幅值/A***000频率/Hz80001000012000

结果分析

观察加噪声后声音的频谱图可知，噪音频率主要在4000Hz处，只要我们设计一个，滤波器滤除大概在4000Hz的频谱即可，回放滤波后的语音信号，可证噪音基本滤除。

六．心得体会：

通过这次课程设计，让我对MATLAB的基本应用有了更深的了解，还有数字信号处理在MATLAB中的一些函数的用法。通过理论推导得出相应结论，并利用MATLAB作为工具进行实现，从而复习巩固课堂所学的理论知识，提高对所学知识的综合应用能力，并从实践上初步实现对数字信号的处理。

在这次实验中，也遇到了很多问题，比如画信号频谱的分贝图时（20*log10(abs(y))）指数为零时的处理。滤波器的设计也花了好大的功夫，刚开始不会设计参数，一头雾水，通过同学的指导和讨论，得知通过观察信号的频谱图，看噪音频率集中在那一部分，设计滤波器把其滤除即可。可反复设置参数直到滤波后语音信号的效果好为止。

七．参考文献：

（1）《MATLAB LabVIEW SystemView》翁剑枫 叶志前 编著, 机械工业出版社；

（2）《MATLAB及在电子信息课程中的应用》陈怀琛 吴大正 高西全编著，电子工业出版社；

（3）《MATLAB在数字信号处理中的应用》（弟2版）薛年喜 编著，清华大学出版社；

（4）《MATLAB扩展编程》何强 何英

编著，清华大学出版社；（5）《MATLAB7简明教程》吴清 曹辉林 编著，清华大学出版社；（6）MATLAB5.3精要.编程及高级应用》程卫国 冯峰 王雪梅 刘艺 编著，机械工程出版社。

第三篇：数字信号处理课程总结（推荐）

数字信号处理课程总结

信息09-1班 陈启祥 金三山 赵大鹏 刘恒

进入大三，各种专业课程的学习陆续展开，我们也在本学期进行了数字信号处理这门课程的学习。

作为信心工程专业的核心课程之一，数字信号处理的重要性是显而易见的。在近九周的学习过程中，我们学习了离散时间信号与系统的时域及频域分析、离散傅里叶变换、快速傅里叶变换、IIR及FIR数字滤波器的设计及结构等相关知识，并且在实验课上通过MATLAB进行了相关的探究与实践。总体来说，通过这一系列的学习与实践，我们对数字信号处理的有关知识和基础理论已经有了初步的认知与了解，这对于我们今后进一步的学习深造或参加实际工作都是重要的基础。

具体到这门课程的学习，应当说是有一定的难度的。课本所介绍的相关知识理论性很强，并且与差分方程、离散傅里叶级数、傅里叶变换、Z变换等数学工具联系十分紧密，所以要真正理解课本上的相关理论，除了认真聆听老师的讲解，还必须要花费大量时间仔细研读课本，并认真、独立地完成课后习题。总之，理论性强、不好理解是许多同学对数字信号处理这门课程的学习感受。

另外，必须要说MATLAB实验课程的开设是十分必要的。首先，MATLAB直观、简洁的操作界面对于我们真正理解课堂上学来的理论知识帮助很大；其次，运用MATLAB进行实践探究，也使我们真正意识到，在信息化的今天，研究数字信号离不开计算机及相关专业软件的帮助，计算机及软件技术的发展，是今日推动信息技术发展的核心动力；最后，作为信息工程专业的学生，在许多学习与实践领域需要运用MATLAB这样一个强大工具，MATLAB实验课程的开设，锻炼了我们的实践能力，也为我们今后在其他领域运用MATLAB打下了基础。

课程的结束、考试的结束不代表学习的结束，数字信号处理作为我们专业的基础之一，是不应当被我们抛之脑后的。

最后感谢老师这几周来的教诲与指导，谢谢老师！

2012年5月7日

第四篇：数字信号处理课程总结(全)

数字信号处理课程总结

以下图为线索连接本门课程的内容：

xa(t)数字信号前置滤波器A/D变换器处理器D/A变换器AF（滤去高频成分）ya(t)x(n)

一、时域分析

1． 信号

 信号：模拟信号、离散信号、数字信号（各种信号的表示及关系） 序列运算：加、减、乘、除、反褶、卷积  序列的周期性：抓定义

njwna、e(n)（可表征任何序列）cos(wn)u(n)、 典型序列：、、RN(n)、x(n)x(m)(nm)

m特殊序列：h(n)2． 系统

 系统的表示符号h(n) 系统的分类：y(n)T[x(n)]

线性：T[ax1(n)bx2(n)]aT[x1(n)]bT[x2(n)] 移不变：若y(n)T[x(n)]，则y(nm)T[x(nm)] 因果：y(n)与什么时刻的输入有关 稳定：有界输入产生有界输出

 常用系统：线性移不变因果稳定系统  判断系统的因果性、稳定性方法  线性移不变系统的表征方法：

线性卷积：y(n)x(n)*h(n)

NMk差分方程： y(n)ak1y(nk)bk0kx(nk)3． 序列信号如何得来？

xa(t)x(n)抽样

 抽样定理：让x(n)能代表xa(t) 抽样后频谱发生的变化？  如何由x(n)恢复xa(t)？

sin[xa(mT)T(tmT)]

xa(t)=mT

(tmT)

二、复频域分析（Z变换）

时域分析信号和系统都比较复杂，频域可以将差分方程变换为代数方程而使分析简化。A． 信号 1．求z变换

定义：x(n)X(z)x(n)znn

收敛域：X(z)是z的函数，z是复变量，有模和幅角。要其解析，则z不能取让X(z)无穷大的值，因此z的取值有限制，它与x(n)的种类一一对应。

 x(n)为有限长序列，则X(z)是z的多项式，所以X(z)在z=0或∞时可能会有∞，所以z的取值为：0z；

 x(n)为左边序列，0zRx，z能否取0看具体情况；

 x(n)为右边序列，Rxz，z能否取∞看具体情况（因果序列）；  x(n)为双边序列，RxzRx 2．求z反变换：已知X(z)求x(n)

 留数法

 部分分式法（常用）：记住常用序列的X(z)，注意左右序列区别。 长除法：注意左右序列 3．z变换的性质：

 由x(n)得到X(z)，则由x(nm)zmX(z)，移位性；  初值终值定理：求x(0)和x()；

 时域卷积和定理：y(n)x(n)*h(n)Y(z)X(z)H(z)；  复卷积定理：时域的乘积对应复频域的卷积；  帕塞瓦定理：能量守恒



nx(n)212X(ejw)dw2

4．序列的傅里叶变换

公式：X(ejw)x(n)enjwn

x(n)12X(ej)ejnd

注意：X(ejw)的特点：连续、周期性；X(ejw)与X(z)的关系 B． 系统

由h(n)H(z)，系统函数，可以用来表征系统。

 H(z)的求法：h(n)H(z)；H(z)=Y(z)/X(z)；  利用H(z)判断线性移不变系统的因果性和稳定性  利用差分方程列出对应的代数方程

MNMy(n)ak1y(nk)kbk0x(nk)kY(z)X(z)bk0Nkzk

k1ak1zk 系统频率响应H(ejw)：以2为周期的的连续函数



H(e)jwh(n)enjwn

H(ejw)h(n)enjwn，当h(n)为实序列时，则有H(ejw)=H*(ejw)

三、频域分析

根据时间域和频域自变量的特征，有几种不同的傅里叶变换对

 时间连续，非周期频域连续（由时域的非周期造成），非周期（由时域的连续造成）； X(j)x(t)ejtdt

x(t)12X(j)ejtd

 时间连续，周期频域离散，非周期

X(jk0)1T0T0/2x(t)ejk0tdt

T0/2x(t)X(jk0)ejk0t

 时间离散，非周期频域连续，周期



X(e)jwx(n)enjwn

x(n)12X(ej)ejnd，wT（数字频率与模拟频率的关系式）

 时间离散，周期频域离散，周期

~X(k)N1n0~x(n)ej2Nkn~x(n)W

knNn0N11~x(n)NN1n0~X(k)ej2Nkn1NN1n0~knX(k)WN

 本章重点是第四种傅里叶变换-----DFS  注意：

x(n)和X(k)都是以N为周期的周期序列； 1）~x(n)和X(k)的定义域都为（,）

2）尽管只是对有限项进行求和，但~；

~~~例如：k0时，X(0)N1x(n)

n0~~k1时，X(1)N1n0~x(n)ej2Nn

2NNnN1~kN时，X(N)N1n0j~x(n)en02N~~x(n)=X(0)

~kN1时，X(N1)N1n0~x(n)ej(N1)n~X(1)

x(n)也有类似的结果。x(n)和X(k)一

同理也可看到~可见在一个周期内，~~一对应。

 比较X(e)jwx(n)enjwn~和X(k)N1n0~x(n)ej2Nkn~x(n)W，当x(n)knNn0N1x(n)的一个周期内有定义时，即x(n)=~x(n)，0nN1，则在只在~N12Nj2Nk时，X(ejw)X(k)。

1,kr 0,kr~ en0(kr)nx(n)和X(k)的每个周期值都只是其主值区间的周期延拓，所以求和 因为~~在任一个周期内结果都一样。

 DFT：有限长序列x(n)只有有限个值，若也想用频域方法分析，它只属于序列的傅里叶变换，但序列的傅氏变换为连续函数，所以为方便计算机处理，也希望能像DFS一样，两个域都离散。将x(n)想象成一个周期x(n)的一个周期，然后做DFS，即 序列~

~X(k)N1n0~x(n)ej2NknN1n0x(n)ej2Nkn

x(n)只有x(n)，不是真正的周期序列，但因为求和只需N注意：实际上~个独立的值，所以可以用这个公式。同时，尽管x(n)只有N个值，但依上式求出的X(k)还是以N为周期的周期序列，其中也只有N个值独立，这样将~X(k)规定在一个周期内取值，成为一个有限长序列，则会引出

N1j2Nkn~DFT X(k)x(n)en0RN(k)

x(n)1NN1n0X(k)ej2NknRN(n)

比较：三种移位：线性移位、周期移位、圆周移位

三种卷积和：线性卷积、周期卷积、圆周卷积

重点：1）DFT的理论意义，在什么情况下线性卷积=圆周卷积 2）频域采样定理：掌握内容，了解恢复

3）用DFT计算模拟信号时可能出现的几个问题，各种问题怎样引起？

混叠失真、频谱泄漏、栅栏效应

 FFT：为提高计算速度的一种算法

1）常用两种方法：按时间抽取基2算法和按频率抽取基2算法，各自的原理、特点是什么，能自行推导出N小于等于8的运算流图。2）比较FFT和DFT的运算量； 3）比较DIT和DIF的区别。

四、数字滤波器（DF）

一个离散时间系统可以用h(n)、H(z)、差分方程和H(ejw)来表征。问题：

1、各种DF的结构

2、如何设计满足要求指标的DF？

3、如何实现设计的DF？

A． 设计IIR DF，借助AF来设计，然后经S---Z的变换即可得到。

1）脉冲响应不变法：思路、特点 2）双线性变换法：思路、特点、预畸变 3）模拟滤波器的幅度函数的设计 B． 设计FIR DF 1）线性相位如何得到？条件是什么？各种情况下的特点。2）窗函数设计法：步骤、特点 3）频率抽样法：步骤、特点 C． 实现DF

Ma

标准形式：H(z)k0Nkzk

bkzk1k1

第五篇：数字信号处理实验-FFT的实现

实

验

报

告

学生姓名：

学 号：

指导教师：

一、实验室名称：数字信号处理实验室

二、实验项目名称：FFT的实现

三、实验原理：

一．FFT算法思想：

1．DFT的定义：

对于有限长离散数字信号{x[n]}，0  n  N-1，其离散谱{x[k]}可以由离散付氏变换（DFT）求得。DFT的定义为：

N1X[k]通常令ej2Nx[n]en0j2Nnk，k=0,1,…N-1 WN，称为旋转因子。

2．直接计算DFT的问题及FFT的基本思想：

由DFT的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N点DFT需要（N-1）2次复数乘法和N（N-1）次加法。因此，对于一些相当大的N值（如1024）来说，直接计算它的DFT所作的计算量是很大的。

FFT的基本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。例如，若N为偶数，将原有的N

22点序列分成两个（N/2）点序列，那么计算N点DFT将只需要约[(N/2)·2]=N/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT（假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。

3．基2按时间抽取（DIT）的FFT算法思想：

设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

将长度为N2L的序列x[n](n0,1,...,N1)，先按n的奇偶分成两组：

x[2r]x1[r]x[2r1]x2[r]，r=0,1,…,N/2-1 DFT化为：

N1N/21N/21X[k]DFT{x[n]}N/21n0x[n]WnkN2rkr0x[2r]W2rkNr0x[2r1]WN(2r1)kN/21r0N/21x1[r]Wx1[r]W2rkNWWkNr0N/21x2[r]WN

r0rkN/2kNr0x2[r]WN/22rkrk上式中利用了旋转因子的可约性，即：WNN/21NrkN/21rkWN/2。又令

rkX1[k]r0x[1r]W,/X2[k]2r0x[r]WN2，则上式可以写成： /2X[k]X1[k]WNX2[k](k=0,1,…,N/2-1)

k可以看出，X1[k],X2[k]分别为从X[k]中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了X[k]前N/2点的DFT值，还需要继续利用X1[k],X2[k]表示X[k]的后半段本算法推导才完整。利用旋转因子的周期性，有：WN/2WN/2X1[N2N/21rkr(kN/2)，则后半段的DFT值表达式：

rkk]r0x1[r]W2N/2r(Nk)N/21r0x1[r]WN/2X1[k]，同样，X2[N2k]X2[k]

（k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表达式获得，中间还利用到WN(N2k)NWN2Wk得到后半段的X[k]值表达式W，k为：X[k]X1[k]WNkX2[k]（k=0,1,…,N/2-1）。

这样，通过计算两个N/2点序列x1[n],x2[n]的N/2点DFTX1[k],X2[k]，可以组合得到N点序列的DFT值X[k]，其组合过程如下图所示：

X1[k] X1[k]WNkX2[k]

X2[k] WNnk-1 X1[k]WNkX2[k]

比如，一个N = 8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示：

x(0)W0x(1)W0x(2)W0x(3)W2W0W1W0x(5)W0x(6)W0x(7)W2X(7)W3X(6)W2X(5)X(3)X(2)X(1)X(0)x(4)X(4)

4．基2按频率抽取（DIF）的FFT算法思想：

设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

在把X[k]按k的奇偶分组之前，把输入按n的顺序分成前后两半：

N1N/21nkNN1X[k]DFT{x[n]}N/21N/21x[n]Wn0(nn0N2)kx[n]WnkNnN/2x[n]WNnkn0N/21x[n]WnkNn0x[nNkN2]WNnk

Nn0[x[n]x[nN2NkN2]W2N]WN,k0,1,...,N1因为W2N1，则有WX[k](1)，所以：

kkN/21n0[x[n](1)x[nN2]]WN,k0,1,...,N1

nk按k的奇偶来讨论，k为偶数时：

N/21X[2r]n0[x[n]x[nN2]]WN,k0,1,...,N1 N22rnN/21k为奇数时：X[2r1]前面已经推导过WNN/21n0[x[n]x[n]]WN(2r1)n,k0,1,...,N1

2rkWN/2，所以上面的两个等式可以写为：

N2]]WN/2,r0,1,...,N/21 N2rnrkX[2r]n0[x[n]x[nN/21X[2r1]n0{[x[n]x[n]]WN}WN/2,r0,1,...,N/21

nnr通过上面的推导，X[k]的偶数点值X[2r]和奇数点值X[2r1]分别可以由组合而成的N/2点的序列来求得，其中偶数点值X[2r]为输入x[n]的前半段和后半段之和序列的N/2点DFT值，奇数点值X[2r1]为输入x[n]的前半段和后半段之差再与WN相乘序列的N/2点DFT值。

令x1[n]x[n]x[nN/21nN2]，x2[n][x[n]x[nN/21N2]]WN，则有：

nX[2r]n0x1[n]WrnN/2,X[2r1]n0x2[n]WrnN/2,r0,1,...,N21

这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下图所示：

x[n] x[n]x[nN2]

x[nN2]-1 WNn [x[n]x[nN2]]WNn

二．在FFT计算中使用到的MATLAB命令：

函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R

四、实验目的：

离散傅氏变换（DFT）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。

本实验通过直接计算DFT，利用FFT算法思想计算DFT，以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。

五、实验内容：

a)计算实数序列x(n)cos516n,0n256的256点DFT。

b)计算周期为1kHz的方波序列（占空比为50％，幅度取为＋/-512，采样频率为25kHz，取256点长度）256点DFT。

六、实验器材（设备、元器件）：

安装MATLAB软件的PC机一台，DSP实验演示系统一套。

七、实验步骤：

（1）先利用DFT定义式，编程直接计算2个要求序列的DFT值。

（2）利用MATLAB中提供的FFT函数，计算2个要求序列的DFT值。（3）(拓展要求)不改变序列的点数，仅改变DFT计算点数（如变为计算1024点DFT值），观察画出来的频谱与前面频谱的差别，并解释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。

（4）利用FFT的基本思想（基2－DIT或基2－DIF），自己编写FFT计算函数，并用该函数计算要求序列的DFT值。并对前面3个结果进行对比。

（5）(拓展要求)尝试对其他快速傅立叶变换算法（如Goertzel算法）进行MATLAB编程实现，并用它来计算要求的序列的DFT值。并与前面的结果进行对比。

（6）（拓展要求）在提供的DSP实验板上演示要求的2种序列的FFT算法（基2－DIT），用示波器观察实际计算出来的频谱结果，并与理论结果对比。

八、实验数据及结果分析：

程序：（1）对要求的2种序列直接进行DFT计算的程序

（2）对要求的2种序列进行基2－DIT和基2－DIF FFT算法程序（3）对要求的2种序列用MATLAB中提供的FFT函数进行计算的程序

结果：（1）对2种要求的序列直接进行DFT计算的频域波形

（2）对2种要求的序列进行基2－DIT和基2－DIF FFT算法频域波形（3）对2种要求的序列用MATLAB中提供的FFT函数计算的频域波形。（4）（拓展要求）分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。

（5）（拓展要求）保持序列点数不变，改变DFT计算点数（变为1024点），观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。

九、实验结论：

十、总结及心得体会：

十一、对本实验过程及方法、手段的改进建议：