哈工大数字信号处理实验报告
数字信号处理实验一 实验报告
数字信号处理实验一1.完成本文档内容的自学阅读和其中各例题后子问题；Q1.1运行程序P1.1，以产生单位样本序列u[n]并显示它。
答: clf;n=-10:20;u=[zeros(1,10) 1 zeros(1,20)];stem(n,u);xlabel('时间序号n');ylabel('振幅');title('单位样本序列');axis([-10 20 0 1.2])Q1.2命令clf，axis，title，xlabel和ylabel的作用是什么？答：clf清除图对象，axis 控制轴刻度和风格的高层指令，title 设置图名，xlabel和ylabel设置横纵坐标轴名称。
Q1.3修改程序P1.1以产生带有延时11个单位样本的延迟单位样本序列ud[n]。
运行修改的程序并显示产生的序列。
答：clf;n=0:30;ud=[zeros(1,11) 1 zeros(1,19)];stem(n,ud);xlabel('时间序号n');ylabel('振幅');title('单位样本序列');axis([0 30 0 1.2])Q1.4修改程序P1.1以产生单位步长序列s[n].运行修改后程序并显示产生的序列。
答：clf;n = 0:30;u = [1.*n];stem(n,u);title('Unit Sample Sequence');axis([0 30 0 30])Q1.5修改程序P1.1，以产生带有超前7个样本的延时单位阶跃序列sd[n]。
运行修改后的程序并显示产生的序列。
答：clf;n = -15:30;s=[zeros(1,8) ones(1,38)];stem(n,s);xlabel('Time index n');ylabel('Amplitude'); title('Unit Sample Sequence');axis([-15 30 0 1.2]);Q1.6 运行程序P1.2，以产生复数值的指数序列。
哈尔滨工程大学数字信号处理实验四
t/T x (n )k X (k)t/T x (n)k X (k )一、 实验原理DFT 的快速算法FFT 利用了WN^(nk)的三个固有特性：(1)对称性，(WN^(nk))*=WN^(-nk),(2)周期性，WN^(nk)=WN^(n+N)k=WN^n(k+K),(3)可约性WN^(nk)=WmN^(nmk)和WN(nk)=WN/m^(nk/m)。
FFT 算法基本上可以分为两大类，即按时间抽选法和按频率抽选法。
MATLAB 中提供了进行快速傅里叶变换的fft 函数，X=fft(x),基2时间抽取FFT 算法，x 是表示离散信号的向量；X 是系数向量； X=fft(x ，N)，补零或截断的N 点DFT ，当x 的长度小于N 时，对x 补零使其长度为N ，当x 的长度大于N 时，对x 截断使其长度为N 。
Ifft 函数计算IDFT ，其调用格式与fft 函数相同，参考help 文件。
例3.1程序及图形文件 >> k=8;>> n1=[0:19];>> xa1=sin(2*pi*n1/k); >> subplot(221) >> stem(n1,xa1)>> xlabel('t/T');ylabel('x(n)'); >> xk1=fft(xa1);>> xk1=abs(xk1);>> subplot(222)>> stem(n1,xk1)>> xlabel('k');ylabel('X(k)'); >> n2=[0:1:15]; >> xa2=sin(2*pi*n2/k); >> subplot(223) >> stem(n2,xa2)>> xlabel('t/T');ylabel('x(n)');>> xk2=fft(xa2);>> xk2=abs(xk2);>> subplot(224)>> stem(n2,xk2)>> xlabel('k');ylabel('X(k)');上两个图为N=20是的截取信号和DFT 结果，由于截取了两个半周期，频谱出现泄漏；下面的两个图为N=16时的截取信号和DFT 结果，由于截取了两个整周期，得到单一谱线的频谱。
哈工大数字信号处理实验2011
实验一 离散傅里叶变换的性质一、 实验目的1、 掌握离散傅里叶变换的性质，包括线性特性、时移特性、频移特性、对称性和循环卷积等性质；2、 通过编程验证傅里叶变换的性质，加强对傅里叶变换性质的认识。
二、 实验原理和方法 1. 线性特性1212DFT[()()]()()ax n bx n aX k bX k +=+2. 时移特性DFT[()]()DFT[()]()km kmx n m W X k x n m W X k -+=-=3. 频移特性()()nl N IDFT X k l IDFT X k W +=⎡⎤⎡⎤⎣⎦⎣⎦4. 对称性设由x(n) 延拓成的周期序列为 ()x n 则()()()e o xn x n x n =+ 共轭对称序列()()()*12e xn x n x N n ⎡⎤=+-⎣⎦ 共轭反对称序列()()()*12o x n x n x N n ⎡⎤=--⎣⎦ 将()e xn 和()o x n 截取主周期，分别得 ()()()ep e N x n x n R n = ()()()o p o N x n x n R n= 则()()()()()N ep op x n xn R n x n x n ==+ x(n)序列的实部和虚部的离散立叶变换(){}()Re ep DFT x n X k =⎡⎤⎣⎦ (){}()Im op DFT j x n X k =⎡⎤⎣⎦当x(n)为实数序列[][]()(())()()(())()()()(())()(())()()(())()(())()()()arg ()arg ()N N N N R R N N R N N I I N N I N N X k X k R k X k X N k R k X N k X k X k R k X N k R k X k X k R k X N k R k X k X N k X k X k *=-≅-=-≅-=-=-=--=--=-=--5. 循环卷积()312312()()()()()x n x n x n X k X k X k =⊗⇒=有限长序列线性卷积与循环卷积的关系 x1(n)和x2(n)的线性卷积：111212()()()()()N l m m x n x m x n m x m x n m -∞=-∞==-=-∑∑1120()()N m x m x n m -==-∑将x1(n)和x2(n)延拓成以N 为周期的周期序列11()()r xn x n rN ∞=-∞=+∑ 22()()q xn x n qN ∞=-∞=+∑ 则它们的周期卷积为14120()()()N p m x n xm x n m -==-∑ 1120()()N m x m xn m -==-∑ 1120()()N m q x m x n m qN -∞==-∞=-+∑∑1120()()N q m x m x n qN m ∞-=-∞=⎡⎤=+-⎢⎥⎣⎦∑∑ ()lq x n qN ∞=-∞=+∑x1(n)和x2(n)周期延拓后的周期卷积等于他们的线性卷积的的周期延拓。
数字信号处理实验报告_五个实验
实验一 信号、系统及系统响应一、 实验目的1、熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解；2、熟悉时域离散系统的时域特性；3、利用卷积方法观察分析系统的时域特性；4、掌握序列傅立叶变换的计算机实现方法，利用序列的傅立叶变换对连续信号、离散信号及系统响应进行频域分析。
二、 实验原理及方法采样是连续信号数字处理的第一个关键环节。
对采样过程的研究不仅可以了解采样前后信号时域和频域特性发生变化以及信号信息不丢失的条件，而且可以加深对傅立叶变换、Z 变换和序列傅立叶变换之间关系式的理解。
对一个连续信号)(t x a 进行理想采样的过程可用下式表示：)()()(^t p t t xx aa=其中)(^t x a 为)(t x a 的理想采样，p(t)为周期脉冲，即∑∞-∞=-=m nT t t p )()(δ)(^t x a的傅立叶变换为)]([1)(^s m a m j X T j a XΩ-Ω=Ω∑∞-∞=上式表明^)(Ωj Xa为)(Ωj Xa的周期延拓。
其延拓周期为采样角频率（T /2π=Ω）。
只有满足采样定理时，才不会发生频率混叠失真。
在实验时可以用序列的傅立叶变换来计算^)(Ωj X a 。
公式如下：Tw jw ae X j X Ω==Ω|)()(^离散信号和系统在时域均可用序列来表示。
为了在实验中观察分析各种序列的频域特性，通常对)(jw e X 在[0，2π]上进行M 点采样来观察分析。
对长度为N 的有限长序列x(n)，有：n jw N n jw k ke m x eX--=∑=)()(1其中，k Mk πω2=，k=0,1,……M-1 时域离散线性非移变系统的输入/输出关系为 ∑∞-∞=-==m m n h m x n h n x n y )()()(*)()(上述卷积运算也可在频域实现)()()(ωωωj j j e H e X eY =三、 实验程序s=yesinput(Please Select The Step Of Experiment:\n 一.(1时域采样序列分析 s=str2num(s); close all;Xb=impseq(0,0,1); Ha=stepseq(1,1,10);Hb=impseq(0,0,3)+2.5*impseq(1,0,3)+2.2*impseq(2,0,3)+impseq(3,0,3); i=0;while(s);%时域采样序列分析 if(s==1) l=1; k=0;while(1)if(k==0)A=yesinput('please input the Amplitude:\n',...444.128,[100,1000]); a=yesinput('please input the Attenuation Coefficient:\n',...222.144,[100,600]); w=yesinput('please input the Angle Frequence(rad/s):\n',...222.144,[100,600]); end k=k+1;fs=yesinput('please input the sample frequence:\n',...1000,[100,1200]); Xa=FF(A,a,w,fs); i=i+1;string+['fs=',num2str(fs)]; figure(i)DFT(Xa,50,string); 1=yesinput 1=str2num(1); end%系统和响应分析else if(s==2)kk=str2num(kk);while(kk)if(kk==1)m=conv(Xb,Hb);N=5;i=i+1;figure(i)string=('hb(n)');Hs=DFT(Hb,4,string);i=i+1;figure(i)string('xb(n)');DFT(Xb,2,string);string=('y(n)=xb(n)*hb(n)');else if (kk==2)m=conv(Ha,Ha);N=19;string=('y(n)=ha(n)*(ha(n)');else if (kk==3)Xc=stepseq(1,1,5);m=conv(Xc,Ha);N=14;string=('y(n)=xc(n)*ha(n)');endendendi=i+1;figure(i)DFT(m,N,string);kk=yesinputkk=str2num(kk);end卷积定理的验证else if(s==3)A=1;a=0.5;w=2,0734;fs=1;Xal=FF(A,a,w,fs);i=i+1;figure(i)string=('The xal(n)(A=1,a=0.4,T=1)'); [Xa,w]DFT(Xal,50,string);i=i+1;figure(i)string =('hb(n)');Hs=DFT(Hb,4,string);Ys=Xs.*Hs;y=conv(Xal,Hb);N=53;i=i+1;figure(i)string=('y(n)=xa(n)*hb(n)');[yy,w]=DFT(y,N,string);i=i+1;figure(i)subplot(2,2,1)plot(w/pi,abs(yy));axis([-2 2 0 2]);xlabel('w/pi');ylabel('|Ys(jw)|');title(FT[x(n)*h(n)]');subplot(2,2,3)plot(w/pi,abs(Ys));axis([-2 2 0 2]);xlabel('w/pi');ylabel('|Ys(jw)|');title('FT[xs(n)].FT[h(n)]');endendend子函数:离散傅立叶变换及X(n),FT[x(n)]的绘图函数function[c,l]=DFT(x,N,str)n=0:N-1;k=-200:200;w=(pi/100)*k;l=w;c=x*Xc=stepseq(1,1,5);子函数：产生信号function c=FF(A,a,w,fs)n=o:50-1;c=A*exp((-a)*n/fs).*sin(w*n/fs).*stepseq(0,0,49); 子函数：产生脉冲信号function [x,n]=impseq(n0,n1,n2)n=[n1:n2];x=[(n-n0)==0];子函数：产生矩形框信号function [x,n]=stepseq(n0,n1,n2) n=[n1:n2];x=[(n-n0>=0)];四、 实验内容及步骤1、认真复习采样理论，离散信号与系统，线性卷积，序列的傅立叶变换及性质等有关内容，阅读本实验原理与方法。
哈尔滨工业大学威海校区_《数字信号处理》实验一
哈尔滨⼯业⼤学威海校区_《数字信号处理》实验⼀数字信号处理实验报告实验名称：实验⼀离散傅⾥叶变换的性质实验⽇期：2011．11．16姓名：尤伟学号：090240328哈尔滨⼯业⼤学（威海）实验⼀离散傅⾥叶变换的性质⼀、实验⽬的1、掌握离散傅⾥叶变换的性质，包括线性特性、时移特性、频移特性、对称性和循环卷积等性质；2、通过编程验证傅⾥叶变换的性质，加强对傅⾥叶变换性质的认识。
⼆、实验原理和⽅法1.线性特性1212D FT [()()]()()ax n bx n aX k bX k +=+ 2.时移特性DFT[()]()DFT[()]()km kmx n m W X k x n m WX k -+=-=3.频移特性()()nlN IDFT X k l IDFT X k W +=4. 对称性设由x(n)开拓成的周期序列为 ()p x n 则()()()p pe po x n x n x n =+ 偶序列()()()*12pe p p x n x n x N n ??=+-?奇序列()()()*12pop p x n x n x N n ??=--?? 将()pe x n 和()po x n 截取主周期，分别得()()()pet pe N x n x n R n = ()()()p o tp oN x n x n R n =则()()()()()p N pet pot x n x n R n x n x n ==+ x(n)序列的实部和虚部的离散⽴叶变换(){}()R e petD FT x n X k = (){}()Im potj x n Xk =[][]()()()()()()()()()()()arg ()arg ()R R R I I I X k X k X N k X k X k X N k X k X k X N k X k X N k X k X k * =-=-=-=-=--=--=-=-- 5.循环卷积()3123121()()()()()x n x n x n X k X k X k N=?=有限长序列线性卷积与循环卷积的关系 X1(n)和x2(n)的线性卷积：11312120()()()()()N m m x n x m x n m x m x n m -∞=-∞==-=-∑∑112()()N m x m xn m -==-∑将X1(n)和x2(n)开拓成以N 为周期的周期序列11()()p r x n x n rN ∞=-∞=+∑22()()p q x n x n qN ∞=-∞=+∑则它们的周期卷积为14120()()()N p p p m x n xm x n m -==-∑12()()N p m x m xn m -==-∑1120()()N m q x m x n m qN -∞==-∞=-+∑∑1120()()N q m x m x n qN m ∞-=-∞=??=+-∑∑ 3()q x n qN ∞=-∞=+∑X1(n)和x2(n)周期开拓后的周期卷积等于他们的线性卷积的的周期开拓。
哈尔滨工程大学数字信号处理实验八
一、实验原理1、频率（周期）检测对周期信号来说，可以用时域波形分析来确定信号的周期，也就是计算相邻的两个信号波峰的时间差或过零点的时间差。
这里采用过零点（ti）的时间差T（周期）。
频率即为f=1/T,由于能够求得多个T的值（ti有多个），故采用它们的平均值作为周期的估计值。
2、幅值检测在一个周期内，求出信号最大值ymax与最小值ymin的差的一半，即A=（ymax-ymin）/2,同样，也会求出多个A值，但第一个A值对应的ymax和ymin不是在一个周期内搜索得到的，故以除以第一个以外的A值的平均作为幅值的估计值。
3、相位检测采用过零法，即通过判断与同频零相位信号过零点时刻，计算其时间差，然后换成相应的相位差。
φ=2п（１－ｔｉ／Ｔ），同样，以φ的平均值作为相位的估计值。
频率、幅值和相位估计的流程图见实验讲义64页图3.13。
4、数字信号统计量估计（1）峰值p的估计在样本数据x中找出最大值与最小值，其差值为双峰值，双峰值的一半即为峰值。
（2）均值估计（3）均方值估计(4) 差估估计5、频谱分析原理时域分析只能反映信号的幅值随时间的变化情况，除单频谱分量的简单波形外，很难明确提示信号的频率组成的各频率分量的大小，而频域分析能很好的解决此问题。
(1)DFT与FFT对于给定的时域信号y，可以通过Fourier变换得到频域信息Y。
采样信号是一个不连续的频谱，不能算出所有点的值，所以采用离散Fourier变换DFT，但其计算效率较低，故实际中常用快速傅立叶变换。
(2)频率、周期的估计对于Y(k△f)，如果当k△f=f~时，Y(k△f)取最大值，则f为频率的估计值，由于采样间隔的误差，f~也存在误差，其误差最大为△f/2。
周期T=1/f，从原理上可以看出，如果在标准信号中混有噪声，用上述方法仍能够精确地估计出原标准信号的频率和周期。
(3)频谱图为了直观地表示信号的频率特性，常常将傅立叶变换的结果用图形的方式表示，及频谱图。
数字信号处理实验报告
数字信号处理实验报告一、实验目的本次数字信号处理实验的主要目的是通过实际操作和观察，深入理解数字信号处理的基本概念和方法，掌握数字信号的采集、处理和分析技术，并能够运用所学知识解决实际问题。
二、实验设备与环境1、计算机一台，安装有 MATLAB 软件。
2、数据采集卡。
三、实验原理1、数字信号的表示与采样数字信号是在时间和幅度上都离散的信号，可以用数字序列来表示。
在采样过程中，根据奈奎斯特采样定理，为了能够准确地恢复原始信号，采样频率必须大于信号最高频率的两倍。
2、离散傅里叶变换（DFT）DFT 是将时域离散信号变换到频域的一种方法。
通过 DFT，可以得到信号的频谱特性，从而分析信号的频率成分。
3、数字滤波器数字滤波器是对数字信号进行滤波处理的系统，分为有限冲激响应（FIR）滤波器和无限冲激响应（IIR）滤波器。
FIR 滤波器具有线性相位特性，而 IIR 滤波器则在性能和实现复杂度上有一定的优势。
四、实验内容与步骤1、信号的采集与生成使用数据采集卡采集一段音频信号，或者在 MATLAB 中生成一个模拟信号，如正弦波、方波等。
2、信号的采样与重构对采集或生成的信号进行采样，然后通过插值算法重构原始信号，观察采样频率对重构信号质量的影响。
3、离散傅里叶变换对采样后的信号进行DFT 变换，得到其频谱，并分析频谱的特点。
4、数字滤波器的设计与实现（1）设计一个低通 FIR 滤波器，截止频率为给定值，观察滤波前后信号的频谱变化。
（2）设计一个高通 IIR 滤波器，截止频率为给定值，比较滤波前后信号的时域和频域特性。
五、实验结果与分析1、信号的采集与生成成功采集到一段音频信号，并在MATLAB 中生成了各种模拟信号，如正弦波、方波等。
通过观察这些信号的时域波形，对不同类型信号的特点有了直观的认识。
2、信号的采样与重构当采样频率足够高时，重构的信号能够较好地恢复原始信号的形状；当采样频率低于奈奎斯特频率时，重构信号出现了失真和混叠现象。
哈工程信号处理实验三
信号处理实验三实验要求：研究抽样过程，分析产生混叠效应的原因，实现不同的重建方案。
3.3.1>> n=0:80;>> f0=300;fs=8000;>> x=sin(2*pi*f0/fs*n+pi/3);>>stem(n,x)>>gtext('f0=300,fs=8k')>>plot(n,x)>>gtext('f0=300,fs=8k')>>gtext('Continue')b.>> n=0:80;>>fs=8000;>>subplot(221)>> x=sin(2*pi*100/fs*n+pi/3); >>stem(n,x);grid;gtext('f0=100') >>subplot(222)>> x=sin(2*pi*225/fs*n+pi/3); >>stem(n,x);grid;gtext('f0=225') >>subplot(223)>>stem(n,x);grid;gtext('f0=350') >>subplot(224)>> x=sin(2*pi*225/fs*n+pi/3); >>stem(n,x);grid;gtext('f0=475')d>> n=0:80;>>fs=8000;>>subplot(221)>> x=sin(2*pi*7525/fs*n+pi/3); >>stem(n,x);grid;gtext('f0=7525') >>subplot(222)>> x=sin(2*pi*7650/fs*n+pi/3); >>stem(n,x);grid;gtext('f0=7650') >>>>subplot(223)>> x=sin(2*pi*7775/fs*n+pi/3);>>stem(n,x);grid;gtext('f0=7675') >>subplot(224)>> x=sin(2*pi*7900/fs*n+pi/3); >>stem(n,x);grid;gtext('f0=7900')e>> n=0:80;fs=8000;>>subplot(221)x=sin(2*pi*32100/fs*n+pi/3); stem(n,x);grid;gtext('f0=32100') subplot(222)x=sin(2*pi*32225/fs*n+pi/3); stem(n,x);grid;gtext('f0=32225') subplot(223)x=sin(2*pi*32350/fs*n+pi/3); stem(n,x);grid;gtext('f0=32350') subplot(224)x=sin(2*pi*32475/fs*n+pi/3); stem(n,x);grid;gtext('f0=32475') >> t=0:1/80000:1000/80000; >> x=cos(2*pi*300*t+pi/6); >>stem(t,x)>> t=0:1/8000:1000/80000;>> x=cos(2*pi*300*t+pi/6);>>stem(t,x)>>fmagplot(x,1/80000)可以看到，当频率变化时，正弦信号并不是严格按信号频率变化的，而是呈周期的，有时变密有时变疏3.3.3>> t=0:1/80000:1000/80000;>> x=cos(2*pi*300*t+pi/6);>>stem(t,x)>> t=0:1/8000:1000/80000;>> x=cos(2*pi*300*t+pi/6);>>stem(t,x)>>fmagplot(x,1/80000)3.3.4>> n=0:100;>> t=n./80000;>> x=cos(2*pi*t*300+pi/6); >>stem(t./10,x)>> [X,W]=dtft(x,900);>>plot(W*10/pi,abs(X))3.3.5fs=8000;fsim=80000;fcut=2*(fs/2)/fsim;>> [b,a]=cheby2(9,60,fcut); %调用雪比滤波器>> [c,d]=freqz(b,a,1000,'whole');%1000点的频率响应>> d(501:1000)=d(501:1000)-2*pi; %d的位置全部左移2pi>> plot(d/2/pi/1000*fsim,abs(c));3.3.6fs=8000;fsim=80000;>> n=0:1000;>> l=length(n);>> x=cos(2*pi*300*n/fsim+pi/6);>> y=x(1:80000/8000:l);>> L=length(y); %取长度避免后面的维度不相等>> y2=zeros(1,length(x)); %置零序列>>for i=0:100;y2(i*10+1)=y(i+1); %取10倍数中的数，相当于抽样置零end>> t=0:0.001:1;>>plot(t,y2)>>fmagplot(y2,0.0000125)3.3.7fs=8000;fsim=80000;k=fsim/fs;n=0:1000;l=length(n);x=cos(2*pi*2000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('2k')fs=8000;fsim=80000;k=fsim/fs;n=0:1000;l=length(n);x=cos(2*pi*6000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)>>gtext('6k')x=cos(2*pi*7000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('7k')x=cos(2*pi*9000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('9k')x=cos(2*pi*10000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('10k')x=cos(2*pi*15000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('15k')x=cos(2*pi*17000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('17k')x=cos(2*pi*18000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('18k')x=cos(2*pi*19000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('19k')x=cos(2*pi*20000*n/fsim+pi/6);y=x(1:k:l);L=length(y);y2=zeros(1,length(x));for i=0:100;y2(i*10+1)=y(i+1);endt=0:0.001:1;y3=filter(b,a,y2); subplot(211);fmagplot(y3,0.0000125) subplot(212);plot(t,y3)gtext('20k')从20k开始出现了混叠。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验一： 用FFT 作谱分析实验目的：(1) 进一步加深DFT 算法原理和基本性质的理解(因为FFT 只是DFT 的一种快速算法， 所以FFT 的运算结果必然满足DFT 的基本性质)。
(2) 熟悉FFT 算法原理和FFT 子程序的应用。
(3) 学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT 。
实验原理：DFT 的运算量：一次完整的DFT 运算总共需要2N 次复数乘法和(1)N N -复数加法运算，因而直接计算DFT 时，乘法次数和加法次数都和2N 成正比，当N 很大时，运算量很客观的。
例如，当N=8时，DFT 运算需64位复数乘法，当N=1024时，DFT 运算需1048576次复数乘法。
而N 的取值可能会很大，因而寻找运算量的途径是很必要的。
FFT 算法原理：大多数减少离散傅里叶变换运算次数的方法都是基于nk N W 的对称性和周期性。
（1）对称性()*()k N n kn kn N N N W W W --==（2）周期性()(mod`)()()kn N kn n N k n k N N N N N W W W W ++===由此可得()()/2(/2)1n N k N n k nk N N N N N k N k N N W W W W W W ---+⎧==⎪=-⎨⎪=-⎩这样：1.利用第三个方程的这些特性，DFT 运算中有些项可以合并；2.利用nk N W 的对称性和周期性，可以将长序列的DFT 分解为短序列的DFT 。
前面已经说过，DFT 的运算量是与2N 成正比的，所以N 越小对计算越有利，因而小点数序列的DFT 比大点数序列的DFT 运算量要小。
快速傅里叶变换算法正是基于这样的基本思路而发展起来的，她的算法基本上可分成两大类，即按时间抽取法和按频率抽取法。
我们最常用的是2M N =的情况，该情况下的变换成为基2快速傅里叶变换。
完成一次完整的FFT 计算总共需要2log 2N N 次复数乘法运算和2log N N 次复数加法运算。
很明显，N 越大，FFT 的优点就越突出。
实验步骤(1) 复习DFT 的定义、 性质和用DFT 作谱分析的有关内容。
(2) 复习FFT 算法原理与编程思想， 并对照DIT-FFT 运算流图和程序框图， 读懂本实验提供的FFT 子程序。
(3) 编制信号产生子程序， 产生以下典型信号供谱分析用：1423()()1,03()8470403()347x n R n n n x n n n n n x n n n =⎧+≤≤⎪=-≤≤⎨⎪⎩-≤≤⎧⎪=-≤≤⎨⎪⎩456()cos4()sin 8()cos8cos16cos 20x n n x n n x t t t t πππππ===++(4) 编写主程序。
(5) 按实验内容要求， 上机实验， 并写出实验报告。
实验程序与结果:（1）对x 1(n)进行FFT 变换（N=8和N=16）clear all;N=4;%读入长度x=[1,1,1,1];%调用信号产生子程序产生实验信号n=0:N-1;subplot(3,1,1);stem(n,abs(x));%调用绘图子程序(函数)绘制时间序列波形图title('原时间序列');N=8;%读入长度y1=fft(x,N);%调用FFT 子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1)));holdon;plot(n,fftshift(abs(y1)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=8');N=16;%读入长度y2=fft(x,N);%调用FFT 子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2)));holdon;plot(n,fftshift(abs(y2)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=16');结果图形实验误差分析：理论FFT频谱为一个抽样信号，由图形看出，随着采样率的提高，得到的FFT 频谱分辨率就越高，当N趋于无限时频谱包络接近理论的抽样函数。
（2）对x2(n)进行FFT变换（N=8和N=16）clear all;x=[1,2,3,4,4,3,2,1];%调用信号产生子程序产生实验信号N=8;%读入长度n=0:N-1;subplot(3,1,1);stem(n,abs(x));%调用绘图子程序(函数)绘制时间序列波形图title('原时间序列');N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1)));holdon;plot(n,fftshift(abs(y1)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=8');N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2)));holdon;plot(n,fftshift(abs(y2)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=16');结果图形实验误差分析：理论FFT频谱为一个抽样函数平方的信号，由图形看出，随着采样率的提高，得到的FFT频谱分辨率就越高，当N趋于无限时频谱包络接近理论的抽样函数。
（3）对x3(n)进行FFT变换（N=8和N=16）clear all;x=[4,3,2,1,1,2,3,4];%调用信号产生子程序产生实验信号N=8;%读入长度n=0:N-1;subplot(3,1,1);stem(n,abs(x));%调用绘图子程序(函数)绘制时间序列波形图title('原时间序列');N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1)));holdon;plot(n,fftshift(abs(y1)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=8');N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2)));holdon;plot(n,fftshift(abs(y2)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=16');结果图形实验误差分析：由图形看出，随着采样率的提高，得到的FFT频谱分辨率就越高，但对于这个函数来说，N仍然过小，因此幅频特性曲线与理论结果有较大差距。
（4）对x4(n)进行FFT变换（N=8和N=16）clear all;n=0:15;%读入长度x=cos(pi/4*n);%调用信号产生子程序产生实验信号subplot(3,1,1);stem(n,abs(x));%调用绘图子程序(函数)绘制时间序列波形图title('cos(pi/4*n)');N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1)));holdon;plot(n,fftshift(abs(y1)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=8');N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2)));holdon;plot(n,fftshift(abs(y2)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=16');结果图形：实验误差分析：这个函数的理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，FFT 都可以与理论图形一样的曲线。
当N=8,16时，满足抽样定理。
（5）对x5(n)进行FFT变换（N=8和N=16）clear all;n=0:15;%读入长度x=sin(pi/8*n);%调用信号产生子程序产生实验信号subplot(3,1,1);stem(n,abs(x));%调用绘图子程序(函数)绘制时间序列波形图title('sin(pi/8*n)');N=8;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,2);stem(n,fftshift(abs(y1)));holdon;plot(n,fftshift(abs(y1)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=8');N=16;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(3,1,3);stem(n,fftshift(abs(y2)));holdon;plot(n,fftshift(abs(y2)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=16');结果图形：实验误差分析：这个函数的理论FFT抽样频谱为单位冲激函数，因此只要满足抽样定理，FFT 都可以与理论图形一样的曲线。
当N=8时，不满足抽样定理，因此产生了混叠失真，当N=8时，满足抽样定理，因此可得到与理论曲线一致的图形。
（6）对x6(n)进行FFT变换（N=16、N=32和N=64）clear all;n=0:64;%读入长度x=cos(pi/8*n)+cos(pi/4*n)+cos(pi/16*5*n);%调用信号产生子程序产生实验信号subplot(4,1,1);stem(n,abs(x));hold on;plot(n,abs(x),'r--');%调用绘图子程序(函数)绘制时间序列波形图title('cos(pi/8*n)+cos(pi/4*n)+cos(pi/16*5*n)');N=16;%读入长度y1=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(4,1,2);stem(n,fftshift(abs(y1)));holdon;plot(n,fftshift(abs(y1)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=16');N=32;%读入长度y2=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(4,1,3);stem(n,fftshift(abs(y2)));holdon;plot(n,fftshift(abs(y2)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=32');N=64;%读入长度y3=fft(x,N);%调用FFT子程序（函数）计算信号的DFTn=0:N-1;subplot(4,1,4);stem(n,fftshift(abs(y3)));holdon;plot(n,fftshift(abs(y3)),'r--');%调用绘图子程序(函数)绘制|X(k)|曲线title('N=64');结果图形：实验误差分析：这个函数的理论FFT抽样频谱为有三个不同频率分量的单位冲激函数，因此只要满足抽样定理，FFT都可以与理论图形一样的曲线。