1、数字信号处理报告数字信号处理课程设计报告设计一 信号的基本运算【一】 设计目的熟悉信号的基本运算，通过仿真结果加深理解。【二】 设计内容1.现有一段以采样频率8192Hz采样得到的语音数据x(k)，为了模仿出回音的效果，可以对数据进行如下处理：x(k)=x(k)+a*x(k-d)，其中d为时延，a为时延信号的衰减幅度。根据上述原理，进行仿真实验，并演示回声的效果。(参数：时延d=0.4秒，衰减幅度a=0.5)2.对上述语音信号进行时间反转x(-k)、上采样x(2k)、下采样x(k/2)操作，并演示这些运算的效果。【三】 设计过程，调试，分析。数字信号处理的主要研究对象是语音信号和图像信号,语音
2、信号的研究可以从时域和频域两个方面来进行。其中时域的分析处理有两种方法: 一种是进行语音信号分析, 这属于线性处理的范畴,主要是通过信号的加减、时移、倍乘、卷积、求相关函数等来实现; 另一种是生成和变换成各种调制信号, 这属于非线性的范畴,主要是对信号平均累加器的动 态范围进行压缩扩张, 用门限方法对噪声的抑制。对频域分析处理, 即对信号的频率特性 在频谱中加以分析研究, 这拓展了信号分析的范围, 是对不确定信号分析的主要方法。在实际应用中, 信号的时域频分析经常同时进行。1 语音信号的采集与时频分析1. 1 语音信号的采集利用 PC 机上的声卡和WINDOWS操作系统可以进行数字信号的采集。
3、将话筒输入计算机的语音输入插口上, 启动录音机。按下录音按钮, 接着对话筒说话/大家好0 ,说完后停止录音, 屏幕左侧将显示所录声音的长度。点击放音按钮,可以实现所录音的重现。以 文件名/ dajiahao0 保存入 g:MATLAB work中。可以看到, 文件存储器的后缀默认wav。在 Matlab 软件平台下可以利用函数wavread对语音信号进行采样,得到了声音数据变 量x,同时把x的采样频率Fs= 8000Hz和数据位Nbits= 8Bit放进了MATALB的工作空间。其程序如下Fs= 8000Hz; % 给出抽样频率Nbits= 8Bit; % 放入数据位 x, fs, Nbits
4、 = wavread(- g: Matlab wor k chaotian.wav ) ; % 把语音信号进行加载入 Matlab 仿真软件平台中sound( x, fs) ; % 对加载的语音信号进行回放stem( x) ; title( 语音信号的时域波形 ) ;1. 2 语音信号的时域分析因录音采样频率无法准确。故用余弦序列表示信号的迟延，反转，上采样下采样。迟延：题目要求时延d为0.4s.取点数k为0.4*Fs。产生一组新序列将原序列依次推后k点。【四】 程序清单与运行结果程序清单：y=load(E:erf1s1t0);fs=8192;n=length(y);m=n+round(0.4
5、*fs);q=(m-1)/fs; l=round(0.4*fs);z=zeros(1,m);x=z;for i=1:n x(i)=y(i);end %原信号for i=1:n z(l+i)=y(i);end %信号的延拖o=x+0.4*z; %回声信号sound(o) %声音t=0:1/fs:q ;subplot(3,1,1);plot(t,x)title(原信号);subplot(3,1,2);plot(t,z)title(时延信号);subplot(3,1,3);plot(t,o)title(回声信号);运行结果：翻转：K为序列长度，将序列第i个点值依次与第k-i点值进行调换产生新序列，即
6、为翻转序列。采样：下采样，也就是对信号的抽取。其实，上采样和下采样都是对数字信号进行重采，重采的采样率与原来获得该数字信号的采样率比较，大于原信号的称为上采样，小于的则称为下采样。上采样的实质也就是内插或插值。程序清单：%对上述语音信号进行时间反转x(-k)、上采样x(2k)、下采样x(k/2)操作data=load(E:erf1s1t0); x=data; %原声信号d=0.4; %d为时延Fs=8192; %Fs为采样频率k=0.4*8192; a=0.5; %a为时延信号的衰减幅度z=linspace(0,0,k); %linspace为线性间隔的向量y=z x;x=x z;f=x+a*
7、y; %回声信号f1=fliplr(x); %翻转信号f2=downsample(x,2); %下采样信号f3=upsample(x,2); %上采样信号subplot(4,1,1);plot(x);title(原声信号)subplot(4,1,2);plot(f1);title(翻转信号)subplot(4,1,3);plot(f2);title(下采样)subplot(4,1,4);plot(f3);title(上采样)sound(f1)运行结果：设计二 正余弦信号的谱分析【一】设计目的用DFT实现对正余弦信号的谱分析；观察DFT长度和窗函数长度对频谱的影响；对DFT进行谱分析中的误差现象
8、获得感性认识。【二】设计内容1.对一个频率为10Hz，采样频率为64Hz的32点余弦序列进行谱分析，画出其频谱图；若将频率改为11Hz，其他参数不变，重新画出该序列的频谱图，观察频率泄漏现象，分析原因。2.考察DFT的长度对双频率信号频谱分析的影响。设待分析的信号为 （1.2）令两个长度为16的正余弦序列的数字频率为及。取N为四个不同值16，32，64，128。画出四个DFT幅频图，分析DFT长度对频谱分辨率的影响。3.在上题中若把两个正弦波的频率取得较近，令，试问怎样选择FFT参数才能在频谱分析中分辨出这两个分量？【三】设计过程，调试，分析。连续信号的频谱分析在工程上有着广泛的应用，计算机的
9、第一步就是把连续信号离散化，一是采样，二是截断。由此会产生频率混叠误差、栅栏效应和截断误差。（一） 对一个频率为10Hz，和11Hz，采样频率为64Hz的32点余弦序列用Matlab计算它的DFT程序为N=32;Fs=64;f=10;n=0:N-1;x=cos(2*pi*10*n/64)y=cos(2*pi*11*n/64)subplot(2,1,1)X=abs(fft(x,N)stem(n,X,.);xlabel(n);title(10HZ 32点变化幅频曲线)subplot(2,1,2)Y=abs(fft(y,N)stem(n,Y,.);xlabel(n);title(11HZ 32点变化
10、幅频曲线)运行结果：通过图可看出：频率为10Hz的余弦曲线DFT只有两个点不等于零，位于k=5和k=27处，k=5对应于频率10Hz，k=27对应于频率54Hz(也就是-10Hz)。这样DFT确实正确的分辨了余弦信号的频率。但是这样的理想结果是恰巧得到的，此时我们借去了五个完整的余弦周期（f*N/Fs=5）.将频率改为11Hz，采样频率和窗长度依然为32点，计算图像可看出：频谱图上k=5和k=27处都有较大的峰值，而其它的点上幅度不再为零。这两个峰值对应的频率为10Hz和12Hz，所以，信号的峰值位于两者之间，本来是单一的11Hz频率的能量会分不到许多DFT频率上，这种现象叫频率泄露，来源于截
11、断效应。（二） 对双信号的频谱分析程序为：16点，32点，64点，128点的频谱分析程序：N=16;n=0:15;f1=0.22;f2=0.34;x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n);subplot(4,2,1),stem(n,x,.);%xlable(n); ylable(x1(n);title(余弦序列);X=abs(fft(x,N);subplot(4,2,2)k=0:N-1;stem(k,X,.)%xlable(k);ylable(X(k)string=num2str(N),点FFT幅频曲线;title(string);N=32;n=0:15;f1=0
12、.22;f2=0.34;x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n);subplot(4,2,3),stem(n,x,.)%xlable(n); ylable(x1(n);title(余弦序列);X=abs(fft(x,N);subplot(4,2,4)k=0:N-1;stem(k,X,.)%xlable(k);ylable(X(k)string=num2str(N),点FFT幅频曲线;title(string);N=64;n=0:15;f1=0.22;f2=0.34;x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n);subplot(4,2,5)
13、,stem(n,x,.)%xlable(n); ylable(x1(n);title(余弦序列);X=abs(fft(x,N);subplot(4,2,6)k=0:N-1;stem(k,X,.)%xlable(k);ylable(X(k)string=num2str(64),点FFT幅频曲线;title(string);N=128;n=0:15;f1=0.22;f2=0.34;x=0.5*sin(2*pi*f1*n)+sin(2*pi*f2*n);subplot(4,2,7),stem(n,x,.)%xlable(n); ylable(x1(n);title(余弦序列);X=abs(fft(x
14、,N);subplot(4,2,8)k=0:N-1;stem(k,X,.)%xlable(k);ylable(X(k)string=num2str(128),点FFT幅频曲线;title(string);运行结果:分析：DFT样本值就是其DTFT在相应位置的采样。 在图一中很难看出两个峰值，因此要提高它的分辨率，故把N增大，逐渐可以看出它有两个峰值，将k换算成数字频率f=w/2*pi=k/N.这样可确定峰值的位置大体在f=0.21和0.35之附近，与信号的给定频率有一定的误差，这也是截断和泄露带来的问题，在这图上还可以看到一些较小的峰，这是很难判断是输入信号固有的，还是由泄露引起的。说明了增加
15、DFT长度N减小了相邻样本间的频率间距，提高频谱的视在分辨率，因而可以提高样本未知的测定精度。（三） 把两个正弦波的频率取得较近，另fr=0.22,0.25应怎样选择FFT参数要能分清两个频率，分辨率至少应达到f=0.03.因为此处的数字频率是对采样频率Fs进行归一化后的，几fr最大为1，因此总的样本数至少要达到1/0.03=33加窗以后可以使频谱函数更加光滑，便于分辨峰值位置和准确的数值，为了提高实际分辨率，应该尽量增加信号的长度n及DF长度N,当受到条件限制不能提高n，则单独提高N可以提高视在分辨率。当，时，16点，32点，64点，128点的频谱运行结果：（四） 谱分析中的误差现象混叠误差
16、。产生混叠误差的原因是：信号的离散化是通过抽样实现的，而抽样频率再高总是有限的。除带限信号外，如果信号的最高频率趋于无穷，则实际器件无法满足抽样定理。而抽样过程如果不满足抽样定理，就会产生频谱的混叠，即混叠误差。要减少或避免混叠误差，应该提高抽样频率，以设法满足抽样定理，或者采用抗混叠滤波这样的信号预处理措施。栅栏效应。对于非周期信号来说，理论上应具有连续的频谱，但数字谱分析是用DFT 来近似的，是用频谱的抽样值逼近连续频谱值，只能观察到有限个频谱值，每一个间隔中的频谱都观察不到。如同通过“栅栏”观察景物一样，一部分被“栅栏”所阻挡，看不见，把这种现象称为栅栏效应。连续时间信号只要采用数字谱分
17、析的方法，就必定产生栅栏效应，栅栏效应只能减小而无法避免。能够感受的频谱最小间隔值，称为频谱分辨率，而频率分辨率与信号截断长度成反比，且频率分辨率越小，插值效果越好。因此即使连续信号是有限时间的，也应选择一个大的截断长度来改善栅栏效应。截断误差。截断误差是由于对信号进行截断，把无限长的信号限定为有限长，即令有限区间以外的函数值均为零值的近似处理而产生的，这种处理相当于用一个矩形(窗)信号乘待分析的连续时间信号。必然会引起吉布斯效应（波动），也会把窗函数的频谱引入信号频谱 ,造成混叠，所以需要考虑其误差的问题。设计三 数字滤波器的设计及实现【一】 设计目的1. 熟悉IIR数字滤波器和FIR数字滤
18、波器的设计原理和方法；2. 学会调用MATLAB信号处理工具箱中的滤波器设计函数设计各种IIR和FIR数字滤波器，学会根据滤波要求确定滤波器指标参数；3. 掌握用IIR和FIR数字滤波器的MATLAB实现方法，并能绘制滤波器的幅频特性、相频特性；4. 通过观察滤波器的输入、输出信号的时域波形及其频谱，建立数字滤波的概念。【二】 设计内容1. 调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，观察st的时域波形和幅频特性曲线；2. 通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻
19、带截止频率，假定要求滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB；3. 编程调用MATLAB滤波器设计函数分别设计这三个数字滤波器，并绘图显示其幅频特性曲线；4. 用所设计的三个滤波器分别对复合信号st进行滤波，分离出st中的三路不同载波频率的调幅信号，并绘图显示滤波后信号的时域波形和频谱，观察分离效果。【三】 设计过程，调试，分析。（一）IIR数字滤波设计方法数字滤波是数字信号处理的一种重要算法，广泛用于对信号的过滤、检测与参数的估计等信号处理中。数字滤波器按其单位脉冲响应的长度可分为有限脉冲响应（FIR）滤波器和无限脉冲响应（IIR）滤波器两类。数字滤波器按其通频带分为低通滤波
20、器，高通滤波器，带通滤波器和带阻滤波器。在本次实验中旨在用椭圆滤波器实现IIR低通，带通，高通滤波。IIR 滤波器基本概念一个 N 阶递归型数字滤波器的差分方程H ( Z) =ME br z- rr = 0N1+K =E akz- k从以上的系统函数可知, 设计 IIR 滤波器的任务就是通过计算寻求一个因果、物理上可实现的系统函数 H ( Z) , 使其频率响应 H ( ejw) 满足所希望得到的频域指标, 即符合给定的通带截止频率、阻带截止、通带衰减和阻带衰减。数字滤波器的设计过程数字滤波器的一般设计过程可分为以下四个步骤:1) 按照实际的需要, 确定滤波器的性能要求或指标; 2) 用一个因
21、果稳定的离散线性时不变系统,去逼近这一性能指标; 3) 用有限精度的运算实现所设计的系统; 4) 通过模拟, 验证所设计的系统是否符合性能要求（二）1：调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st,画出st的时域波形和幅频特性，三路信号在时域混叠无法在时域分离，但在频域分离，所以可以通过滤波的方法在频域分离。抑制载波单频调幅信号的数学表达式为 这三路调幅信号的载波频率分别为250Hz、500Hz和1000Hz2：通过观察st的幅频特性，分别可以确定三路信号的低通，带通，高通滤波器。的通带截止频率和阻带截止频率，要求滤波器的通带最大衰减为0.1db,阻带最小衰减为60
22、db 分离混合信号st的滤波器指标参数选取如下：载波为250Hz的调幅信号，用低通滤波器分离。指标为：通带截止频率fp=300HZ,阻带截止频率fs=400Hz,通带最大衰减ap=0.1db,阻带最小衰减as=60db对载波频率为500Hz的调幅信号，用带通滤波器分离，指标为：通带截止频率fl=400HZ,fu=500Hz,阻带截止频率fl=600Hz,fu=700Hz,通带最大衰减ap=0.1db,阻带最小衰减as=60db对载波频率为1000Hz的调幅信号，用高通滤波器分离，其指标为：通带截止频率fp=900HZ,阻带截止频率fs=800Hz,通带最大衰减ap=0.1db,阻带最小衰减，6
23、0db3：用椭圆滤波器实现椭圆滤波器特点：幅值响应在通带和阻带内都是等波纹的。对于给定的阶数和给定的波纹要求，椭圆滤波器能获得较其它滤波器更窄的过渡带宽，就这点而言，椭圆滤波器是最优的。通带和阻带内波纹固定时，阶数越高，过渡带越窄阶数固定，通带和阻带波纹越小，过渡带越宽【四】 程序清单与运行结果程序清单：%1. 调用信号产生函数mstg 产生由三路抑制载波调幅信号相加构成的复合信号st，观察%st 的时域% 波形和幅频特性曲线；function st=mstg%产生信号序列st，并显示st的时域波形和频谱%st=mstg返回三路调幅信号相加形成的混合信号，长度N=800N=800; %信号长度
24、N=800Fs=10000; %采样频率Fs=10kHz,Tp为采样时间T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=Fs/10; %第一路调幅信号载波频率fc1=1000Hzfm1=fc1/10; %第一路调幅信号的调制信号频率fm1=100Hzfc2=Fs/20; %第二路调幅信号载波频率fc2=500Hzfm2=fc2/10; %第二路调幅信号的调制信号频率fm2=50Hzfc3=Fs/40; %第三路调幅信号载波频率fc3=250Hzfm3=fc3/10; %第三路调幅信号的调制信号频率fm3=25Hzxt1=cos(2*pi*fm1*t
25、).*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的频谱%以下为绘图命令figure(1)subplot(2,1,1)plot(t,st);grid;xlabel(t/s);ylabel(s(t);axis(0,Tp,min(st),max(st);title(a)s(t)的波形)subplo
26、t(2,1,2)stem(f,abs(fxt)/max(abs(fxt),.);grid;title(b)s(t)的频谱)axis(0,Fs/8,0,1.2);xlabel(f/Hz);ylabel(幅度)%低通滤波器fp=320;fs=350;wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）N,wp=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wpB,A=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B和Ay1t=filter(B,A
27、,st); %滤波器软件实现y1=fft(y1t);% 绘图部分figure(2);subplot(3,3,1);H,W=freqz(B,A,800);plot(W*Fs/2/pi,abs(H);xlabel(Hz);ylabel(H(w);title(低通滤波器);axis(0,2000,0,1.2);grid;subplot(3,3,2);plot(t,y1t);xlabel(t);ylabel(y(t);title(分离出的250Hz的波形);subplot(3,3,3);stem(f,abs(y1)/max(abs(y1),.);xlabel(Hz);ylabel(|H|);title
28、(频谱);axis(0,1200,0,1);grid;%带通滤波器fpl=420;fpu=580;fsl=300;fsu=700;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); %调用ellipord计算椭圆DF阶数N和通带截止频率wpB,A=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B和Ay2t=filter(B,A,st); %滤波器软件实现y2=fft(y2t);%绘图部分figure(2);subplot(3,3,4);H,
29、W=freqz(B,A,800);plot(W*Fs/2/pi,abs(H);xlabel(Hz);ylabel(H(w);title(带通滤波器);axis(0,2000,0,1.2);grid;subplot(3,3,5);plot(t,y2t);xlabel(t);ylabel(y(t);title(分离出的500Hz的波形);subplot(3,3,6);stem(f,abs(y2)/max(abs(y2),.);xlabel(Hz);ylabel(|H|);title(频谱);axis(0,1200,0,1);grid;%高通滤波器fp=870;fs=780;wp=2*fp/Fs;w
30、s=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）N,wp=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wpB,A=ellip(N,rp,rs,wp,high); %调用ellip计算椭圆带通DF系统函数系数向量B和Ay3t=filter(B,A,st); %滤波器软件实现y3=fft(y3t);%绘图部分figure(2);subplot(3,3,7);H,W=freqz(B,A,800);plot(W*Fs/2/pi,abs(H);xlabel(Hz);ylabel(H(w);title(高通滤波器);axis(0,2000,0,1.2);grid;
