matlab 快速傅立叶变换 FFT 实验报告.doc

南京邮电大学实 验 报 告实验名称_熟悉MATLAB环境_ 快速傅里叶变换（FFT）及其应用 IIR数字滤波器的设计 FIR数字滤波器的设计 课程名称 数字信号处理A 班级学号 B08020221_ 姓 名_孟祥熙_开课时间 2009 /2010 学年， 第 一学期实验一一、观察采样引起的混叠。（来源：p57 习题1-3）设模拟信号为EMBED Equation.3，t的单位为毫秒(ms)。设采样频率为3kHz，确定与EMBED Equation.3混叠的采样重建信号EMBED Equation.3。画出EMBED Equation.3和EMBED Equation.3在EMBED Equation.3范围内的连续波形。（因数字计算机无法真正画出连续波形，可用较密的离散点的连线来近似。）分别用EMBED Equation.3和EMBED Equation.3在两信号波形上标记出3kHz采样点。两信号波形是否相同？采样后的两序列是否相同？ 程序与分析：% problem 1% =%作用与汇编中的；以及C+编程中的/作用一样，起注释作用Clear%清屏% estimate x(t) and xa(t) with a much higher sampling freq. fs1%因数字计算机无法真正画出连续波形，可用较密的离散点的连线来近似。%因此，为画出x(n)与xa(t)的波形，取fs1=50KHZ。time_period=6; % unit: msfs1=50; % unit: kHzT1=1/fs1; % unit: msn1=0:fix(time_period/T1);%用两个冒号组成等增量语句，其格式为：t=初值：增量：终值。%如：t=0：0.02：0.08，t=0 0.02 0.04 0.06 0.08%当增量为1时，这个增量值可以略去，因而该语句只有一个冒号。%如：k=1：6，k=1 2 3 4 5 6%fix 向0舍入为整数，舍尾法x=cos(5*pi*n1*T1)+4*sin(2*pi*n1*T1).*sin(3*pi*n1*T1);xa=cos(pi*n1*T1);% obtain x(nT) and xa(nT) with given sampling freq. fsfs=3; T=1/fs; n=0:fix(time_period/T);x_sample=cos(5*pi*n*T)+4*sin(2*pi*n*T).*sin(3*pi*n*T);xa_sample=cos(pi*n*T);figure,plot(n1*T1,x,r,n1*T1,xa,b,n*T,x_sample,ro),%figure指令可以打开一个新的图形绘制窗口，在下一条这样的指令出现之前，%所有的绘图指令都会在这个窗口上绘制。%plot函数可以同时画多组数值图像，即以（X坐标数组1，Y坐标数组1，X%坐标数组2，Y坐标数组2，）的形式书写。每组坐标书写后，可用一个%表示前面的一组坐标用怎样的形式绘制：线型，点型和颜色。如r就表示%用绿色的线绘制，b蓝色。ro红色圆圈等等。hold on, stem(n*T,xa_sample,b:x)%hold保持当前图形窗的内容即：在画完一张图后用hold命令保持住，继续在%当前窗画图。%stem离散序列绘图%b:x表示蓝色虚线x号legend(x(t),xa(t),x(nT),xa(nT),xlabel(t(ms)%legend这个指令可以在图像的右上角绘出一个图例，表示每条曲线代表什么，%名称需编程者指定。%xlabel x轴标注，ylabel：y轴标注，title：三维坐标标注。上面程序结果如下图所示：波形分析：x(t)与xa(t)两波形并不相同，但采样后，波形一致。二、判别离散时间系统的时不变性。（来源：p105 例3.2.2）设输入序列为 EMBED Equation.3 ，系统 EMBED Equation.3 实现对 EMBED Equation.3 的抽取。设 EMBED Equation.3 。取延迟量D（例如D30）。记 EMBED Equation.3 ，画出 EMBED Equation.3 、 EMBED Equation.3 的序列波形。编程求出系统对 EMBED Equation.3 的响应 EMBED Equation.3 以及对 EMBED Equation.3 的响应 EMBED Equation.3 画出 EMBED Equation.3 、 EMBED Equation.3 的波形。该系统是否为时不变的？程序与分析：% =% problem 2% =clear% plot x(n) and x(n-D)D=30;N=500;n=1:N;%增量为1的等增量语句x=sin(2*pi/100*n);for n=1:N+D, if (n-D)=0, xD(n)=0; else xD(n)=x(n-D); endend%求xD(n)的表达式figure,subplot(2,1,1),%subplot(n,m,p),将图形分成n*m个子图，在第p个子图内绘图plot(1:N,x,r:,1:length(xD),xD,b),legend(x(n),xD(n),xlabel(n)% plot y(n) and yD(n)for n=1:fix(N/2)%fix,对N/2去尾取整 y(n)=x(2*n);endfor n=1:length(y)+D, if (n-D)=0, y_delay(n)=0; else y_delay(n)=y(n-D); endend%输入x(n)，输出延迟Dfor n=1:fix(length(xD)/2) yD(n)=xD(2*n);end%输入延迟D，对应的此时的输出subplot(2,1,2),%在第二个子图内绘图plot(1:length(y),y,r:,1:length(y_delay),y_delay,r.:,1:length(yD),yD,b.),legend(y(n),y(n-D),yD(n),xlabel(n)axis(0 530 -1 1)%axis(v),v四元向量xmin,xmax,ymin,ymax，坐标轴设定在v规定的范围内上面程序结果如下图所示：波形分析：根据系统时不变的性质，当xD(n)=x(n-D)时，有yD(n)= y(n-D)。由图2可以看出，y(n-D)与yD(n)的波形并不一致，则：系统为时变系统三、利用卷积计算信号通过FIR滤波器的输出，并观察输出信号的input-on暂态、input-off暂态和稳态阶段。（来源：p144 例4.1.8）考虑两个滤波器， EMBED Equation.3 ， EMBED Equation.3 ；输入 EMBED Equation.3 为周期方波，第一个周期内 EMBED Equation.3 。1.分别画出 EMBED Equation.3 通过两个滤波器的输出 EMBED Equation.3 、 EMBED Equation.3 的波形 EMBED Equation.3 ，并与书上p144例4.1.8的两幅图比较是否一致。2. 计算图中稳态部分的响应值。程序与分析：% =% problem 3% =clearh1=0.25*0.75.(0:14);h2=1/5*1 -5 10 -10 5 -1;%矩阵表示方法N=200;n=0:N-1;%增量为1的等增量语句x1=ones(1,25) zeros(1,25);% one period of x(n)%ones(m,n),全1矩阵(m行，n列);zeros(m,n),全0矩阵(m行，n列)%eye(n),单位矩阵，n阶方阵x=x1 x1 x1 x1;y1=conv(x,h1);y2=conv(x,h2);%y=conv(u,h)卷积语句figure,subplot(2,1,1),plot(n,x,r:,n,y1(1:N),b),%将图形分成2个子图，在第一个子图内绘图axis(0 200 -0.5 2.5),grid on,%grid：图上加坐标网格legend(input,output),xlabel(n)subplot(2,1,2),plot(n,x,r:,n,y2(1:N),b),axis(0 200 -1.5 2.5),grid on,%图上加坐标网格legend(input,output),xlabel(n)上面程序结果如下图所示：与书上P144图形一致2.1 EMBED Equation.3 Its DC gain is almost unity:(M=14,m=0)Ydc=h(m)=0.250.75m=1-0.75(14+1)=0.987 2.2 EMBED Equation.3 H(Z)=1/5(1-1/Z)5 ydc=H(Z)|z=1ydc=0实 验 二一、观察序列频谱，观察信号通过系统后波形与频谱的变化已知输入信号 EMBED Equation.3 ，其中 EMBED Equation.3 ， EMBED Equation.3 ， EMBED Equation.3 ，N可取5000点。（1）画出 EMBED Equation.3 的前100点波形（2）画出 EMBED Equation.3 的DTFT频谱 EMBED Equation.3 （ EMBED Equation.3 ） 由于计算机无法画出连续频谱，所以可在 EMBED Equation.3 内均匀取足够密的点数，如M5000个频率点 EMBED Equation.3 EMBED Equation.3 ，求出这些频率点上的频谱值 EMBED Equation.3 EMBED Equation.3 ，并画出 EMBED Equation.3 随 EMBED Equation.3 变化的曲线。（3）某LTI系统 EMBED Equation.3 ，画出系统的幅度频响 EMBED Equation.3 （4）求系统对 EMBED Equation.3 的响应 EMBED Equation.3 （可以自己编程也可利用卷积函数 EMBED Equation.3 ）。画出 EMBED Equation.3 的波形,并与 EMBED Equation.3 的波形比较（各画100点）；画出 EMBED Equation.3 的幅度谱 EMBED Equation.3 ，并与 EMBED Equation.3 比较。 你从中观察到什么？Promble 1二、系统函数 EMBED Equation.3 ，根据正准型结构（canonical form）编写样本处理算法。内部状态的初始值设为零，输入信号 EMBED Equation.3 采用逐个样本手动输入的方式，求输出信号 EMBED Equation.3 。Promble 2w1=0; for i=1:10, x=input(input x =); w0=0.8*w1+x; y=w0-0.5*w1; w1=w0; yend input x =0y = 0input x =1y = 1input x =2y = 2.3000input x =3y = 3.8400input x =4y = 5.5720input x =5y = 7.4576input x =6y = 9.4661input x =7y = 11.5729input x =8y = 13.7583input x =9y = 16.0066实验三理解例8.3.1（p393）IIR平滑器的作用和特点 （来源： p468 习题8.30）（1）重现p394图8.3.4 中a=0.90 和 a=0.98 两情形下L=200点的输出响应； EMBED Equation.3 ，其中s=5，噪声方差 EMBED Equation.3 相关程序：% - Problem 1 -clearL=200; %（用来确定信号的长度）s=5*ones(L,1); %（ones（）这个函数用来产生单位序列，高度为5）v=randn(L,1); %（randn（）是用来产生随机序列，喜爱MATLAB中这个函数是以均值为零方差为一的方式产生随机序列。）v=v-mean(v); %（下面这两句话是对产生的长度为200的序列进行修正，因为随机序列v虽然是以高斯白噪声产生的，不过由于其随机性，有可能会造成均值和方差不是0，1所以要分别减去均值mean(v)，方差std(v)）v=v/std(v); x=s+v;%（此时有用信号加上模拟出来的白噪声信号v，即是信道传输出来的信号）%（下面一步就是对滤波器的模拟，注意filter的用法：filter（b，【1 -a】，x）这是一种差分方程的表示方法y(n)-a*y(n-1)=b*x(n)用这种方法可以比较方便的表示一个差分方程及其响应）a=0.90;b=1-a;y1=filter(b,1 -a,x);%（注意！方框中不是1-a 而是1 -a）%（接下来就是对a=0.98的响应y2,并且进行画图）a=0.98;b=1-a;y2=filter(b,1 -a,x);subplot(2,1,1),plot(1:L),x,:,(1:L),y1,(1:L),s,r),legend(x(n),y(n),title(a=0.90)subplot(2,1,2),plot(1:L),x,:,(1:L),y2,(1:L),s,r),legend(x(n),y(n),title(a=0.98) Ps：图中的噪声信号和老师的信号有所不同原因是randn产生的噪声信号不同，不过总体趋势是一样的！（2）估算NRR。估算公式为： 输入噪声方差估值： EMBED Equation.3 ，输出噪声方差估值： EMBED Equation.3 NRR估值： EMBED Equation.3 mx=mean(x); varx=sum(x-mx).2)/L; % (请大家注意：（x-mx）后面有.是对矩阵运算时才愈要加上！根据公式估算出方差)start=1;%（从y=1这一点开始计算方差，会发现计算值与理论值有很大的误差）a=0.90;y1=y1(start:L);my=mean(y1);vary=sum(y1-my).2)/length(y1);NRResti=vary/varx;NRR=(1-a)/(1+a);a NRResti NRRa=0.98;y2=y2(start:L);my=mean(y2);vary=sum(y2-my).2)/length(y1);NRResti=vary/varx;NRR=(1-a)/(1+a);a NRResti NRR 根据以上的第一问的噪声信号所计算的理论值和估算值：（3）计算NRR理论值，与估计值比较。为何相差较大？如果想估计得准确，该怎么办？ （提示：增大L，增大程序中变量 start）Start=100时：Start=100；L=2000时：可以看出增大start可以减少与L的差别，而增大L则会使NRR减小！2（来源p469习题8.31 ）仍用上述一阶IIR平滑器，只是输入高频信号 EMBED Equation.3 b取何值时，信号部分通过系统后不变（指稳态）？如果你已经完成了作业，验证一下。代码如下：% - Problem 2 -clearL=200;s=5*(-1).(1:L); % clean signalv=randn(1,L); v=v-mean(v); v=v/std(v);% （与上题一样v是一个高斯白噪声信号）x=s; % （首先输入一个无杂波的信号，看系统响应，是否能够产生出输入信号）a=0.98;b=1+a; y=filter(b,1 -a,x);figure,plot(1:L),s,r:,(1:L),y),legend(s(n),y(n),axis(1 L -10 6),title(no noise) 从此可以看出，系统的冲击响应在与输入信号S(n)卷积后经过一段时间能还原出输入信号，即此滤波器可以实现原始信号的通过。（2）噪声部分是否放大了？然后再把噪声加进去：x=s+v;a=0.98;b=1+a; y=filter(b,1 -a,x);figure,plot(1:L),x,r:,(1:L),y),legend(x(n),y(n)这时，我们看到通过滤波器的信号已经无法恢复出输入信号了，证明滤波器无法抑制噪声信号，我们在进行编程看系统是否对噪声有放大作用：NRR=b.2/(1-a.2) NRR = 99.0000 系统对信号放大了99倍！M=2500;w=pi/M*(0:M-1); j=sqrt(-1);k=1:M,Hw=b./(1-a*exp(-j*w(k);%（求出H(w)即使用公式 EMBED Equation.KSEE3 * MERGEFORMAT 从而得到系统的频谱）H=abs(Hw).2;figure,plot(w/pi,H,w/pi,ones(M,1),r),legend(|H(w)|2,input noise spectrum),xlabel(omega(pi),ylabel(|H(omega)|2) 小结：通过这次试验，和自己的对照测试，查help觉得收获最大的就是了解了filter()这个函数的用法。很方便的根据差分方程或者z域里的方程直接编程写出，而且效率很高！实 验 四观察窗函数的影响。信号为 EMBED Equation.3 ， EMBED Equation.3 KHz， EMBED Equation.3 KHz， EMBED Equation.3 KHz， 采样频率 EMBED Equation.3 KHz。写出 EMBED Equation.3 （ EMBED Equation.3 ）的频谱 EMBED Equation.3 ；分别画出窗长度 EMBED Equation.3 ， EMBED Equation.3 ， EMBED Equation.3 ， EMBED Equation.3 的矩形窗频谱和Hamming窗频谱。观察L的变化对窗谱的主瓣宽度、旁瓣密集度、相对旁瓣水平的影响。 c) 时域采样点数分别取 EMBED Equation.3 ， EMBED Equation.3 ， EMBED Equation.3 ， EMBED Equation.3 ，分别画出 EMBED Equation.3 加矩形窗及加Hamming窗时DTFT频谱 EMBED Equation.3 ；你能否从频谱上分辨出信号的三个频率分量？若能分辨出，它们的位置和相对大小是否准确？L的大小对频率的物理分辨率（physical frequency resolution）有何影响？a):clear;N = 4000;fs = 10000;t = 0:1/fs:(N-1)/fs;w = 0:2*pi/N:(N-1)/N*2*pi;f1 = 2000;f2 = 2500;f3 = 3000;x = cos(2*pi*f1*t) + cos(2*pi*f2*t) + cos(2*pi*f3*t);X = fft(x);plot(w/pi,abs(X)/N);xlabel(omega in units of pi);ylabel(|X(omega)|) L = 10;subplot(2,2,1);windows_spectrum(L);L = 20;subplot(2,2,2);windows_spectrum(L);L = 40;subplot(2,2,3);windows_spectrum(L);L = 100;subplot(2,2,4);windows_spectrum(L); L=100;x=cos(0.4*pi*(0:(L-1)+cos(0.5*pi*(0:(L-1)+cos(0.6*pi*(0:(L-1);L = 10;subplot(2,2,1);windowed_spectrum(x(1:L);L = 20;subplot(2,2,2);windowed_spectrum(x(1:L);L = 40;subplot(2,2,3);windowed_spectrum(x(1:L);L = 100;subplot(2,2,4);windowed_spectrum(x(1:L); 由图形可看出，当L40时，汉明窗才可以分辨出三个峰值，而矩形窗只要L20即可。理解频率的物理分辨率和计算分辨率的区别信号同上，加矩形窗。时域采样点数分别取 EMBED Equation.3 ， EMBED Equation.3 ， EMBED Equation.3 ， EMBED Equation.3 。画出以上各种时长情况下，频域采样点数分别为 EMBED Equation.3 ， EMBED Equation.3 时的DFT（在同一个图上用虚线画出相应的DTFT频谱，用于比较）。离散频谱DFT和连续频谱DTFT有什么关系？L一定的情况下，能否通过增加N改善频率的物理分辨率？N的作用是什么？N0 = 32;L = 10;subplot(2,2,1);windowed_dft_spectrum(x(1:L),N0);L = 20;subplot(2,2,2);windowed_dft_spectrum(x(1:L),N0);L = 40;subplot(2,2,3);windowed_dft_spectrum(x(1:L),N0);L = 100;subplot(2,2,4);windowed_dft_spectrum(x(1:L),N0); DFT相当于DTFT在频域的采样。将两张图对比可发现，当L较小（如10），此时N增加也无效果，因为物理分辨率由L决定。N决定在频域采样的密度，也就是决定了计算分辨率。N0 = 64;L = 10;subplot(2,2,1);windowed_dft_spectrum(x(1:L),N0);L = 20;subplot(2,2,2);windowed_dft_spectrum(x(1:L),N0);L = 40;subplot(2,2,3);windowed_dft_spectrum(x(1:L),N0);L = 100;subplot(2,2,4);windowed_dft_spectrum(x(1:L),N0); 另附三个函数的源代码。由于时间所限，我减少了老师提供的程序之中的输入参数，使得这几个函数只提供对汉明窗和矩形窗的分析。以及附上的功能函数的有关解释：function windows_spectrum(L)% M = 1000;w = -pi:2*pi/M:(M-1)*pi/M;h_rect = rectwin(L);%这一句是以输入的参数L为长度进行矩形窗的截取H_rect_w = h_rect * exp(-j * (1:L)* w );%将上一句的得到的长度为L的数组进行频域的转换h_hamming = hamming(L);%同理是用hamming窗进行截取H_hamming_w = h_hamming * exp(-j * (1:L)* w );%再次进行转换plot(w/pi,abs(H_rect_w),w/pi,abs(H_hamming_w),:);%将hamming窗和矩形窗显示到图像上xlabel(omega in units of pi);title(L=,num2str(L);legend(rect,ham);function windowed_spectrum(x)%功能：对函数X进行加窗取样% L = length(x);%窗函数长度等于输入的函数长度M = 1000;w = -pi:2*pi/M:(M-1)*pi/M;h_rect = rectwin(L);%以长度L进行矩形窗函数的截取x_rect = h_rect.*x;%对于X乘以窗函数值的到加窗取样值X_rect_w = x_rect * exp(-j * (1:L)* w );%将函数值转换成频域的值h_hamming = hamming(L);%同样，取L长的hamming窗值x_hamming = h_hamming.*x;%进行相乘的到hamming窗的取样值X_hamming_w = x_hamming * exp(-j * (1:L)* w );%转换到频域中plot(w/pi,abs(X_rect_w),w/pi,abs(X_hamming_w),:);xlabel(omega in units of pi);title(L=,num2str(L);legend(rect,ham);function windowed_dft_spectrum(x,N)%以N长进行DFT% L = length(x);%加窗取样长度为L % 进行DTFT得到的X_rect_w为连续频域值M = 1000;w = 0:2*pi/M:(M-1)*pi/M;h_rect = rectwin(L);x_rect = h_rect.*x;X_rect_w = x_rect * exp(-j * (1:L)* w ); % N长的DFT对连续频域取样dw = 2*pi/N;dft_w = 0:dw:pi;X_rect_dft_w = x_rect * exp(-j * (1:L)* dft_w); plot(w/pi,abs(X_rect_w),:);%同时也画出连续频谱图进行对比hold on;stem(dft_w/pi,abs(X_rect_dft_w);xlabel(omega in units of pi);title(L = ,num2str(L), ,N = ,num2str(N);实验五窗口法设计FIR数字滤波器（来源：p571 习题10.18）(a) 用矩形窗设计长度分别为N=11、41、81、121的低通FIR滤波器，要求截止频率为 EMBED Equation.3 。画出滤波器的单位冲激响应 EMBED Equation.3 和幅度频响 EMBED Equation.3 曲线。问题：理想滤波器的频响是怎样的？当N增大时，FIR滤波器在 EMBED Equation.3 附近的最大纹波幅度是否降低？其余纹波的幅度是否减小？纹波的密度怎样变化？过渡带宽度怎样变化？（如有必要可增大N值观察）。在N=11时，画出滤波器的相频曲线。它是否是线性的？（b）用汉明窗再次设计同样的滤波器。用汉明窗设计出的滤波器与用矩形窗相比有什么特点？%=% Design FIR DF with window method% =wc=0.3*pi;% - rectangular windowwin=1; p=1; % observe magnitude responseN=11; figure,FIR_design(N,wc,win,2 2 1 2 2 2,p)N=41; FIR_design(N,wc,win,2 2 3 2 2 4,p)N=81; figure,FIR_design(N,wc,win,2 2 1 2 2 2,p)N=121; FIR_design(N,wc,win,2 2 3 2 2 4,p) 可以看出随着N的增大FIR滤波器在 EMBED Equation.3 处最大纹波幅度降低，其余纹波的幅度减小纹波的密度增高，过渡带宽度减少。p=2; % observe phase responseN=11; figure,FIR_design(N,wc,win,1 1 1 1 1 1,p)这就是对应的频率和相位响应。% - hamming windowwin=2; p=1;N=11; figure,FIR_design(N,wc,win,2 2 1 2 2 2,p) N=41; FIR_design(N,wc,win,2 2 3 2 2 4,p)N=81; figure,FIR_design(N,wc,win,2 2 1 2 2 2,p)N=121; FIR_design(N,wc,win,2 2 3 2 2 4,p)可见hamming窗的十分平滑，过渡带宽度却比矩形窗宽，符合我们理论讨论的结果。以Butterworth 模拟低通滤波器为原型，设计IIR数字滤波器。（a）截止频率 EMBED Equation.3 。设计11阶IIR数字低通滤波器，画出幅频、相频曲线。（b）截止频率 EMBED Equation.3 。设计11阶IIR数字高通滤波器，画出幅频、相频曲线。 问题：所设计的IIR滤波器与FIR滤波器的频率特性有何区别？N=11;Wn=wc/pi;b,a = butter(N,Wn); %默认为 lowpass， Wn为截止角频率为最后的输出形式，是以Z域表示滤波器的形式。H,w=freqz(b,a,8192);%同样使用freqz（）函数b为分子多项式，a为分母多项式合起来就是滤波器的表达式,H得到的是频率响应，w得到是相位变化，序列长度均为8192。figure,subplot(2,1,1),plot(w/pi,abs(H),xlabel(omega(pi),ylabel(|H(omega)|