数字信号实验报告

中国地质大学（武汉）数字信号处理实验报告1实验一：快速傅立叶变换的谱分析一、实验目的：学会利用matlab中的FFT函数，即进行信号的谱分析。二、实验题目：1．已知：t1=[0:0.001:0.3];t2=[0.301:0.001:0.6];t3=[0.601:0.001:0.9];t4=[0.901:0.001:1.199];x1=sin(2*pi*100*t1);x2=sin(2*pi*50*t2);x3=sin(2*pi*25*t3);x4=sin(2*pi*10*t4);信号s1为t=0:0.001:1.199;s1=sin(2*pi*100*t)+sin(2*pi*50*t)+sin(2*pi*25*t)+sin(2*pi*10*t);信号s2为t=[t1,t2,t3,t4];s2=[x1,x2,x3,x4];信号s3为t=[t1,t2,t3,t4];s3=[x1,x4,x2,x3];a.编写程序分别画出信号s1,s2,s3的时域波形和幅频图（参考图1，2，3）。b.观察信号s1,s2,s3的时域波形和频谱图，思考其幅频图的差别及其原因。c.信号s1,s2,s3的抽样频率fs为多少？对于1200个点的时域离散序列，对其FFT后仍为长度为1200个点的序列，即周期N=1200，试分析N与离散时间信号频谱的周期fs的对应关系。各频谱图的频率分辨率为多少Hz？d.若有信号s4的幅频图与s2的幅频图完全相同（如图2所示），问s4的时域波形和s2是相同的吗，为什么?实验程序清单：t1=[0:0.001:0.3];t2=[0.301:0.001:0.6];t3=[0.601:0.001:0.9];t4=[0.901:0.001:1.199];x1=sin(2*pi*100*t1);x2=sin(2*pi*50*t2);x3=sin(2*pi*25*t3);x4=sin(2*pi*10*t4);%信号s1的时域图和频域图t5=0:0.001:1.199;x5=sin(2*pi*100*t5)+sin(2*pi*50*t5)+sin(2*pi*25*t5)+sin(2*pi*10*t5);y5=abs(fft(x5));中国地质大学（武汉）数字信号处理实验报告2f=1000*(0:(1/1200):0.5);figure(1)subplot(2,1,1),plot(t5,x5);gridon;axistight;%波形图title('时域波形');xlabel('t(s)');subplot(2,1,2),plot(f,y5(1:601));gridon;axistight;%幅频图title('幅频图');xlabel('f(Hz)');%信号s2的时域图和频谱图t6=[t1,t2,t3,t4];x6=[x1,x2,x3,x4];y6=abs(fft(x6));figure(2)subplot(2,1,1),plot(t6,x6);gridon;axistight;%波形图title('时域波形');xlabel('t(s)');subplot(2,1,2),plot(f,y6(1:601));gridon;axistight;%幅频图title('幅频图');xlabel('f(Hz)');%信号s3的时域图和频域图x7=[x1,x4,x2,x3];y7=abs(fft(x7));figure(3)subplot(2,1,1),plot(t6,x7);gridon;axistight;%波形图title('时域波形');xlabel('t(s)');subplot(2,1,2),plot(f,y7(1:601));gridon;axistight;%幅频图title('幅频图');xlabel('f(Hz)');结果分析a.编写程序分别画出信号s1,s2,s3的时域波形和幅频图（参考下图）。中国地质大学（武汉）数字信号处理实验报告3b.观察信号s1,s2,s3的时域波形和频谱图，思考其幅频图的差别及其原因。因为信号s1、s2、s3这三个信号中都有频率为100hz，50hz，25hz，10hz的中国地质大学（武汉）数字信号处理实验报告4频率分量，从波形图可以看出s1只是单纯的含有这四种分量，而s2，s3的频谱图相同，除了这四种分量还有其他分量，因为s1中这四种分量一直存在，但s2，s3中只是特定时间段出现。c.信号s1,s2,s3的抽样频率fs为多少？对于1200个点的时域离散序列，对其FFT后仍为长度为1200个点的序列，即周期N=1200，试分析N与离散时间信号频谱的周期fs的对应关系。各频谱图的频率分辨率为多少Hz？频谱图横坐标所以两相邻点的频率间隔为1000/1200=5/6hz,所以各频谱图的频率分辨率F=5/6HZ.又因为F=1/tp=1/NT=fs/N,因为N=1200,所以s1，s2，s3的抽样频率fs=NF=1200*5/6=1000HZ.d.若有信号s4的幅频图与s2的幅频图完全相同（如图2所示），问s4的时域波形和s2是相同的吗，为什么?若有信号s4的幅频图与s2的幅频图完全相同，s4的时域波形和s2可能是不同的，因为幅频特性图只表明了幅度和频率的关系，即对应频率分量的幅值，而其相位特性未知，时域图是幅度和相位图完全一致时才相同。实验小结虽然第一次实验看起来比较简单，但是由于平时不怎么接触Matlab，导致很多对Matlab语法以及命令操作不是很熟悉，但是Matlab本省提供的帮助文档很全，所以遇到不会的命令或者函数可以直接调用帮助文档，可以很快的解决问题。其次就是通过实验，反复研究代码，可以加深对快速傅立叶变换的谱分析的理解。实验二：IIR数字滤波器设计一、实验目的：1、熟悉滤波器的基本概念；2、了解滤波器的分类；3、熟悉应用matlab设计各种滤波器的方法。二、实验相关原理：1.数字滤波器分类实现方法有限长冲激响应----FIR滤波器无限长冲激响应----IIR滤波器功能低通（LP）高通（HP）带通（BP）带阻（BS）2、数字滤波器设计步骤一、给出所需要的滤波器的技术指标二、设计一个H(Z)使其逼近所需要的技术指标三、实现所设计的H(Z)（H(Z)为系统的数学模型：传递函数模型、状态方程模型和零极点增益模型中国地质大学（武汉）数字信号处理实验报告5等）3、IIR数字滤波器设计步骤一、按一定规则将给出的数字滤波器的技术指标转换为模拟滤波器的技术指标二、根据转换后的技术指标使用滤波器阶数选择函数，确定最小阶数N和固有频率Wn三、运用最小阶数N、固有频率Wn产生模拟滤波器原型四、运用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器函数1、Buttord(Wp,Ws,Rp,Rs)butterworth滤波器阶数选择函数,返回符合要求性质的滤波器最小阶数N以及固有频率WnWp通带截至频率Ws阻带截至频率Rp通带衰减（不超过）Rs阻带衰减（不小于）Wp,Ws是归一化频率，范围是[0，1]，对应д弧度归一化处理：Wp（或Ws）/(fs/2)例：确定数字低通Butterworth滤波器的阶数和固有频率。要求：wp=500hz,ws=550hz,rp=1db,rs=50db,fs=2000hz程序清单：Wp=500;ws=550;rp=1;rs=50;fs=2000;[N,Wn]=buttord(wp/(fs/2),ws/(fs/2),rp,rs)2、Butter(N,Wn)设计数字滤波器函数格式：[B,A]=Butter(N,Wn)[B,A]是滤波器系数矩阵。如果Wn=[w1,w2]是一个二元向量,则该函数将设计出带通、带阻数字滤波器3、freqz(B，A)求线性离散时间系统的传输函数的频率响应格式:[H,W]=freqz(B，A)H:频率响应W:所要计算响应的频率点组成的向量4、filter()一维数字滤波函数格式:y=filter(B，A,x)x:待滤波信号y:滤波以后的信号x1=wavread('pb8k.wav');%读取语音信号实验题目：用双线性变换法设计Butterworth滤波器对加噪后的语音信号进行滤波。绘制频率响应曲线，画出滤波前后的时域图和频谱图。实验程序清单：中国地质大学（武汉）数字信号处理实验报告6%双线性变换法设计Butterworth滤波器fs=8000;x1=wavread('pb8k.wav');t=(0:length(x1)-1)/8000;f=fs*(0:1023)/2048;A1=0.05;A2=0.10;d=[A1*cos(2*pi*3800*t)+A2*sin(2*pi*3600*t)]';x2=x1+d;%用双线性变换法设计一数字巴特沃斯低通器fp=3200;fs=3400;Fs=8000;Ts=1/Fs;wp=2*pi*fp*Ts;ws=2*pi*fs*Ts;Rp=1;Rs=15;wp1=2/Ts*tan(wp/2);%将模拟截止频率进行预畸变ws1=2/Ts*tan(ws/2);[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');%选择滤波器的最小阶数[Z,P,K]=buttap(N);%创建butterworth模拟滤波器%输入N：滤波器的阶数%输出Z:零点；P：极点；K：增益[Bap,Aap]=zp2tf(Z,P,K);[b,a]=lp2lp(Bap,Aap,Wn);%求模拟滤波器的传递函数[bz,az]=bilinear(b,a,Fs);%用双线性变换法实现模拟滤波器到数字滤波器的转换（即求的数字滤波器的系统函数）[H,W]=freqz(bz,az);%绘制频率响应曲线figure(1)plot(W*Fs/(2*pi),abs(H))gridon;axistight;xlabel('频率（Hz）')ylabel('频率响应')title('Butterworth')f1=filter(bz,az,x2);figure(2)subplot(2,1,1)plot(t,x2)%画出滤波前的时域图gridon;axistight;title('滤波前的时域波形');subplot(2,1,2)plot(t,f1);%画出滤波后的时域图gridon;axistight;中国地质大学（武汉）数字信号处理实验报告7title('滤波后的时域波形');y3=fft(f1,2048);figure(3)y2=fft(x2,2048);subplot(2,1,1);plot(f,abs(y2(1:1024)));%画出滤波前的频谱图gridon;axistight;title('滤波前的频谱')xlabel('Hz');ylabel('幅度');subplot(2,1,2)plot(f,abs(y3(1:1024)));%画出滤波后的频谱图gridon;axistight;title('滤波后的频谱')xlabel('Hz');ylabel('幅度');实验波形图：中国地质大学（武汉）数字信号处理实验报告8实验小结对于IIR数字滤波器的设计还真是遇到很多问题，刚开始是找不到音频文件，寻求助教和同学的帮助后，自己慢慢一步一步操作才得到实验结果。主要是刚开始对于IIR数字滤波器原理不是很理解，导致设计滤波器来不够熟练。经过这次实验，可以说从新认识了IIR滤波器原理以及滤波器的设计原理，日后还需多加中国地质大学（武汉）数字信号处理实验报告9练习，熟练掌握IIR数字滤波器的设计。实验三：综合实验一、实验目的对之前的学习内容做到学以致用，通过实验能够设计滤波器并解决具体问题，能够利用数字信号处理语音或图像等实际问题。二、实验要求及内容给一段原始的语音信号（可以是自己录制的一段语音），加上一频率为3.8kHz的高频余弦噪声和频率为3.6kHz的高频正弦噪声（幅度自己可以选择），用窗函数设计一滤波器（要求最小阻带衰减为50dB）对加噪后的语音信号进行滤波，画出滤波器的频率响应曲线，画出滤波前后的时域图和频谱图。需要用到的函数：fir1用窗函数设计FIR滤波器的