dsp研究性学习报告

1 / 24DSP 研究性学习报告 数字信号处理课程研究性学习报告姓名 王佳学号 09211013同组成员 李禹霏指导教师 黄琳琳时间 XX 年 4 月 9 日DFT 近似计算信号频谱专题研讨【目的】掌握利用 DFT 近似计算不同类型信号频谱的原理和方法。理解误差产生的原因及减小误差的方法。培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】1.利用 DFT 分析 x=Acos+Bcos 的频谱，其中f1=100Hz，f2=120Hz。A=B=1; A=1,B=要求选择不同的 DFT 参数及窗函数，并对实验2 / 24结果进行比较,总结出选择合适 DFT 参数的原则.【题目分析】 频率分辨率在信号处理领域是否达到要求,会直接影响分析的结果。在最小采样率的条件下,可以用补零 DFT 来提高信号可视频率分辨率,但信号真正的频率分辨率并没有得到改善。通过对补零 DFT 和信号采样点数改变的分析研究,可以区分可视分辨率和真正频率分辨率两个重要的概念【结果分析】 对于窗函数的选择，应考虑被分析信号的性质与处理要求。如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用主瓣宽度比较窄而便于分辨的矩形窗，例如测量物体的自振频率等；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。【自主学习内容】 不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。信号的截断产生了能量泄漏，而用 FFT 算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的，但是我们可以通过选择不同的窗函数对它们的影响进行抑制。【阅读文献】【发现问题】3 / 24【问题探究】【仿真程序】N=30;L=512;f1=100;f2=120;fs=600;T=1/fs;t=*T;x=cos+cos;y=fft;mag=abs;f=-1)*fs/length;plot,mag);N=600;L=512;f1=100;f2=120;fs=600;T=1/fs;t=*T;x=cos+cos;4 / 24y=fft;mag=abs;f=-1)*fs/length;plot,mag);N=20;L=512;f1=100;f2=120;fs=600;T=1/fs;t=*T;x=cos+cos;y=fft;mag=abs;f=-1)*fs/length;plot,mag);A=1;B=1 用哈明窗N=30; %数据的长度L=512;%DFT 的点数f1=100;f2=120;fs=600;%抽样频率5 / 24T=1/fs; %抽样间隔ws=2*pi*fs;t=*T;x=cos+cos;wh=);X=x*wh;X=fftshift);w=*ws/L)/;plot)使用矩形窗N=30;L=512;f1=100;f2=120;fs=600;T=1/fs;ws=2*pi*fs;t=*T;x=cos+cos;wh=);X=x*wh;X=fftshift);6 / 24w=*ws/L)/;plot)A=1，B=矩形窗N=30; %y?Y?3?L=512;%DFT?yf1=100;f2=120;fs=600;%3?T=1/fs; %3?ws=2*pi*fs;t=*T;x=cos+*cos;wh=);X=x*wh;X=fftshift);w=*ws/L)/;plot)用哈明窗N=30; %y?Y?3?L=512;%DFT?yf1=100;7 / 24f2=120;fs=600;%3?T=1/fs; %3?ws=2*pi*fs;t=*T;x=cos+*cos;wh=);X=x*wh;X=fftshift);w=*ws/L)/;plot)【研讨题目】2.分析矩形窗、 Hann 窗、 Hamming 窗、 Blackman 窗、 Kaiser 窗的频谱，并进行比较。【题目分析】矩形窗矩形窗属于时间变量的零次幂窗。矩形窗使用最多，习惯上不加窗就是使信号通过了矩形窗。这种窗的优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。Hanning 窗汉宁窗又称升余弦窗，汉宁窗可以看作是 3 个8 / 24矩形时间窗的频谱之和，或者说是 3 个 sinc Hamming窗型函数之和，而括号中的两项相对于第一个谱窗向左、右各移动了 /T，从而使旁瓣互相抵消，消去高频干扰和漏能。可以看出，汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降Hamming 窗海明窗也是余弦窗的一种，又称改进的升余弦窗。海明窗与汉宁窗都是余弦窗，只是加权系数不同。海明窗加权的系数能使旁瓣达到更 凯泽窗小。分析表明，海明窗的第一旁瓣衰减为一42dB海明窗的频谱也是由 3 个矩形时窗的频谱合成，但其旁瓣衰减速度为 20dB，这比汉宁窗衰减速度慢。海明窗与汉宁窗都是很有用的窗函数。Kaiser 窗是一种可以任意调整窗谱主瓣宽度与旁瓣衰减比例的函数。以瞬时无功功率理论为基础，提出一种谐波电流的实时检测新方法，用 Kaiser 滑动时窗截取谐波电流信号，通过对窗函数参数的选定，提高了谐波电流检测的准确性。1矩形窗：利用 w=boxcar 的形式得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，它的元素由窗函数的值组成。 w=boxcar 等价于9 / 24w=ones.2.三角窗：利用 w=triang 的形式得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，它的元素由窗函数的值组成。 w=triang 等价于 bartlett。3汉宁窗：利用 w=hanning 得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，包含了窗函数的 n 个系数。4海明窗：利用 w=hamming 得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，包含了窗函数的 n 个系数。它和汉宁窗的主瓣宽度相同，但是它的旁瓣进一步被压低。5.布拉克曼窗：利用 w=blackman 得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，包含了窗函数的 n 个系数。它的主瓣宽度是矩形窗主瓣宽度的 3 倍，为 12*pi/N，但是它的最大旁瓣值比主瓣值低57dB。6切比雪夫窗：它是等波纹的，利用函数 w= 方式设计出 N 阶的切比雪夫 2 窗函数，函数的主瓣值比旁瓣值高 RdB,且旁瓣是等波纹的。7.巴特里特窗：利用 w=bartlett 的形式得到窗函数，其中 n 为窗函数的长度，而返回值 w 为一个 n 阶的向量，包含了窗函数的 n 个系数。10 / 248凯塞窗：利用 w=kaiser 的形式得到窗函数。【阅读文献】【发现问题】【问题探究】【仿真程序】1.矩形窗 调用格式：w=boxcar，根据长度 n 产生一个矩形窗 w。N=30;K=300;t=0:K-1;w=boxcar;Y=abs);plot;2.三角窗 调用格式：w=triang，根据长度 n 产生一个三角窗 w。 N=30;K=300;t=0:K-1;w=triang;Y=abs);plot;3.汉宁窗 调用格式：w=hanning，根据长度 n 产生一个汉宁窗 w。 N=30;11 / 24K=300;t=0:K-1;w=hanning;Y=abs);plot;4.海明窗 调用格式：w=hamming，根据长度 n 产生一个海明窗 w。 N=30;K=300;t=0:K-1;w=hamming;Y=abs);plot;5.布拉克曼窗调用格式：w=blackman，根据长度 n 产生一个布拉克曼窗 w。N=30;K=300;t=0:K-1;w=blackman;Y=abs);plot;6.恺撒窗 调用格式：w=kaiser，根据长度 n 和影响窗函数旁瓣的 参数产生一个恺撒窗 w。12 / 24N=30;K=300;t=0:K-1;w=kaiser;Y=abs);plot;【研讨题目】3.利用 DFT 分析 x=e -|t| u 的频谱。【题目分析】首先利用 CTFT 公式计算其模拟频谱的理论值；然后对其进行等间隔理想采样，得到 序列，利用 DTFT 公式计算采样序列的数字连续频谱理论值，通过与模拟频谱的理论值对比，理解混叠误差形成的原因及减小误差的措施. 接下来是对 序列进行加窗处理，得到有限长加窗序列 ，再次利用 DTFT 公式计算加窗后序列 的数字连续频谱，并与加窗前 的数字连续频谱进行对比，理解截断误差形成的原因及减小误差的措施；最后是对加窗序列进行 DFT 运算，得到加窗后序列 的 DFT 值，它是对 数字连续频谱进行等间隔采样的采样值.【仿真结果】该程序过程清晰、容易理解，程序运行结果如图 2 所示。图 2 是信号的时域波形，图 2 是对应的幅度谱图。由于在归一化频率为的地方还有较大幅度，所以对信13 / 24号进行理想采样后存在较大的混叠失真，表现在图 2 中归一化频率为-范围内的波形与图 2 的波形明显不同。图 2是理想采样序列加矩形窗得到的图形，对应的幅度谱线如图 2 中实线所示。该谱线与图 2 相比有明显的不同，这是加窗截断的结果。窗口长度越短截断效应会越明显。对加窗序列进行 DFT 运算，只要 DFT 点数大于等于窗口长度，算出的 DFT 值就是对加窗序列的连续频谱在一个周期内进行的等间隔采样的采样值。在图 2 中表现为 DFT 幅值在加窗序列连续幅谱的谱线上，是连续频谱曲线上的有限个数据点，即所谓栅栏效应。图 2 中用虚线画出了连续信号的幅谱理论值，与 DFT 的幅值对比，存在一定误差，但只要采样频率足够高，时窗长度足够长，FFT 点数足够大，得到的 DFT 值越逼近实际频谱。【自主学习内容】【阅读文献】【发现问题】【问题探究】【仿真程序】clc;close all;clear;%CTFT 程序，以 x=exp t=0 为例%利用数值运算计算并绘制连续信号波形L=4, %定义信号波形显示时间长度14 / 24fs=4,T=1/fs; %定义采样频率和采样周期t_num=linspace;%取若干时点,点数决定作图精度xt_num=exp;%计算信号在各时点的数值subplot;plot,%绘信号波形xlabel),ylabel),%加标签grid,title 信号时域波形),% 加网格和标题%利用符号运算和数值运算计算连续信号幅度谱的理论值 syms t W %定义时间和角频率符号对象xt=exp*heaviside,%连续信号解析式XW=fourier,%用完整调用格式计算其傅氏变换 %在0 两边取若干归一化频点，点数决定作图精度w1=linspace,linspace;XW_num=subs;%利用置换函数求频谱数值解 mag1=abs;%计算各频点频谱的幅值subplot;plot,%绘制归一化频率幅值谱线 xlabelrad/s),ylabel,%加标签 grid,title 连续信号幅频理论值),%加网格和标题 %DTFT 程序，以 x=exp n=0 为例%利用数值运算计算并绘制离散信号图形N=L*fs+1;n_num=0:N-1;%生成信号波形采样点15 / 24xn_num=exp;%计算信号理想采样后的序列值 subplot;stem,%绘序列图形 xlabel,ylabel),%加标签grid,title 理想采样图形),% 加网格和标题%利用符号运算和数值运算计算离散信号幅度谱的理论值 syms n z w %定义符号对象xn=exp, %定义离散信号Xz=ztrans,%用完整调用格式计算其 Z 变换%利用复合函数计算序列傅里叶变换的解析解Z=exp;Hejw=compose;Hejw_num=subs;%求频谱数值解mag2=abs;%计算各频点频谱的幅度subplot;plot,%绘制频谱幅度曲线 xlabelrad),ylabel,%加标签 grid,title 离散信号幅频理论值),%加网格和标题 %序列加窗图示及频谱幅值绘制%利用数值运算计算并绘制加窗后序列 xw 的图形M=8;win=);%定义窗点和窗型xwn=xn_num.*win,zeros;%给离散信号加窗subplot;stem,%加窗序列图示 xlabel,ylabel),%加标签grid,title 加窗序列图形),% 加网格和标题%利用符号运算和数值线16 / 24%利用 DFT 计算加窗序列 xw 的离散谱幅值Ndft=16, Xk=fft;%定义 DFT 点数和 DFT 运算 Xk0=fftshift*T;%将 DFT 值 0 对称和幅值加权处理 if mod=0; N1=Ndft; else N1=Ndft-1; end; k=0:-N1/2;wk=k/Ndft;%0 对称取值并归一化 hold on;stem,r.),%绘制 DFT 图形legend,%加响应图例，位置自动最佳 xlabel,ylabel,%加标签grid,title 加窗序列幅谱及其 DFT 幅值),plot,%与连续信号幅谱的理论值比较【研讨题目】4.利用 DFT 分析计算 x=cos2pit 的频谱,设w0=pi / 2【题目分析】 因为 f=w/2pi=1Hz , 所以抽样频率需要大于 2Hz .【序列频谱计算的基本方法】【仿真结果】信号时域波形 连续信号幅频理论值120-104617 / 24时间理想采样图形 28 幅度 x10频率 rad/s 离散信号幅频理论值20n加窗序列图形幅度频率xw【结果分析】【自主学习内容】已知一连续信号为x?cos?cos其中 f1?100Hz,f2?120Hz。试利用 DFT 分析其频谱。解：信号 x 的最高频率 fm?f2?120Hz，抽样频率 fs?2fm?240Hz，取抽样频率 fs?600Hz；最低的频率分辨率为信号样点数为 N?fs600?30。 18 / 24?fc20?fc?f2?f1?20Hz，最少的 N?30 的 MATLAB 程序如下：% program exa_1_,利用矩形窗计算有限长余弦信号频谱N=30; %数据的长度L=512;%DFT 的点数f1=100;f2=120;fs=600;%抽样频率T=1/fs;%抽样间隔t=*T;x=cos+cos;y=fft;mag=abs;f=-1)*fs/length;plot,mag);xlabel)ylabel程序运行结果如下图所示。由图可见，频谱图显示出两个较为明显的峰值。结论：当截取信号 N?30 样点时，频率分辨率?fc?fs600?20Hz，刚好能够 N3019 / 24分辨出 f1 和 f2 两个频谱分量，但频谱泄漏较严重。若取 N?20。N?20 的 MATLAB 程序如下：% program exa_1_,利用矩形窗计算有限长余弦信号频谱N=20; %数据的长度L=512;%DFT 的点数f1=100;f2=120;fs=600;%抽样频率T=1/fs;%抽样间隔t=*T;x=cos+cos;y=fft;mag=abs;f=-1)*fs/length;plot,mag);xlabel)ylabel程序运行结果如下图所示。由图可见，频谱图显示不出对应20 / 24f1?100Hz,f2?120Hz 的两个明显峰值。 结论：当截取信号N?20 样点时，频率分辨率?fc?fs600?30Hz，达不到最N20低的频率分辨率 f2?f1?20Hz，故分辨不出 f1和 f2 两个频谱分量，且频谱泄漏更为严重。若取频率分辨率?fc?1Hz，则对应的信号样点数为 N?N?600 的 MATLAB 程序如下： fs600?600。 ?fc1% program exa_1_,利用矩形窗计算有限长余弦信号频谱N=600; %数据的长度f1=100;f2=120;fs=600;%抽样频率T=1/fs;%抽样间隔t=*T;x=cos+cos;y=fft;mag=abs;f=-1)*fs/length;21 / 24plot/2),mag/2);xlabel)ylabel程序运行结果如下图所示。由图可见，频谱图显示出两个特别明显的峰值。结论：当截取信号 N?600 样点时，频率分辨率?fc?fs600?1Hz，高分辨率 N600的频谱图具有较高的质量，频谱分析时必须保证获取足够的信号数据长度。【阅读文献】【发现问题】【问题探究】【仿真程序】clc;close all;clear;%CTFT 程序，以 x=exp t=0 为例%利用数值运算计算并绘制连续信号波形L=8, %定义信号波形显示时间长度fs=4,T=1/fs; %定义采样频率和采样周期t_num=linspace;%取若干时点,点数决定作图精度xt_num=cos;%计算信号在各时点的数值subplot;plot,%绘信号波形22 / 24xlabel),ylabel),%加标签grid,title 信号时域波形),% 加网格和标题%利用符号运算和数值运算计算连续信号幅度谱的理论值syms t W %定义时间和角频率符号对象xt=cos*heaviside,%连续信号解析式XW=fouri