信号处理实验四离散傅里叶变换

哈尔滨工程大学实验报告实验名称：实验四：离散傅里叶变换班级：电子信息工程4班学号：姓名：实验时间：2016年10月19日成绩：________________________________指导教师：栾晓明实验室名称：数字信号处理实验室哈尔滨工程大学实验室与资产管理处制实验四离散傅里叶变换一、实验原理1.由DFT定义式：10()[()]()NNknnXkDFTxnxnWk=0,1,…,N-1,将其写成矩阵方程表示为NXWx利用MATLAB的矩阵运算功能，可编写出计算DFT的函数文件。function[Xk]=dft(xn,N)%计算离散傅里叶变换%Xk=在0=k=N-1间的DFT系数数组%xn=N点有限长序列%N=DFT的长度n=[0:1:N-1];%n的行向量k=[0:1:N-1];%k的行向量WN=exp(-j*2*pi/N);%Wn因子nk=n'*k;%产生一个含bk值的N乘N维矩阵WNnk=WN.^nk;%DFT矩阵Xk=xn*WNnk;%DFT系数的行向量由IDFT定义式：101()[()]()NnkkxnIDFTXkXkWN，n=0,1,2,…,N-1，利用MATLAB的矩阵运算功能，可编写出计算傅里叶反变换的函数文件。function[xn]=idft(Xk,N)%计算离散傅里叶反变换%-----------------%xn=在0=n=N-1%Xk=N点有限长序列%N=IDFT的长度k=[0:1:N-1];%k的行向量n=[0:1:N-1];%n的行向量WN=exp(-j*2*pi/N);%Wn因子nk=n'*k;%产生一个含bk值的N乘N维矩阵WNnk=WN.^nk;%DFT矩阵xn=Xk*WNnk;%傅里叶反变换计算序列值DFT的快速算法FFT利用了WnkN的三个固有特性：(1)对称性，(W)WnknkNN，(2)周期性，()()，(3)可约性，WnknmkNmNW和//WnknkmNNmW。FFT算法基本上可以分为两大类，即按时间抽选法(DIT，Decimation-In-Time)和按频率抽选法(DIF，Decimation-In-Frequency)。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文件。2．利用DFT做连续信号的频谱分析DFT(实际中用FFT计算)可用来对连续信号和数字信号进行谱分析。在实际分析过程中，要对连续信号采样和截断，由此可能引起分析误差。(1)混叠效应对连续信号进行频谱分析时，首先要对其采样，变成时域离散信号后才能用DFT（FFT）进行谱分析。采样速率fs必须满足采样定理，否则会在w=π（对应模拟频率f=fs/2）附近发生频谱混叠现象。(2)截断效应处理实际信号序列x(n)时，一般总要将它截断为一有限长序列，长为N点，相当于乘以一个矩形窗形成有限长序列y(n)=x(n)w(n)。矩形窗函数其频谱有主瓣，也许许多副瓣，窗口越大，主瓣越窄，当窗口趋于无穷大时，就是一个冲击函数。时域的乘积对应于频域的卷积，所以，加窗后的频域实际是原信号频谱与矩形窗函数频谱的卷积，卷积的结果使频谱延伸到了主瓣以外，且一直延时到无穷。当窗口无穷大时，与冲击函数的卷积才是其本身，这是无需畸变。由此可见，阶段后频谱Y(ejw)与原序列频谱X(ejw)必然有差别，这种差别表现在：a.频谱泄露。原来序列x(n)的频谱是离散谱线，经截断后，原来离散谱线向附近展宽，成为泄露。显然，泄露使频谱变模糊，使谱分辨率降低。b.谱间干扰。主谱线两边又很多旁谱，引起不同频率分量间干扰，这使谱分析产生较大偏差。程度与窗函数幅度谱主瓣宽度直接相关。(3)栅栏效应N点DFT是频率区间[0,2π]上对时域离散信号的频谱进行N点等间隔采样，而采样点之间的频谱函数是看不见的。这就好像从N个栅栏缝隙中观看信号的频谱情况，仅得到N个缝隙中看到的函数值。由于栅栏效应，有可能漏掉大的频谱分量。例3.1对连续的单一频率周期信号按采样频率𝑓𝑠=8𝑓𝑎采样，截取长度N分别选N=20和N=16，观察其DFT结果的幅度谱。解：此时离散序列()sin(2)sin(28)asxnnffn，用MATLAB计算并作图，函数fft()用于计算离散傅里叶变换DFT，程序如下：k=8;n1=[0:19];%序列点数20xa1=sin(2*pi*n1/k);subplot(221)stem(n1,xa1)xlabel('t/T');ylabel('x(n)');%图一：序列图像xk1=fft(xa1);%调用fft函数xk1=abs(xk1);%取xk1绝对值subplot(222)stem(n1,xk1)xlabel('k');ylabel('X(k)');图二：序列DFT图像n2=[0:1:15];%序列点数16xa2=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)');二、实验内容1．已知连续周期信号()cos(10)2sin(18)xttt(1)确定信号的基频和基本周期.(2)当信号长度取0.5Tp,1.5Tp,2Tp,时,对x(t)采样，利用FFT计算其幅度谱；对所得结果进行比较，总结应如何选取分析长度。解：(1)基频=2，基本周期T=1s(2)采样频率𝑓𝑠=10𝑓𝑎采样，0.5Tp程序：f=10;%设置fs=10fan=[0:4];05101520-1-0.500.51t/Tx(n)0510152002468kX(k)051015-1-0.500.51t/Tx(n)0510150510kX(k)00.511.522.533.54-3-2-101t/Tx(n)采样结果00.511.522.533.540246kX(k)幅度谱02468101214161820-4-2024t/Tx(n)采样结果0246810121416182005101520kX(k)幅度谱%设置分析长度为0.5Tpxa=cos(10*pi*n/f)+2*sin(18*pi*n/f);%设置采样函数subplot(211)stem(n,xa)xlabel('t/T');ylabel('x(n)');%图一：采样的结果title('采样结果')xk=fft(xa);%进行FFT变换xk=abs(xk);%取X(k)的绝对值subplot(212)stem(n,xk)xlabel('k');ylabel('X(k)');%绘制出DFT结果title('幅度谱')采样频率𝑓𝑠=10𝑓𝑎采样，1.5Tp时将上述程序中n=[0,4]改为n=[0,14]即可采样频率𝑓𝑠=10𝑓𝑎采样，2Tp时将上述程序中n=[0,4]改为n=[0,19]即可结果分析：有三图比较可知，当长度为2Tp时，采样信号频谱不发生泄露，故应取2Tp2．对模拟信号()2sin(4)5cos(8)axttt以采样间隔T=0.01s采样，分别取N=40，N=50，N=60得到x(n)，用N点DFT得到𝑥𝑎(t)幅度谱的估计并比较结果解：T=0.01s,N=40时程序：T=0.01;%采样间隔02468101214-4-2024t/Tx(n)采样结果02468101214051015kX(k)幅度谱0510152025303540-10-50510t/Tx(n)0510152025303540050100kX(k)幅度谱N=40;%采样点数n=[0:N-1];xa=2*sin(4*pi*n*T)+5*cos(8*pi*n*T);%对信号进行采样xk=fft(xa);%进行FFT变换xk=abs(xk);%取X(k)的绝对值stem(n,xk)N=50时，将程序中N=40改为N=50即可N=60时，将程序中N=40改为N=60即可结果分析：比较三图可以看出，由于采样点数不同的栅栏效应，采样点越多栅栏效应对频谱分量的影响越小，谱分析误差越小。05101520253035404550-10-50510t/Tx(n)05101520253035404550050100150kX(k)幅度谱0102030405060-10-50510t/Tx(n)0102030405060050100150kX(k)幅度谱