北交大数字信号处理研究性学习报告
- 格式:doc
- 大小:7.56 MB
- 文档页数:67
《数字信号处理》课程研究性学习报告姓名学号同组成员指导教师时间多速率信号处理专题研讨【目的】(1) 掌握序列抽取运算与内插运算的频谱变化规律。
(2) 掌握确定抽取滤波器与内插滤波器的频率指标。
(3) 掌握有理数倍抽样率转换的原理及方法。
(4) 培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨题目】 基本题1.抽取、内插信号特征的时域/频域分析对于给定的单频模拟信号y (t )=sin(1000πt ),确定一个合适的采样频率f sam ,获得离散信号y [k ],试进行以下问题的分析:(1) 对离散信号y [k ]进行M=2倍抽取,对比分析y [k ]和y [M k ]在时域/频域的关系; (2) 对离散信号y [k ]进行L=2倍内插,对比分析y [k ]和y [k /L]在时域/频域的关系。
【温磬提示】在多速率信号分析中,离散序列的抽取和内插是多速率系统的基本运算,抽取运算将降低信号的抽样频率,内插运算将提高信号的抽样频率。
两种运算的变换域描述中,抽取运算可能出现频谱线性混叠,而内插运算将出现镜像频谱。
【设计步骤】1、 已知y (t )=sin(1000πt )频率为500Hz ,周期为0.002s ,可取时间范围T 为0到0.004秒,两个周期,根据抽样定理取Hz f sam 8000=,每个周期抽取16个点。
2、 用函数xD=x(1:M:end)对离散信号进行M=2倍的抽取,用fft 计算频谱。
3、 用函数xL=zeros(1,L*length(x));xL(1:L:end)=x;对离散信号进行L=2的内插,用fft 计算频谱。
【仿真结果】对离散信号y [k ]抽取和内插的时域/频域对比分析结果。
抽取:内插:【结果分析】1、抽取运算在频域描述:对x[k]进行M 倍抽取后得到][k x D 的频谱为∑-=-ΩΩ=102)(1)(M l Ml jj D eX Me X π,即将x[k]的频谱)(Ωj eX 扩展M 倍,得到)(Mj e X Ω,再以π2为周期进行周期化并乘以因子M1。
《信号与系统》课程研究性学习手册专题一信号时域分析1.基本信号的产生,语音的读取与播放【研讨内容】1) 生成一个正弦信号,改变正弦信号的角频率和初始相位,观察波形变化; 2) 生成一个幅度为1、基频为2Hz 、占空比为50%的周期方波,3) 观察一定时期内的股票上证指数变化,生成模拟其变化的指数信号, 4) 录制一段音频信号,进行音频信号的读取与播放 【题目分析】(1) 正弦信号的形式为Acos (ω0t +ψ)或Asin (ω0t+ψ),分别用MATLAB 的内部函数cos 和sin 表示,其调用形式为)*0cos(*phi t w A y +=、)*0sin(*phi t w A y += 。
生成正弦信号为y=5sin(t),再依次改变其角频率和初相,用matlab 进行仿真。
(2) 幅度为1,则方波振幅为0.5,基频w0=2Hz ,则周期T=pi ,占空比为50%,因此正负脉冲宽度比为1。
(3) 将波形相似的某一段构造成一个指数函数,在一连续时间内构造不同的2~3个不同指数函数即可大致模拟出其变化。
(4) 录制后将文件格式转化为wav ,再用wavread 函数读取并播放,用plot 函数绘制其时域波形。
【仿真】(1) 正弦信号 正弦信号1:A=1;w0=1/4*pi;phi=pi/16; t=-8:0.001:8;xt1=A*sin(w0*t+phi); plot(t,xt1)title('xt1=sin(0.25*pi*t+pi/16)')-8-6-4-22468-1-0.8-0.6-0.4-0.200.20.40.60.81xt1=sin(0.25*pi*t+pi/16)正弦信号2(改变1中频率)A=1;w1=1/4*pi;w2=1*pi;phi=pi/16; t=-8:0.001:8;xt1=A*sin(w1*t+phi); xt2=A*sin(w2*t+phi); plot(t,xt1,t,xt2)-8-6-4-22468-1-0.8-0.6-0.4-0.200.20.40.60.81正弦信号3(改变1中相位)A=1;w=1/4*pi;phi1=pi/16;phi2=pi/4; t=-8:0.001:8;xt1=A*sin(w*t+phi1); xt3=A*sin(w*t+phi2) plot(t,xt1,t,xt3)-8-6-4-22468-1-0.8-0.6-0.4-0.200.20.40.60.81(2) 方波信号 t=-100:0.01:100; T=0.5; f=1/T;y=square(2*pi*f*t,50); plot(t,y);axis([-2 2 -3 3]);-2-1.5-1-0.500.51 1.52-3-2-1123(3) 模拟股票上证指数变化的指数信号 x1=0:0.001:5;y1=2500+1.8*exp(x1); x2=5:0.001:10;y2=2847-1.5*exp(0.8*x2); x3=10:0.001:15;y3=2734+150*exp(-0.08*x3); x4=15:0.001:20;y4=2560-156*exp(-0.08*x4); x=[x1,x2,x3,x4]; y=[y1,y2,y3,y4]; plot(x,y);2468101214161820-2000-1500-1000-500050010001500200025003000(4) 音频信号的读取与播放[x,Fs,Bits]=wavread('C:\Users\Ghb\Desktop\nansheng.wav') sound(x,Fs,Bits) plot(x)00.51 1.52 2.53 3.54 4.55x 105-0.8-0.6-0.4-0.200.20.40.60.8[x,Fs,Bits]=wavread('C:\Users\Ghb\Desktop\nvsheng.wav') sound(x,Fs,Bits) plot(x)123456x 105-1-0.8-0.6-0.4-0.200.20.40.60.812.信号的基本运算(语音信号的翻转、展缩)【研讨内容】1) 将原始音频信号在时域上进行延展、压缩, 2) 将原始音频信号在频域上进行幅度放大与缩小, 3) 将原始音频信号在时域上进行翻转, 【题目分析】用matlab 的wavread 函数读取录制的音频,用length 函数计算出音频文件的长度,最后计算出时间t ,然后用plot 函数输出录制的音频信号 (1) 延展与压缩分析把时间t 变为原来的一半,信号就被延展为原来的2倍,把时间他变为原来的2倍,信号就被压缩为原来的一半。
《数字信号处理》课程研究性学习报告试点班专用姓名学号同组成员指导教师陈后金时间小波分析专题研讨【目的】(1) 掌握正交小波分析的基本原理。
(2) 学会Haar 小波分解和重建算法,理解小波分析的物理含义。
(3) 学会用Matlab 计算小波分解和重建。
(4) 了解小波压缩和去噪的基本原理和方法。
【研讨题目】 基本题【题目目的】:(1)掌握小波变换分解和重建算法的基本原理和计算方法; (2)掌握小波变换中Haar 基及其基本特性;8-1 (1)试求信号=T x [2, 2, 2, 4, 4, 4]T的Haar 小波一级变换系数]|[11T T d c 。
(2)将Haar 小波一级变换系数中的细节分量 1d 置零,试计算由系数]0|[1TT c 重建的近似信号1a , 求出x 与1a间的最大误差ε。
解:(1)}4,4,4,2,2,2{][][1↓+==k x k c j}21,21{][0↓=k h ,}21,21{][1-=↓k h }21,21{][0↓=-k h ,}21,21{][1↓-=-k h }5,4,3,2,1,0,1;22,24,24,23,22,22,2{][][][010-==-*=+k k h k c k y j }5,4,3,2,1,0,1;22,0,0,2,0,0,2{][][][111-=--=-*=+k k h k c k y j }2,1,0;24,23,22{]2[][0===k k y k c j }2,1,0;0,2,0{][=-=k k d j}0,2,0|24,23,22{]|[11-=→→TTd c(2)}4,3,2,1,0;24,0,23,0,22{]2[][0===k k c k y j}5,4,3,2,1,0;4,4,3,3,2,2{][][][00==*=k k h k y k x1=∴ε8-2 (1) 试求信号=T x [2, 2, 4, 6,−2,−2,−2, 0]T 的Haar 小波三级变换系数]|||[1233TT T T d d d c 。
《数字信号处理》课程研究性学习报告姓名学号同组成员指导教师时间DFT近似计算信号频谱专题研讨【目的】(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法。
(2) 理解误差产生的原因及减小误差的方法。
(3) 培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨题目】基本题1.D大调和弦由频率为293.66, 369.99, 440Hz的正弦信号组成,即x(t)= cos(2πf1t)+ cos(2πf2t)+ cos(2πf3t)的频谱,其中f1=293.66Hz,f2=369.99Hz,f3=440Hz。
(1)确定合适的DFT参数。
(2)利用DFT分析其频谱,并比较不同窗函数对谱分析结果的影响。
(3)若乐曲全音符的持续时间为0.2s,该和弦为16分音符,利用DFT分析其频谱会出现什么问题?【题目分析】【仿真结果】【结果分析】对实验结果进行比较,总结出选择合适DFT参数的原则。
【自主学习内容】【阅读文献】【发现问题】(专题研讨或相关知识点学习中发现的问题):【问题探究】【仿真程序】【研讨题目】 基本题2.已知某离散序列为 x [k ]=cos(Ω0k )+0.75cos(Ω1k ), 0≤ k ≤ 63 其中Ω0=0.4π, Ω1=Ω0+π/64(1) 对x [k ]做64点FFT, 画出此时信号的频谱。
(2) 如果(1)中显示的谱不能分辨两个谱峰,是否可对(1)中的64点信号补零而分辨出两个谱峰。
通过编程进行证实,并解释其原因 。
(3) 给出一种能分辨出信号中两个谱峰的计算方案,并进行仿真实验。
(M2-4)【温磬提示】在计算离散非周期序列频谱时常用Ω/π作为横坐标,称Ω/π为归一化频率(normalized frequency)。
在画频谱时需给出横坐标。
每幅图下都需给出简要的文字说明。
由于离散非周期序列频谱是周期的,所以在计算时不需要用fftshift 函数对fft 计算的结果进行重新排列。
《信号与系统》课程研究性学习手册专题一信号时域分析1. 基本信号的产生,语音的读取与播放【研讨内容】1) 生成一个正弦信号,改变正弦信号的角频率和初始相位,观察波形变化;2) 生成一个幅度为1、基频为2Hz 、占空比为50%的周期方波,3) 观察一定时期内的股票上证指数变化,生成模拟其变化的指数信号,4) 录制一段音频信号,进行音频信号的读取与播放【题目分析】⑴正弦信号的形式为Acosg o t+书)或Asin (3 o t+,分别用MATLAB 的内部函数cos 和sin 表示,其调用形式为y A* cos(w0* t phi)、y A*sin(wo*t phi)。
生成正弦信号为y=5sin(t), 再依次改变其角频率和初相,用matlab 进行仿真。
⑵幅度为1 ,则方波振幅为0.5 ,基频wO=2Hz ,则周期T=pi ,占空比为50% , 因此正负脉冲宽度比为 1 。
(3) 将波形相似的某一段构造成一个指数函数, 在一连续时间内构造不同的2~3 个不同指数函数即可大致模拟出其变化。
(4) 录制后将文件格式转化为wav ,再用wavread 函数读取并播放,用plot 函数绘制其时域波形。
【仿真】( 1 ) 正弦信号正弦信号 1 :A=1;w0=1/4*pi;phi=pi/16;t=-8:0.001:8;xt 仁A*si n(w0*t+phi);plot(t,xt1)title('xt 仁si n( 0.25*pi*t+pi/16)')正弦信号2 (改变1中频率)A=1;w1=1/4*pi;w2=1*pi;phi=pi/16; t=-8:0.001:8; xt 1= A*si n(w1*t+phi);xt2=A*si n(w2*t+phi);plot(t,xt1,t,xt2)正弦信号3 (改变1中相位)A=1;w=1/4*pi;phi仁pi/16;phi2=pi/4; t=-8:0.001:8; xt 1=A*si n(w*t+phi1);xt3=A*si n(w*t+phi2) plot(t,xt1,t,xt3)0.4 -0.2 -0 --0.2 --0.4 --0.6 --0.8 〜(2) 方波信号t=-100:0.01:100;T=0.5;f=1/T;y=square(2*pi*f*t,50);Plot(t,y);axis([-2 2 -3 3]);-3 1—--------- [ ------------ ■ ----------- 1- ---------- 1 ----------- 1 ----------- 1 ----------- 1 -------------------------t-2 -1.5 -1 -0.5 0 0.5 1 1.520.80.6-1 ------------- [ ---------- L-8 -6 -4(3) 模拟股票上证指数变化的指数信号x1=0:0.001:5;y1=2500+1.8*exp(x1);x2=5:0.001:10;y2=2847-1.5*exp(0.8*x2);x3=10:0.001:15;y3=2734+150*exp(-0.08*x3);x4=15:0.001:20;y4=2560-156*exp(-0.08*x4);x=[x1,x2,x3,x4];y=[y1,y2,y3,y4];plot(x,y);30002500200015001000500-500-1000-1500(4) 音频信号的读取与播放 [x,Fs,Bits]=wavread( sou nd(x,Fs,Bits) plot(x)-2000 ---------- [-------- [---------- L0 2 4 6 8 10 1214 16 18 20 'C:\Users\Ghb\Desktop\na nsheng.wav'C\Users\Ghb\Desktop\nvshe ng.wav' [x,Fs,Bits]=wavread(sou nd(x,Fs,Bits)plot(x)2. 信号的基本运算(语音信号的翻转、展缩)【研讨内容】1)将原始音频信号在时域上进行延展、压缩,2)将原始音频信号在频域上进行幅度放大与缩小,3)将原始音频信号在时域上进行翻转,【题目分析】用matlab 的wavread 函数读取录制的音频,用length 函数计算出音频文件的长度,最后计算出时间t ,然后用plot 函数输出录制的音频信号(1)延展与压缩分析把时间t 变为原来的一半,信号就被延展为原来的 2 倍,把时间他变为原来的 2 倍,信号就被压缩为原来的一半。
《数字信号处理》课程研究性学习报告试点班专用姓名学号同组成员指导教师陈后金时间DFT近似计算信号频谱专题研讨【目的】(1) 掌握利用DFT近似计算不同类型信号频谱的原理和方法。
(2) 理解误差产生的原因及减小误差的方法。
(3) 培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨题目】基本题1.利用DFT分析x(t)=A cos(2πf1t)+B cos(2πf2t)的频谱,其中f1=100Hz,f2=120Hz。
(1)A=B=1; (2)A=1,B=0.2。
【题目分析】分析题目,给出合适的DFT参数由取样定理知,要使信号频谱不混叠,则抽样频率不小于最高频率的两倍。
而要满足信号分辨率的要求,抽样点数N≧f sam/△f。
在对信号做DFT时,由于对信号进行截短,因此会产生频谱泄漏,要想从频谱中很好的分辨出个频率分量,需要考虑时域抽样频率,所加的窗函数,窗函数的长度,以及DFT的点数等参数对结果的影响。
(1)A=B=1,即x(t)=cos(2πf1t)+cos(2πf2t)矩形窗1:条件:fsam=240Hz;N=20;L=512矩形窗2:条件:fsam=600Hz;N=40;L=512矩形窗3:fsam=1200Hz;N=80;L=512Hamming窗1:N=40;L=512;fs=600;Hamming窗2:N=60;L=512;fs=600;Hamming 窗3:N=120;L=512;fs=600;(2)A=1,B=0.2,即x(t)=cos(2πf1t)+0.2cos(2πf2t)矩形窗:N=100;L=512;fs=600Hamming窗:N=100;L=512;fs=600【仿真结果】【结果分析】对实验结果进行分析比较,回答:加窗对谱分析有何影响?如何选择合适的窗函数?选择合适DFT 参数的原则?在(1)中进行矩形窗仿真时,我们选择了不同的fsam ,分别为240,600,1200它们均满足抽样定理,但是我们在实验中却发现,在240hz 时出现了混叠现象。
《数字信号处理》课程研究性学习报告姓名陶威东学号14212129同组成员张欣悦14212138指导教师刘留时间2015 年6月1日数字滤波器设计专题研讨【目的】(1) 掌握IIR和FIR数字滤波器的设计方法及各自的特点。
(2) 掌握各种窗函数的时频特性及对滤波器设计的影响。
(3) 培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨题目】基本题1.分析矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的频域特性,并进行比较。
【题目分析】不同窗函数的主瓣宽度,主瓣高度以及过渡带长度和阻带波动都有明显差别。
下面我们将通过对时域长度相同的窗函数进行分析来比较不同窗函数的频域特性。
取时域长度为50【仿真结果】【结果分析】各种窗有何特点?通带衰减由小到大的排列为:blackman,hamming,kaiser;阻带衰减由小到大的排列为:kaiser,hamming,blackman;过渡带由小到大的排列为:kaiser,hamming,blackman;在设计滤波器的时候,我们希望得到最大的通带衰减,最小的阻带衰减及最小的过渡带。
但在实践中,这三者是不可能同时满足的。
所以,我们要根据实际需求,选择合适的窗函数,以达到最佳的滤波效果。
最后我们还可以看出,主瓣宽度与主瓣幅值的乘积是定值,即主瓣宽度的增加必将导致主瓣幅值的降低。
最后我们会发现凯泽窗幅度响应曲线和矩形窗的完全吻合。
这是因为凯泽窗是是一个可以变幻形状的窗函数,当系统没有给值时,默认的为矩形窗,所以是完全符合的。
【自主学习内容】几种不同窗函数的设计【阅读文献】[1] 陈后金.数字信号处理[M].北京:高等教育出版社.2008.11【发现问题】(专题研讨或相关知识点学习中发现的问题):对于不同窗函数的选择上,我们需要根据需要在主瓣宽度和过渡带长短之间做出选择【问题探究】在谱分析中如何选择窗函数,在滤波器设计中如何选择窗函数?提高FIR滤波器阻带衰减是以增加过渡带宽度为代价的。
《数字信号处理》课程研究性学习报告DSP基本概念和技能的训练姓名张然学号13211074同组成员蔡逸飞13211078朱斌指导教师陈后金时间2015/6DSP 基本概念和技能研究性学习报告【目的】(1) 掌握离散信号和系统时域、频域和z 域分析中的基本方法和概念; (2) 学会用计算机进行离散信号和系统时域、频域和z 域分析。
(3) 培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨内容】问题一(1)阅读教材1.9节及MATLAB 中的Help ,学会MA TLAB 函数filter 的使用方法;(2)利用filter 函数,求出下列系统的单位脉冲响应,并判断系统是否稳定。
讨论实验所获得的结果。
211850586.0845.111)(--+-=z z z H21285.085.111)(--+-=z z z H 【题目目的】 1. 掌握LTI 系统单位脉冲响应的基本概念、系统稳定性与单位脉冲响应的关系; 2. 学会filter 函数的使用方法及用filter 函数计算系统单位脉冲响应; 3. 体验有限字长对系统特性的影响。
【仿真结果】 极点10.9430 0.9020 极点21.0000 0.8500051015202530354045502468y 1[k ]051015202530354045502468y 2[k ]【结果分析】我们所使用的计算机的是有限字长的,当我们用计算机对系统的各项参数进行量化,计算离散时,这些量化误差会使实际系统的极点值偏离理论值,导致系统的特性发生变化,甚至会使稳定系统变为非稳定系统。
【问题探究】已知LTI 系统的系统函数)(z H ,有哪些计算系统单位脉冲响应方法,比较这些方法的优缺点。
Filter 函数,可计算出差分方程的零状态响应,既可以用来求y[k],也可以求出h[k]; Impulse 函数,只是用来实现冲击响应的;Conv 函数,是用来计算卷积的,可以用来求y[k] 【仿真程序】 b1=[1 0 0]; b2=[1 0 0];a1=[1 -1.845 0.850586]; a2=[1 -1.85 0.85]; x=0:50;y1=filter(b1,a1,x); subplot(2,1,1); stem(y1);axis([0 50 0 8])[r1,p1,m1]=residuez(b1,a1); disp('极点1'); disp(p1');y2=filter(b2,a2,x); subplot(2,1,2); stem(y2);axis([0 50 0 8])[r2,p2,m2]=residuez(b2,a2); disp('极点2'); disp(p2');b1=[1 0 0]; b2=[1 0 0];a1=[1 -1.845 0.850586]; a2=[1 -1.85 0.85]; n=0:512;x=[1 zeros(1,512)] y1=filter(b1,a1,x); subplot(2,1,1); stem(n,y1); axis([0 50 0 8]) axis([0 50 0 8]) ylabel('y1[k]')[r1,p1,m1]=residuez(b1,a1); disp('极点1'); disp(p1');y2=filter(b2,a2,x); subplot(2,1,2); stem(n,y2); axis([0 50 0 8]) ylabel('y2[k]')[r2,p2,m2]=residuez(b2,a2); disp('极点2'); disp(p2');当取下列值时a1=[1 -1.8506 0.850586]; a2=[1 -1.85 0.906];极点11.0001 0.8505 极点20.9250 - 0.2244i 0.9250 + 0.2244i051015202530354045502468y 1[k ]5101520253035404550-505y 2[k ]问题二(1)阅读教材1.9节及MATLAB 中的Help ,学会MA TLAB 函数freqz 的使用方法; (2)利用MATLAB 语句x=firls(511,[0 0.4 0.404 1],[1 1 0 0]产生一个长度为512的序列x [k ],用plot 函数画出序列x [k ]的波形,用freqz 函数画出该序列的幅度频谱。
《数字信号处理》课程研究性学习报告数字信号处理综合应用专题研讨【目的】(1) 能够灵活应用DFT(FFT)分析实际信号的频谱。
(2) 熟悉通过IIR和FIR数字滤波器进行实际系统设计的方法。
(3) 培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
【研讨题目】1.附件myheart_noise1和myheart_noise2是一段含有噪声的音频信号。
(1)选择其中一个信号,分析其频谱特点。
(2)根据谱分析的结果,确定滤波器的指标。
(3)设计一个IIR DF对含有噪声的音频信号进行滤波。
(4)设计一个FIR DF对含有噪声的音频信号进行滤波。
(5)比较IIR DF和FIR DF对含有噪声音频信号的去噪效果。
(6) 请尝试采用其它的音频信号,混入不同的噪声,利用所学的滤波方法进行分析,会得到什么样的效果?【题目分析】本题讨论用IIR和FIR数字滤波器进行实际系统设计的方法。
对于第一小题,我们通过wavread函数播放音频,同时对该序列做DFT,求出其频谱图像。
对于第二小题,我们修改1中程序[s,f]=wavread('C:\Users\Lin\Desktop\matlab\myheart_noise1.wav');wavplay(s,f);L=length(s)m=fft(s);plot(abs(m));f=f得:L = 168873,f = 22050。
从图一可得,40000×22050/168873≈5223Hz,60000×22050/168873≈7834Hz。
取整可得,fp=6000Hz,fs=8000Hz.而后,设定AP =0.3dB,AS=40dB。
对于第三小题,先设计一个低通模拟滤波器,这里使用Butterworth型的低通滤波器,之后通过脉冲响应不变法和双线性变换法求得IIR数字低通滤波器。
对于第四小题,我们采用窗函数法设计FIR数字滤波器。
《数字信号处理》课程研究性学习报告姓名学号同组成员指导教师申艳时间2015年6月10日星期三题目一:基本概念和技能学习报告【目的】(1) 掌握离散信号和系统时域、频域和z 域分析中的基本方法和概念;(2) 学会用计算机进行离散信号和系统时域、频域和z 域分析。
(3) 培养学生自主学习能力,以及发现问题、分析问题和解决问题的能力。
M1-1 已知1()cos(6)g t t π=,2()cos(14)g t t π=,3()cos(26)g t t π=,以抽样频率10sam f Hz =对上述三个信号进行抽样,在同一张图上画出1()g t ,2()g t 和3()g t 及其抽样点,对所得结果进行讨论。
【题目目的】1. 掌握抽样的基本概念;2. 学会MATLAB 中对信号抽样的方法。
【仿真结果】【结果分析】通过仿真抽样前后的信号图像我们可知,g1=cos(6*pi*t),g2=cos(14*pi*t),g3=cos(26*pi*t)虽然它们信号时域表达式不同,但在以抽样频率fm=10Hz抽样后,各离散信号具有相同的时域表达式,所以当抽样频率过小时,所得离散信号可能存在较大失真,无法反映原信号的特征,所以在给定一个信号后,我们应给出合适的抽样频率,抽样频率过小时失真较多,抽样频率过大时又造成了浪费,抽样信号符合我们给定的要求即可,这样得到的信号才能在最大程度保留原信号的特征的基础上节约资源。
【仿真程序】k1=0:0.01:2;k2=0:0.1:2;g1=cos(6*pi*k1);g2=cos(6*pi*k2);subplot(3,2,1);plot(k1,g1);axis([0,2,-1.5,1.5]);title('g1=cos(6*pi*t) 原信号');xlabel('t');ylabel('g1(t)');subplot(3,2,2);stem(k2,g2);axis([0,2,-1.5,1.5]);title('g1=cos(6*pi*t) 抽样频率fm=10Hz');xlabel('k');ylabel('g1[k]');k3=0:0.01:2;k4=0:0.1:2;g3=cos(14*pi*k3);g4=cos(14*pi*k4);subplot(3,2,3);plot(k3,g3);axis([0,2,-1.5,1.5]);title('g2=cos(14*pi*t) 原信号');xlabel('t');ylabel('g2(t)');subplot(3,2,4);stem(k4,g4);axis([0,2,-1.5,1.5]);title('g2=cos(14*pi*t) 抽样频率fm=10Hz');xlabel('k');ylabel('g2[k]');k5=0:0.001:2;k6=0:0.1:2;g5=cos(26*pi*k5);g6=cos(26*pi*k6);subplot(3,2,5);plot(k5,g5);axis([0,2,-1.5,1.5]);title('g3=cos(26*pi*t) 原信号');xlabel('t');ylabel('g3(t)');subplot(3,2,6);stem(k6,g6);axis([0,2,-1.5,1.5]);title('g3=cos(26*pi*t) 抽样频率fm=10Hz');xlabel('k');ylabel('g3[k]');M1-2 利用MATLAB 的filter 函数,求出下列系统的单位脉冲响应,并判断系统是否稳定。
讨论本题所获得的结果。
211850586.0845.111)(--+-=z z z H 21285.085.111)(--+-=z z z H 【题目目的】1. 掌握LTI 系统单位脉冲响应的基本概念、系统稳定性与单位脉冲响应的关系;2. 学会filter 函数的使用方法及用filter 函数计算系统单位脉冲响应;3. 体验有限字长对系统特性的影响。
【仿真结果】下面是两个系统的单位脉冲响应:下面是两个系统的零极点分布图:极点10.9430 0.9020极点21.0000 0.8500【结果分析】根据以上仿真结果可知,两个系统函数H1,H2虽然只是分母系数差了一点,但最终实现时系统的单位脉冲响应会出现很大的不同,系统的极点也发生了不小的变化,这就造成了系统的失真。
比如,完成数字滤波器设计后,具体实现时,可能会出现较大的误差,甚至使设计的稳定系统变成不稳定系统,从而无法达到设计的要求。
引起这些误差的根本原因在于寄存器(存储单元)的字长有限。
误差的特性与系统的类型、结构形式、数字的表示法、运算方式及字的长短有关。
在通用计算机上,字长较长,量化步很小,量化误差不大。
但在专用硬件,如FPGA ,实现数字系统时,其字长较短,就必须考虑有限字长效应了。
【问题探究】已知LTI 系统的系统函数)(z H ,有哪些计算系统单位脉冲响应方法,比较这些方法的优缺点。
filter 函数,可计算出差分方程的零状态响应,既可以用来求y[k],也可以求出h[k],是最基本的方法,也是最复杂的方法;impulse 函数,只是用来实现冲击响应的,是最方便的方法;conv 函数,是用来计算卷积的,可以用来求y[k],这种方法也较麻烦,间接求出单位脉冲响应,不能直接求出单位脉冲响应。
【仿真程序】b1=[1];a1=[1,-1.845,0.850586];b2=[1];a2=[1,-1.85,0.85];x=[1,zeros(1,49)];y1=filter(b1,a1,x);subplot(2,1,1);stem(y1);xlabel('x');ylabel('h1[k]');axis([0 50 0 8]);y2=filter(b2,a2,x);subplot(2,1,2);stem(y2);xlabel('x');ylabel('h2[k]');subplot(2,1,1);[r1,p1,m1]=residuez(b1,a1);disp('极点1');disp(p1');zplane(b1,a1);title('H1');subplot(2,1,2);[r2,p2,m2]=residuez(b2,a2);disp('极点2');disp(p2');zplane(b2,a2);title('H2');M1-3 (1)利用MATLAB 语句x=firls(511,[0 0.4 0.404 1],[1 1 0 0])产生一个长度为512的序列x [k ],并画出该序列的幅度频谱。
(2) 已知序列)cos(][][0k k x k y Ω=,分别画出ππ,9.0π,8.0π,4.00=Ω时序列y [k ]的幅度频谱。
解释所得到的结果。
【题目目的】1. 学会用MATLAB函数freqz计算序列频谱;2. 掌握序列频谱的基本特性及分析方法。
【温磬提示】只需知道MATLAB语句x=firls(511,[0 0.4 0.404 1],[1 1 0 0])产生一个长度为512的序列x[k],该序列满足=-=kkkxx[,255[,1,0],]511不需知道其他细节。
用函数freqz计算该序列的频谱,在画幅度频谱时,建议用归一化频率。
【仿真结果】改进后:【结果分析】y[k]=)2/1))(ex p()](ex p([00k j k j k x Ω-+Ω对信号进行Fourier 变换:)]()()[2/1()()0()0(Ω-ΩΩ+ΩΩ+=j j j e X e X e Y所以,y[k]的幅度频谱是x[k]幅度频谱左移0Ω与右移0Ω之和的一半,在平移的过程中,会出现混叠现象,然而这样的理论分析与实际的结果不同,比如w0=0.8pi 时,本来值较高的地方实际的值却很低,原因如下:在x[k]变化较剧烈的点的值处与cos 函数值为零的值相乘后消去,所以将高频分量丢失,所以频谱会出现理论与实际结果不同的现象。
所以,是512的偶数倍的点均会出现这种现象,所以可以选择将x[k]变成3000个点的序列,即可解决这种现象。
【问题探究】有部分的计算结果可能与理论分析的结果不一致,分析出现该现象的原因,给出解决问题方法并进行仿真实验。
【仿真程序】x=firls(511,[0 0.4 0.404 1],[1 1 0 0]);b=x;a=1;w=linspace(-pi,pi,512);y=freqz(b,a,w);subplot(3,2,1:2);plot(w/pi,abs(y));axis([-2,2,-1,2]);title('生成序列x[k]的幅度频谱');k=0:511;y1=x.*cos(0.4*pi*k);g1=freqz(y1,a,w);subplot(3,2,3);plot(w/pi,abs(g1));xlabel('w0=0.4pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱'); y2=x.*cos(0.8*pi*k);g2=freqz(y2,a,w);subplot(3,2,4);plot(w/pi,abs(g2));xlabel('w0=0.8pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱'); y3=x.*cos(0.9*pi*k);g3=freqz(y3,a,w);subplot(3,2,5);plot(w/pi,abs(g3));xlabel('w0=0.9pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱'); y4=x.*cos(pi*k);g4=freqz(y4,a,w);subplot(3,2,6);plot(w/pi,abs(g4));xlabel('w0=pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱'); 改进后的程序:x=firls(2999,[0 0.4 0.404 1],[1 1 0 0]);b=x;a=1;w=linspace(-pi,pi,512);y=freqz(b,a,w);subplot(3,2,1);plot(w/pi,abs(y));axis([-2,2,-1,2]);title('生成序列x[k]的幅度频谱'); k=0:2999;y1=x.*cos(0.4*pi*k);subplot(3,2,2);plot(k,x);title('x[k]的时域图像');g1=freqz(y1,a,w);subplot(3,2,3);plot(w/pi,abs(g1));xlabel('w0=0.4pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱'); y2=x.*cos(0.8*pi*k);g2=freqz(y2,a,w);subplot(3,2,4);plot(w/pi,abs(g2));xlabel('w0=0.8pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱'); y3=x.*cos(0.9*pi*k);g3=freqz(y3,a,w);subplot(3,2,5);plot(w/pi,abs(g3));xlabel('w0=0.9pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱'); y4=x.*cos(pi*k);g4=freqz(y4,a,w);subplot(3,2,6);plot(w/pi,abs(g4));xlabel('w0=pi');ylabel('y[k]=x[k]cos(w0k)的幅度频谱');M1-4 已知111()(1)2H z z -=+,211()1H z z αα--=-,1α<,312()()()H z H z H z =,当0.8α=±,比较123(),(),()H z H z H z 的幅度响应。