通信1103范进文

1、天津商业大学通信工程专业 2011 级数字信号处理课程设计报告专业班级： 通信 1103 班 姓 名： 范进文 学 号： 20114879 总评成绩：指导教师签章：设计一：利用 FFT 做连续信号的频谱分析一、设计要求用 FFT计算下列连续时间信号的频谱， 并分析选择不同的抽样频率 Fs 和抽样点 数 N 对频谱特性的影响(是否有混叠现象、泄露现象、栅栏现象) 。xa (t) e 0.01t sin2 f1t sin2 f2t sin2 f3t , f1 2Hz, f2 2.1Hz, f3 2.2Hz1)Fs5Hz,N1282)Fs5Hz,N323)Fs5Hz,N10244)Fs4Hz,N10

2、24二、设计所用 MATLAB 函数( 1) fft 函数 功能：一维快速离散傅立叶变换。调用格式：y=fft(x) ：利用 fft 算法计算向量 x 的快速离散傅立叶变换。 y=fft(x,n) ：利用 fft 算法计算向量 x 的 n 点快速离散傅立叶变换。当 x 的长 度小于 n时，在 x 的尾部补零；当 x 的长度大于 n时， fft 函数会截断 x。 ( 2) plot 函数功能：绘制二维图形。调用格式：y=plot(x) ：当 x 是实向量时，以该元素的下标为横坐标，元素值为纵坐标画 出一条连续曲线，实际是绘制折线图。y=plot(x,y) ：以 x 元素为横坐标值， y 元素为纵

3、坐标值绘制曲线。y= plot(x,y1,x,y2.)：以公共的 x 元素为横坐标值，以 y1， y2 等纵坐标绘制多条曲线。( 3) axis 函数 功能：横纵坐标控制 调用格式： axis(xmin xmax ymin ymax zmin zmax)( 4) title 函数功能：加图形标题调用格式 ;title(' ')( 5) xlabel 函数功能：加 X 轴上的说明调用格式： xlabel(' ')( 6) hold on 函数hold on/off 命令控制是保持原有图形还是刷新原有图形，不带参数的 hold 命令在 两种状态之间进行切换。三、设计

4、内容程序如下：Fs=4;Ts=1/Fs;N=1024;n=0:N-1;t=n*Ts;xa=exp(-0.01*t).*(sin(2*pi*2*t)+ sin(2*pi*2.1*t)+ sin(2*pi*2.2*t);Xa=fft(xa,N); % 计算 xa 的 N 点 ifftk=0:N-1;D=Fs/N;%计算模拟频率分辨率plot(k*D,abs(Xa),'k');% 横轴转换成模拟频率axis(0,5,0,1.1*max(abs(Xa) ;% 画出折叠频率之前的随频率变化的频域 振幅title(' 幅度频谱 ');xlabel('Hz')

5、;hold on仿真波形如下 :(1) Fs 5Hz,N 12860幅度 频 谱波形分析：无混叠现象；无泄露现象；有栅栏现象（2） Fs 5Hz, N 32波形分析：无混叠现象；无泄露现象；无栅栏现象（3） Fs 5Hz, N 1024幅 度 频 谱波形分析：无混叠现象；无泄露现象；有栅栏现象（4） Fs 4Hz, N 1024波形分析：有混叠现象；无泄露现象；有栅栏现象设计二：数字滤波器的设计一、设计要求利用双线性法实现巴特沃思 IIR 数字低通、高通、带通、带阻滤波器，具体指 标如下：1）低通 IIR 数字滤波器性能指标， fp=1.2kHz ，fs=2kHz ，Rp=1dB，As=40d

6、B， Fs=8kHz；2）高通 IIR 数字滤波器性能指标， fp=3kHz ,fs=1.2k Hz，Rp=1dB，As=100dB， Fs=8kHz；3）带通 IIR 数字滤波器性能指标， fp 1=1.2kHz，fp 2=3kHz，fs 1=1kHz，fs 2=3.2kHz， Rp=1dB，As=100dB，Fs=8kHz；4）带阻 IIR 数字滤波器性能指标， fp 1=1kHz，fp 2=3.2kHz，fs 1=1.2kHz，fs 2=3kHz， Rp=1dB，As=52dB，Fs=8kHz。画出所设计滤波器的幅频响应图，观察图形，验证所设计滤波器是否满足指标 要求。二、设计所用 MA

7、TLAB 函数（1）buttord 函数功能：确定巴特沃思滤波器的阶数 N 和 3dB截至频率 c调用格式： n,wc=buttord（wp,ws,Rp,Rs, ' s'） n,wc=buttord（wp,ws,Rp,Rs）（2）buttap 函数 功能：设计巴特沃思模拟低通原型滤波器。调用格式： z,p,k=buttap（n）（3）ploy 函数 功能：求某向量指定根所对应的多项式系数向量。调用格式： a=ploy（p）（4）lp2lp ，lp2hp ，lp2bp ，lp2bs 函数功能：将模拟低通原型滤波器分别转化为实际的低通、高通、带通、带阻滤波调用格式： bt,at=l

8、p2lp(b,a,wc)bt,at=lp2hp(b,a,wc)bt,at=lp2bp(b,a,w0,bw)bt,at=lp2bs(b,a,w0,bw)(5) bilinear 函数功能：将模拟滤波器 H ( s)变成数字滤波器 H(z) 。调用格式： bz,az=bilinear(bs,as,fs)(6) butter 函数。功能：设计低通、高通、带通和带阻的数字和模拟巴特沃思滤波器调用格式：b,a=butter(n,wc)b,a=butter(n,wc,'ftype ')b,a=butter(n,wc,'s')b,a=butter(n,wc,'ftyp

9、e ', 's')(7)freqs函数功能：计算模拟滤波器的频率响应 调用格式： h=freqs(b,a,w) h,w=freqs(b,a) h,w=freqs(b,a,n) freqs(b,a)(8)freqz 函数 功能：计算数字滤波器的频率响应 调用格式： h,w=freqz(b,a,n) h,f=freqz(b,a,n,Fs) h=freqz(b,a,w) h=freqz(b,a,f,Fs) freqz(b,a)三、设计内容1、低通 IIR 数字滤波器 fp=1200;fs=2000;Fs=8000;wp=fp/Fs*2*pi;%数字滤波器的通带截止频率ws=f

10、s/Fs*2*pi;%数字滤波器的阻带截止频率Rp=1;As=40;%输入滤波器的通阻带衰减指标n,wc=buttord(wp/pi,ws/pi,Rp,As) ;%计算数字巴特沃思滤波器的阶数 n 和3dB 截至频率 wc。其中， wp和 ws 0,1 ，1 对应抽样频率 0.5Fs 或bd,ad=butter(n,wc) ;%设计 3dB截至频率为 wc 的 n 阶巴特沃思数字滤波器, wc 为一元矢量，产生数字低通滤波器H,w=freqz(bd,ad); % 得到数字滤波器的 n 点频率响应值，这 n 个点均匀 地分布在0 ， 上，并将这 n 个频点的频率记录在 w中，相应的频响值记录在

11、H 中，缺省时 n=512。dbH=20*log10(abs(H)/max(abs(H); % 化为分贝值plot(w/2/pi*Fs,dbH,'k');title(' 幅度响应 (dB)');axis(0,Fs/2,-40,5);ylabel('dB');xlabel(' 频率 (Hz)');set(gca,'XTick',0,fp,fs,Fs/2);set(gca,'YTick',-50,-20,-3,0);grid低通 IIR 数字滤波器012002000 频 率 (Hz)4000幅 度 响

12、应 (dB)2、高通 IIR 数字滤波器 fp=3000;fs=1200;Fs=8000;wp=fp/Fs*2*pi; %数字滤波器的通带截止频率ws=fs/Fs*2*pi; %数字滤波器的阻带截止频率%wp=0.4; %ws=0.25; %数字滤波器的通带截止频率数字滤波器的阻带截止频率Rp=1;As=100;% 输入滤波器的通阻带衰减指标n,wc=buttord(wp/pi,ws/pi,Rp,As) % 计算阶数 n 和截止频率b,a=butter(n,wc,'high') % 直接求数字高通滤波器系数 h,w=freqz(b,a); % 求数字系统的频率特性dbH=20*

13、log10(abs(h)/max(abs(h); % 化为分贝值 plot(w/pi,dbH);title('幅度响应 (dB)');ylabel('dB'); xlabel('频率 (pi)');axis(0,1,-40,5);set(gca,'XTick',0,ws/pi,wp/pi,1);set(gca,'YTick',-40,-20,-1,0);grid高通 IIR 数字滤波器幅度 响 应(dB)dB00.30.75频率()3、带通 IIR 数字滤波器fp1=1200;fs1=1000; wp1=fp1/F

14、s*2*pi;%ws1=fs1/Fs*2*pi;%fp2=3000;fs2=3200;wp2=fp2/Fs*2*pi;%ws2=fs2/Fs*2*pi;%数字滤波器的通带截止频率数字滤波器的阻带截止频率数字滤波器的通带截止频率数字滤波器的阻带截止频率%ws1=0.15;ws2=0.85; % 数字滤波器的阻带截止频率 ws=ws1,ws2;%wp1=0.25;wp2=0.75; % 数字滤波器的通带截止频率 wp=wp1,wp2;Rp=1;As=100; % 输入滤波器的通阻带衰减指标n,wc=buttord(wp/pi,ws/pi,Rp,As) % 计算阶数 n 和截止频率 b,a=butt

15、er(n,wc) % 直接求数字带通滤波器系数H,w=freqz(b,a); % 求数字系统的频率特性 dbH=20*log10(abs(H)/max(abs(H); % 化为分贝值plot(w/pi,dbH,'k');幅度响应 (dB)') ;axis(0,1,-40,5);ylabel('dB)');title(' set(gca,'XTick',0,ws1/pi,wp1/pi,wp2/pi,ws2/pi,1);set(gca,'YTick', -40,-20,-1,0);grid带通 IIR 数字滤波器00.

16、250.30.750.84、带阻 IIR 数字滤波器 fp1=1000;fs1=1200; Fs=8000; wp1=fp1/Fs*2*pi;%ws1=fs1/Fs*2*pi;%fp2=3200;fs2=3000; wp2=fp2/Fs*2*pi;%ws2=fs2/Fs*2*pi;%数字滤波器的通带截止频率数字滤波器的阻带截止频率数字滤波器的通带截止频率数字滤波器的阻带截止频率%wp1=0.15;wp2=0.85; % 数字滤波器的通阻带截止频率 wp=wp1,wp2;%ws1=0.25;ws2=0.75; % 数字滤波器的阻带截止频率ws=ws1,ws2;Rp=1;As=52; % 输入滤波

17、器的通阻带衰减指标 n,wc=buttord(wp/pi,ws/pi,Rp,As) % 计算阶数 n 和 3dB截止频率 b,a=butter(n,wc,'stop') % 直接求数字带通滤波器系数 H,w=freqz(b,a); % 求数字系统的频率特性 dbH=20*log10(abs(H)/max(abs(H); % 化为分贝值 plot(w/pi,dbH,'k');ylabel('dB)');title(' 幅度响应 (dB)') ;axis(0,1,-40,5);set(gca,'XTick',0,wp1

18、/pi,ws1/pi,ws2/pi,wp2/pi,1);set(gca,'YTick', -40,-20,-1,0);grid带阻 IIR 数字滤波器幅 度 响 应 (dB)01设计三：语音信号的消噪一、设计要求利用 windows 下的录音机录制一段自己的语音并储存，时间控制在 3 秒左右； 对语音信号叠加正弦噪声 0.2sin(2 ft)， f 2500 ；分析原语音信号和加噪语音信号 的频谱；设计一个巴特沃斯型 IIR 数字滤波器对正弦噪声进行滤除，并比较滤波前 后语音信号的频谱的变化；最后通过 MATLAB中的函数 sound 对声音进行回放，感 觉滤波前后声音的变化。

19、二、设计所用 MATLAB 函数(1) wavread 函数功能：把 .wav 的语音赋给数据变量。调用格式： x,fs,bits=wavread('filename') ，将名为 filename 的.wav 的语音转换 成数组格式的数据 x，同时把 x 的采样频率 fs 和数据位数 bits 放进 MATLAB 的工 作空间。(2) wavwrite函数功能：将声音数据存成 .wav格式。调用格式： wavwrite(x,fs,bits,'filename') ，将数组格式的数据 x按指定采样频率 fs 和数据位数 bits存入微软的 .wav文件中，并命名

20、为 filename。(3) sound函数功能：对声音进行回放。调用格式： sound(x,fs,bits)，将数据 x 通过声卡转化为声音放出。(4) filter 函数 功能：用数字滤波器对输入信号进行滤波。调用格式： y=filter(b,a,x) ，对信号 x 进行滤波，结果存放于 y 中，b、a 分别为 是数字滤波器传递函数分子、分母多项式的系数。(5) ploy 函数功能：求某向量指定根所对应的多项式系数向量。调用格式： a=ploy(p)(6) axis 函数功能：横纵坐标控制调用格式： axis（xmin xmax ymin ymax zmin zmax）（7）title 函

21、数功能：加图形标题 调用格式 ;title（' '）（8）hold on 函数hold on/off 命令控制是保持原有图形还是刷新原有图形，不带参数的 hold 命 令在两种状态之间进行切换。三、设计内容设计步骤：（1）语音信号的采集（2）语音信号的频谱分析（给出 0-4000Hz 之间的）（3）设计数字滤波器和画出频率响应（4）用滤波器对信号进行滤波（5）比较滤波前后语音信号的波形及频谱（6）回放语音信号1. 音频采集在 windows 操作系统下点击 开始程序 附件 娱乐录音机 ，将 出现图 2 所示的录音机面板。windows 操作系统下的录音机面板2. 程序设计 x,

22、Fs,bits=wavread('E:liuchao.wav'); N=length(x);Ts=1/Fs;n=0:N-1;Xa=fft(x,N); k=0:N-1;D=Fs/N;%计算模拟频率分辨率 subplot(2,3,1);plot(x);title(' 原信号幅度谱 ');%原语音信号 subplot(2,3,2);stem(n,x,'k');plot(k*D,abs(Xa),'k');% 横轴转换成模拟频率 axis(0,4000,0,0.2*max(abs(Xa) ;%画出折叠频率之前的随频率变化的频域振幅 titl

23、e(' 原信号频谱 ');xlabel('Hz');hold on%(转换模拟的语音信号 ) y=0.2*sin(2*pi*2500*n*Ts);% 加噪 z=y+x'wavwrite(z,Fs,bits,'e:1.wav'); Nx=length(z);Ts=1/Fs;a=0:N-1;Xa=fft(z,Nx); k=0:N-1;D=Fs/N;%计算模拟频率分辨率 subplot(2,3,3);stem(a,x,'k');plot(k*D,abs(Xa),'k');% 横轴转换成模拟频率 axis(0,40

24、00,0,0.1*max(abs(Xa) ;%画出折叠频率之前的随频率变化的频域振幅 title(' 加噪信号频谱 ');xlabel('Hz');hold on% ( 叠加正弦信号后的频谱 )%直接法设计数字带阻滤波器 fp1=1800;fs1=2000;fp2=2800;fs2=2600;Fs=22050;wp1=fp1/Fs*2;%数字滤波器的通阻带截止频率wp2=fp2/Fs*2;wp=wp1,wp2; ws1=fs1/Fs*2;%数字滤波器的阻带截止频率ws2=fs2/Fs*2;ws=ws1,ws2;Rp=1;As=36;%输入滤波器的通阻带衰减指标n

25、,wc=buttord(wp,ws,Rp,As) % 计算阶数 n 和 3dB 截止频率 b,a=butter(n,wc,'stop') %直接求数字带通滤波器系数 H,w=freqz(b,a);%求数字系统的频率特性dbH=20*log10(abs(H)/max(abs(H); % 化为分贝值 subplot(2,3,4) plot(w/pi,dbH,'k');ylabel('dB)');title(' 带阻频谱响应 (dB)') ;axis(0,1,-40,5); set(gca,'XTick',0,wp1,w

26、s1,ws2,wp2,5);set(gca,'YTick', -50,-20,-1,5);grid% y=2*sin(2*pi*2500*n*Ts);z=y+x'y=filter(b,a,x)wavwrite(y,Fs,bits,'e:2.wav');N=length(y);a=0:N-1;Xa=fft(y,N);k=0:N-1;D=Fs/N;%计算模拟频率分辨率 subplot(2,3,5);stem(a,x,'k');plot(k*D,abs(Xa),'g');% 横轴转换成模拟频率axis(0,4000,0,0.1*

27、max(abs(Xa) ;%画出折叠频率之前的随频率变化的频域振幅 title(' 滤波后频谱 ');xlabel('Hz');hold on sound(x,Fs,bits);3. 音频结果400200滤波后频谱00加噪信号频谱2000Hz400010203002000 4000Hz如图可见，在滤波前，原信号的频谱分布较为密集， 在 2500HZ频率附近比较平坦， 而再加入噪声以后，整体频谱起伏较大，在 2500HZ 频率附近有明显的噪声信号， 在通过了带阻滤波器以后，可以看到，在 0-2200HZ 的频段内，波形基本没有改变， 但是从 2200HZ-2800HZ 之间内的频谱图明显趋于平坦， 而原来的在 2500HZ 附近的 尖锐噪声信号已经消失。滤波后的图像和原图像基本相似。设计总结为期一周的课程设计结束了，这次课程设计虽然遇到很多问题 , 很多困难，但也学到了很多东西， 不仅学到了书本上的东西， 而且学到了很多书本上没有的东西，很多程序了的东西，特别是程序语法，总是有错误，但是总是不知道错在哪里，在细心的检查， 和小组成员帮助之下， 终于找出了错误和警告， 排除困难后，程序编译就通过了，