数字信号处理实验―用MATLAB设计循环卷积系统仿真

用MATLAB设计循环卷积系统仿真专业：电子信息工程学号：姓名：2014年12月一、实验目的（1）熟悉使用MATLAB软件。（2）学会调用MATLAB信号处理工具的设计函数。（3）对循环卷积有更深的认识和理解。二、实验原理和步骤1、卷积的定义：任意信号都可以根据不同需要进行不同的分解。如信号可以分解为直流分量和交流分量，也可以分解为奇分量和偶分量，或分解为实部分量和虚部分量。如果信号费解为冲击信号，那么信号分解为一系列不同强度，不同时延的冲击信号的叠加，这个过程称为卷积积分。一般而言，如果有两个函数tf1和tf2，则它们的积分dtffty21称为tf1与tf2的卷积积分，简称卷积，表达式为：tftfty21，即:dtfftftfty21212、线性卷积的运算：卷积运算是线性时不变系统分析的重要工具，很多滤波器的设计中都要用到卷积运算。给出线性卷积运算的定义，设有离散信号x(n)和y(n)，其线性卷积为：线性卷积有四步运算：①卷积运算时，y(n)要先反折得到y(-n);②m0表示y(-n)序列右移，m0表示左移，不同的m表示不同的值。线性卷积运算简洁表示为:式中“∗”表示线性卷积运算符。由线性卷积的定义，等式右边是乘积求和形式，，因而考虑能否用矩阵相乘的形式来表示线性卷积。假设序列x(n)长度为4点，y(n)长度为3点，x(n)除区间之外皆为零，y(n)除区间之外皆为零，用矩阵的形式来表达线性卷积Z:00000Z=0000000x(n),y(n)序列长度不同，则将短序列补0使两者相同。3、循环卷积的运算有限长序列的循环移位是指y((m-n))，也就是先让序列y(n)以N为周期进行周期延拓，再进行反折，然后朝右移位，只朝一个方向移位的原因是：对周期序列向右移动一个位置，也就相当于向左移动了N－1个位置,最后取（0，N－1）的N个值就得到了循环移位后的N个序列值。设有序列x(n)和y(n)，其N点循环卷积为：由于循环移位的关系最后得到的循环卷积的长度就是N点，m取[0，1，2，…，N-1]。循环卷积的简洁表示为：式中表示循环卷积运算符。例如N=4的循环卷积如下：Z=其中，N≥length(y(n))。值得说明的是，当N≥length(y(n))+length(x(n))-1时，圆周卷积的值等于线性卷积。三、MATLAB设计循环卷积1、循环卷积的分析：两个序列的循环卷积可以分三个步骤完成：(1)初始化：确定循环点数N，测量输入2个序列的长度，长度小于N的在后面补0。(2)循环右移函数：将序列x(n)循环右移，一共移N次（N为循环卷积的循环次数），最后将每次循环成的新序列组成一个矩阵V。(3)相乘：将x(n)移位后组成的矩阵V与第二个序列h(n)对应相乘，即得循环卷积结果。2、根据循环卷积分析设计流程图：循环卷积流程图如图1所示:（1）、主流程图(2)、循环函数流程图K=1图4.1循环卷积流程图循环卷积输入序列X1,X2输入循环点数N调用右移循环函数V=circlel（X2）Z=X1*V循环卷积结果ZENDX2(k)=X2(K-1)X2(1)=L;实现循环右移功能循环函数V=circlel（X2）A=length（X2）V=zeros(AA)i=1:Aj=1:AV(i.j)=X2(j)L=X2(A)K=A:-1:2当i=A（3）、据循环卷积流程图设计matlab源代码：functiony=myconv(x1,x2)x1=input('x1=');x2=input('x2=');N=input('N=');x1=[x1,zeros(1,N-length(x1))];x2=[x2,zeros(1,N-length(x2))];V=circlel(x2)Z=x1*V;stem(Z');xlabel('n');ylabel('Z');gridon;title('循环卷积结果Z')运行程序，输入序列x1，x2x1=[-123-5]x2=[67-10412]循环卷积结果：[10-5542-33-6986]运行图形如图2所示图4.2循环卷积运行结果四、实验总结总结本次数字信号处理实训，我受益匪浅。首先就是方案的确定。由于这个学期我学习了数字信号处理这门课程，课程中我了解到要实现两信号的卷积，可以通过定义来实现，也可以通过DFT来计算线性卷积。对于有限长序列，存在两种形式的卷积：线性卷积与圆周卷积。由于圆周卷积可以采用DFT的快速算法——快速傅里叶变换进行运算，运算速度上有很大的优越性。其中，设计线性卷积有4个步骤，反折、移位、相乘、求和，而循环卷积则是通过循环移位后得到的矩阵与序列相乘。根据课上老师讲的求法，画出了思路的流程图，然后根据流程图写出程序，事半功倍。将结果与直接调用matlab自带的函数比较，结果显示，自己设计的程序是正确的。通过这次实验，我对卷积和循环的定义、原理、以及实现方法都有了深入的认识。同时也对MATLAB软件产生了更加浓厚的兴趣。本次实训，检验了自己的能力，加强了逻辑思维的能力，不过我也发现了自身存在的一些问题，比如在MATLAB软件的应用上还有一些功能不懂如何运用的地方，但是在老师和同学的帮助下，我认真学习，并且懂得了许多以前不懂的matlab的运用。还有很多matlab的强大功能，希望能在日后好好学习，取得更好的成绩，也希望日后老师能不厌其烦的指导我，给予我更大的支持。五、参考文献[1]丁玉美，高西全。《数字信号处理》，第三版，西安电子科技大学出版社。[2]楼顺天，李博函.基于MATLAB的系统分析与设计一信号处理[M].西安:西安电子科技大学出版社，1998.81--88.用FFT对信号作频谱分析专业：电子信息工程学号：1208040230姓名：许先举2014年12月一、实验目的学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便正确用FFT。二、实验原理用FFT对信号作频谱的分析是学习信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变化区间N有关，因为FFT能够实现的频率分辨率是2π/N≤D。可以根据此式选择FFT的变换区间N。误差主要来自于用时，离散谱的包络才能逼近与连续谱，因此N要适当选择大一些。周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量的选择信号的观察时间长一点。对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的时间长度，经过采样后形成周期序列，按照周期序列的谱分析进行。三、实验步骤及内容（1）、对以下序列进行谱分析：x1(n)=R4(n)x2(n)=x3(n)=选择FFT的变换区间N为8和16两种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论、分析与比较。（2）、对以下周期序列进行谱分析：x4(n)=cos(π/4)*n；x5(n)=cos[(π/4)*n]+cos[(π/8)*n]选择FFT的变换区间N为8和16两种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论、分析与比较。n+10≤n≤38-n4≤n≤70其它n4-n0≤n≤3n-34≤n≤70其它n（3）、对模拟周期信号进行频谱分析：x6(n)=cos(8πt)+cos(16πt)+cos(20πt)选择采样频率Fs=64Hz，FFT的变换区间N为16、32、64三种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论和分析。四、思考题（1）、对于周期序列，如果周期不知道，如何用FFT进行谱分析？答：周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差要求就继续将截取长度加倍,重复比较,直到结果满足要求。（2）、如何选择FFT的变换区间（包括周期信号与非周期信号）答：一、对于非周期信号：有频谱分辨率F，而频谱分辨率直接和FFT的变换区间有关，因为FFT能够实现的频率分辨率是2π/N...因此有最小的N2π/F。就可以根据此式选择FFT的变换区间。二、对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。五、实验结果如下：x1(n)=R4(n)x2(n)=x3(n)=FFT的变换区间N为8和16两种情况进行频谱分析n+10≤n≤38-n4≤n≤70其它n4-n0≤n≤3n-34≤n≤70其它n实验结果图形与理论分析相符。（2）对以下周期序列进行谱分析：x4(n)=cos[(π/4)*n]x5(n)=cos[(π/4)*n]+cos[(π/8)*n]选择FFT的变换区间N为8和16两种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论、分析与比较。实验结果如下：（3）对模拟周期信号进行频谱分析：x6(n)=cos(8πt)+cos(16πt)+cos(20πt)选择采样频率Fs=64Hz，FFT的变换区间N为16、32、64三种情况进行频谱分析，分别打印出幅频特性曲线，并进行讨论、分析与比较。实验结果如下：六、实验中代码：x1n=[ones(1,4)];%产生R4(n)序列向量X1k8=fft(x1n,8);%计算x1n的8点DFTX1k16=fft(x1n,16);%计算x1n的16点DFT%以下绘制幅频特性曲线N=8;f=2/N*(0:N-1);figure(1);subplot(1,2,1);stem(f,abs(X1k8),'.');%绘制8点DFT的幅频特性图title('(1a)8点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');N=16;f=2/N*(0:N-1);subplot(1,2,2);stem(f,abs(X1k16),'.');%绘制8点DFT的幅频特性图title('(1a)16点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度');%x2n和x3nM=8;xa=1:(M/2);xb=(M/2):-1:1;x2n=[xa,xb];%产生长度为8的三角波序列x2(n)x3n=[xb,xa];X2k8=fft(x2n,8);X2k16=fft(x2n,16);X3k8=fft(x3n,8);X3k16=fft(x3n,16);figure(2);N=8;f=2/N*(0:N-1);subplot(2,2,1);stem(f,abs(X2k8),'.');%绘制8点DFT的幅频特性图title('(2a)8点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');subplot(2,2,3);stem(f,abs(X3k8),'.');%绘制8点DFT的幅频特性图title('(3a)8点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');N=16;f=2/N*(0:N-1);subplot(2,2,2);stem(f,abs(X2k16),'.');%绘制8点DFT的幅频特性图title('(2a)16点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度');subplot(2,2,4);stem(f,abs(X3k16),'.');%绘制8点DFT的幅频特性图title('(3a)16点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度');%x4n和x5nN=8;n=0:N-1;x4n=cos(pi*n/4);x5n=cos(pi*n/4)+cos(pi*n/8);X4k8=fft(x4n,8);X4k16=fft(x4n,16);X5k8=fft(x5n,8);X5k16=fft(x5n,16);figure(3);N=8;f=2/N*(0:N-1);subplot