数字信号处理论文

1、数字信号处理实验报告班级：1402014姓名：陈晓霞学号系方式安电子科技大学电子工程学院摘要本试验报告按照数字信号处理课程实验2017的内容及要求而进行的MATLAB仿真实验。共有三个实验，每个实验又有若干实验要求。其中第一个实验主要是实现常见序列的产生，序列运算的应用，以及复杂序列的产生和序列的线性卷积卷积；实验二主要是通过编程实现DFT ,IDFT,DIT-FFT,DIF-FFT,IFFT并利用DFT对连续信号分析；实验三则要求设计IIR以及FIR低通，带通，高通滤波器并对所给音频信号进行相应处理。报告主要包含解决上述内容的程序源代码，m

2、atlab处理图像，以及实验结论三部分。同时在报告所处的文件夹中含有相应.m文件以及处理后的音频信号。绪论 数字信号处理起源于十八世纪的数学，随着信息科学和计算机技术的迅速发展，数字信号处理的理论与应用得到迅速发展，形成一门极其重要的学科。当今数字信号处理的理论和方法已经得到长足的发展，成为数字化时代的重要支撑，其在各个学科和技术领域中的应用具有悠久的历史，已经渗透到我们生活和工作的各个方面。 数字信号处理相对于模拟信号处理具有许多优点，比如灵活性好，数字信号处理系统的性能取决于系统参数，这些参数很容易修改，并且数字系统可以分时复用，用一套数字系统可以分是处理多路信号；高精度和高稳定性，数字系

3、统的运算字符有足够高的精度，同时数字系统不会随使用环境的变化而变化，尤其使用了超大规模集成的DSP芯片，简化了设备，更提高了系统稳定性和可靠性；便于开发和升级，由于软件可以方便传送，复制和升级，系统的性能可以得到不断地改善；功能强，数字信号处理不仅能够完成一维信号的处理，还可以试下安多维信号的处理；便于大规模集成，数字部件具有高度的规范性，对电路参数要求不严格，容易大规模集成和生产。 从20世纪90年代末以来，数字信号处理课程几乎在国内所有大学的电气信息类等学科专业的本科生和研究生中开设，且是本科生的必修课和研究生的学位课。数字信号处理课程在经历了不断丰富发展的过程 后日臻完善。课程体系主要有

4、信号分析与处理，离散系统设计构成，同时增加了近代信号处理的理论和方法，并将MATLAB与数字信号处理有机结合，作为信号处理的仿真分析手段，从而将理论分析与计算机仿真融为一体。 信号处理在生物医学工程，地震学，声呐，雷达，通信，控制等领域都日益显示其重要作用。例如在医学信号或地震信号分析中，我们需要提取某些重要的特征参数，在雷达和通信信号处理中，我们希望提出信号中的噪声或干扰。 数字信号处理用途广泛，对其进行一系列学习与研究也是非常必要的。同时数字信号处理技术正在不断地克服自身的一些缺点，以惊人的速度发展，不断扩大应用领域，毫无疑问其前景非常好。实验一实验目的加深对序列基本知识的掌握理解实验原理

5、与方法1.几种常见的典型序列：2.序列运算的应用：数字信号处理中经常需要将被加性噪声污染的信号中移除噪声，假定信号 s(n)被噪声d(n)所污染，得到了一个含噪声的信号。我们的目的是对 x(n)进行运算，产生一个合理逼近 s(n)的信号 y(n)，从而实现去噪的目的。通常做法是对时刻 n 的样本附近的一些 样本进行求平均，产生输出信号，这是一种简单有效的办法，例如三点滑动平均算法的表达式如下：。3.复杂信号的产生：调幅广播（AM）是一种重要广播形式，使用的是振幅调制信号，其产生可用低频调制信号（如声音信号）来调制高频正弦信号（如载 波信号），不妨假设低频调制信号和高频正弦信号均为正弦信号形式，

6、具体如下： ，振幅调制信号为： 4.线性卷积：实验内容及步骤1.序列的产生：根据所学的知识产生单位阶跃序列，复指数序列，实指数序列，正弦序列，随机信号序列。 2.序列运算：设计一个混有白噪声的正弦信号，并通过三点滑动平均算法实现去噪，并通过时域图对比原正弦信号与去噪后获得的结果，观察去噪效果。3.复杂信号的产生：AM1323kHz是陕西交通广播，张旭东正在发出一个频率为 800Hz 的声音，试以上述信息为基础，利用matlab 工具产生此时交通广播台此时发出的广播信号。4.序列线性卷积验证实验：已知：试利用对位相乘法或者你熟悉的方法手动计算 y(n)，并通过matlab对手动计算结果进行验证，

7、画图说明。实验结果分析及结论总结 1序列的产生：a) 单位阶跃序列：n0=0;n1=-10;n2=10;n=n1:n2; %建立序列自变量向量x=(n-n0)>=0; %建立单位阶跃序列向量stem(n,x); %绘制单位阶跃序列曲线axis(-10 10 -1 1.5); %固定x与y轴刻度最大值及最小值xlabel('n');ylabel('x(n)'); %标注x与y轴代表参数title('单位阶跃序列') %添加曲线名称 图1.1 单位阶跃序列 b) 复指数序列：n=-10:1:10;A=1;alpha=-0.1+0.3j; x=A

8、*exp(alpha*n); %建立复指数序列向量 subplot(221);stem(n,real(x);title('实部');xlabel('n') subplot(222);stem(n,imag(x);title('虚部');xlabel('n')subplot(223);stem(n,abs(x);title('振幅 ');xlabel('n')subplot(224);stem(n,x);title('复指数序列');xlabel('n') %依次绘制复

9、指数序列的实部，虚部，振幅以及自身曲线图1.2 复指数序列c) 实指数序列：clear all;close all;clc;n=0:30;a1=0.8;a2=-0.8;a3=1.5;a4=-1.5; %取不同底数x1=a1.n;x2=a2.n;x3=a3.n;x4=a4.n; %建立实指数序列subplot(221);stem(n,x1);title('x=0.8n');xlabel('n') %绘制曲线subplot(222);stem(n,x2);title('x=(-0.8)n');xlabel('n')subplot(22

10、3);stem(n,x3);title('x=1.5n');xlabel('n')subplot(224);stem(n,x4);title('x=(-1.5)n');xla 图1.3 实指数序列 d) 正弦序列：n=-50:50;A=5;B=0;w=2*pi*0.01;x=A*sin(w*n+B); %建立正弦序列stem(n,x); %绘制曲线xlabel('n');ylabel('x(n)');title('正弦序列')图1.4 正弦序列e) 随机信号序列：n=1:20;A=1;x=A*ran

11、d(1,20); %建立随机序列stem(n,x); %绘制曲线xlabel('n');ylabel('x(n)');title('随机序列') 图1.5 随机信号序列2.序列运算：clear all；close all；clcR=51; d=2*(rand(1,R)-0.5); %建立加性噪声向量n=0:R-1; s=sin(pi/6)*n); %建立原始离散信号x=s+d; %建立被噪声影响后的信号subplot(2,1,1); plot(n,d,'r-',n,s,'g-',n,x,'b-.')

12、; %绘制加性噪声曲线xlabel('Time index n');ylabel('Amplitude');legend('dn','sn','xn'); %添加图标d添加标注x1=0 0 x;x2=0 x 0;x3=x 0 0;y=(x1+x2+x3)/3; %三点滑动平均算法subplot(2,1,2);plot(n,y(2:R+1),'r-',n,s,'g-'); %绘制滤出噪声后的信号legend('yn','sn');xlabel('

13、Time index n');ylabel('Amplitude');图1.6 除噪声前后信号 3.复杂信号的产生：clear allclose allclcR=30;n=1:R-1;A=2;m=0.5;x_l=cos(2*pi*0.8*n); %建立调制信号x_h=cos(2*pi*1323*n); %建立高频载波信号y(n)=A*(1+m.*x_l).*x_h; %建立调幅信号stem(n,y(n) %绘制调幅曲线title('调幅广播')图1.7 调幅广播4.序列线性卷积：clear all;N=5;M=6;L=M+N-1;x=1,2,3,4,5;

14、h=6,2,3,6,4,2; %建立卷积序列y=conv(x,h); %进行线性卷积运算nx=0:N-1;nh=0:M-1;ny=0:L-1; subplot(131); stem(nx,x);xlabel('n');ylabel('x(n)');grid on;title('x(n)');subplot(132);stem(nh,h);xlabel('n');ylabel('h(n)');grid on;title('h(n)');subplot(133);stem(ny,y);xlabel(&#

15、39;n');ylabel('y(n)');grid on;title('y(n)=x(n)*h(n)');图1.8 线性卷积验证结果：计算结果与对位相乘法计算所得结果一致。 思考题1) 可以通过什么参数控制序列的增长率与衰减率？哪个参数可以控制序列的幅值？图1.9 复指数序列振幅控制因素图1.10 复指数序列变化率控制因素通过MATLAB绘制图形结果得出：对于复指数序列，可以控制序列的幅值，控制序列的增长率与衰减率。2) 复指数序列的实部和虚部分别是什么？图1.11 复指数序列的实部与虚部特点 通过MATLAB绘制图形结果得出： 当时，复指数序列的实部

16、和虚部分别是按指数规律增长的正弦振荡序列；当时，复指数序列的实部和虚部分别是按指数规律衰减的正弦振荡序列；当时，复指数序列即为虚指数序列，其实部和虚部分别是按照等幅的正弦振荡序列。3) 正弦序列的频率是如何控制的？如何可以改变其频率？图1.12 正弦序列频率控制因素 通过MATLAB绘制图形结果可以得出： 正弦序列的频率是由参数控制的，并且可通过修改的取值来改变正弦序列的频率。4) 对于一个幅值为5v,频率为20Hz,初始相位为60度正弦连续信号，如何实现对其不失真采样？ 图1.13 正弦连续信号采样 通过MATLAB绘制图形结果可以看出：实现采样不失真的条件是采样周期满足：。（以上代码见qu

17、es1.m,ques12.m,ques2.m,ques3.m,ques4.m）实验二实验目的:加深对DFT和FFT的了解实验原理与方法:1. DFT的定义： 2.利用DFT对连续时间信号进行频谱分析： 连续时间信号的傅里叶变换所得信号的频谱函数是频率的连续函数，而序列的离散傅里叶变换是一种时域与频域均离散化了的变换，适合做数值运算，成为分析离散时间信号和系统的有力工具。而由于DFT具有选频特性，所以常用它对连续时间信号进行频谱分析。 利用DFT对连续时间信号进行频谱分析的过程如下图所示：图2.1 利用DFT对连续时间信号进行频谱分析的过程 从图中可以看出，利用DFT对连续时间信号进行频谱分析是

18、近似的，其近似程度与信号带宽，采样频率和截取长度等参数有关。 此外，利用DFT对连续时间信号进行频谱分析过程中会出现混叠现象，截断效应，栅栏效应以及信号的频谱分辨力等问题，与之对应分别有相应的解决方法，对此需要进行学习和理解。 3.FFT快速离散傅里叶变换： 序列的FFT算法就是把长序列的DFT逐步分解成几个短序列的DFT，并利用的周期性和对称性来减少DFT的运算次数。最常用的是基2FFT算法，它又分为时域抽取法和频域抽取法。 实验内容及步骤现有一信号s(t)，有四个频率分量1050KHz、1100KHz、1200KHz、1500KHz，请完成以下实验。注意：（1）s(t)请用实信号表示；（2

19、）不允许使用matlab中自带的DFT及相关函数，除非需要对所编写函数进行必要的验证。根据DFT的定义编写一个DFT函数和IDFT函数，并对s(n)进行DFT运算，对s(k)，即DFTs(n)进行IDFT运算，验证所编写函数的正确性。要求：s(k)与s(t)的频谱形式一致，即在-范围内对表示。对信号s(t)进行数字信号分析时，需要对DSP系统参数尽心谨慎的选择。（1）如果我们只检测到了5ms的信号，通过DFT操作是否可以将四个频率分量全部检测到，如果可以检测到，请作图说明，如果未检测到，有漏掉的频率分量，请继续完成以下问题；（2）如果你决定进行第2步实验，说明我们需要在参数选择上进行做出改变了

20、，通常我们会有以下措施：a）补零；b）改变采样率；c）增加被分析信号的有效长度。通过实验说明上述三种措施，哪些可以将被观测信号s(t)中的频率分量全部分析出来，哪些不可以。要提供理论分析和相关的实验仿真。1965年，J.W.Cooley与J.W.Tukey提出了DFT的快速算法，即FFT，大幅度的降低了DFT的运算量。（1）请自行编写一个DIT-基2FFT函数和DIF-基2FFT函数，对自定义的一个N=8的序列x(n)进行FFT操作，并对比两个函数的输出结果，最后通过matlab自带的FFT函数对自己所编写的两个函数进行验证。（2）基于（1）中的两个FFT函数，分别通过两个方法实现IFFT操作

21、，并通过与原x(n)的对比验证IFFT操作的正确性。实验结果分析及结论总结 据DFT的定义编写一个DFT函数和IDFT函数，并对s(n)进行DFT运算，对DFTs(n)进行IDFT运算，验证所编写函数的正确性。A. 主函数：clear all;close all;clc;dt=0.00001;tf=300;t=0:dt:tf; %tf决定t的取值范围，可更改st=sin(2*pi*1050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %连续信号s(t)subplot(221)plo

22、t(t,st);title('连续信号s(t)');xlabel('t');ylabel('s(t)'); %绘制连续信号图像Ts=1; %采样周期，可更改n=0:1:tf/Ts; %采样点n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+sin(2*pi*1500*n/20000); %采样信号s(n)subplot(222)stem(n,sn);title('采样信号s(n)');xlabel('n');yla

23、bel('s(n)'); %绘制采样信号图像Sk=dft(sn); %s(n)做DFT得到S(k)s1n=idft(Sk); %S(k)做IDFT得到s1(n)B. DFT函数：function Xk= dft(xn)%计算离散傅里叶变换%Xk=dft(xn)%Xk为在n在0,N-1区间的DFT系数数组%xn为N点有限持续时间序列N=length(xn);disp('输出N');disp(N); %有效长度N的输出n=0:1:N-1; %n的行向量k=0:1:N-1; %k的行向量W=exp(-j*2*pi/N); %Wn因子nk=n'*k; %产生一个

24、含nk值的N乘N维矩阵Xk=xn*(W.nk); %DFT系数的行向量disp('输出Xk');disp(Xk); %Sk的输出subplot(223)stem(k,Xk);title('s(n)做DFT所得S(k)');xlabel('k');ylabel('S(k)');%DFT的图像绘制endC. IDFT函数：function xn= idft(Xk)%计算逆离散傅里叶变换% xn=idft(Xk,N)%xn为n在0,N-1区间的N点有限持续时间序列%Xk为k在0,N-1区间的DFT系数数组N=length(Xk);dis

25、p('序列有效长度');disp(N);n=0:1:N-1; %n的行向量k=0:1:N-1; %k的行向量W=exp(-j*2*pi/N); %Wn因子nk=n'*k; %产生一个含nk值的N乘N维矩阵xn=Xk*(W.-nk)/N; %IDFT系数的行向量disp('输出xn');disp(xn); %sn的输出subplot(224)stem(n,xn);title('S(k)做IDFT所得s1(n)');xlabel('n');ylabel('s1(n)'); %IDFT图像绘制endD. 实验结果

26、：图2.2 DFT和IDFT图像E. 由图像可以看出对s(n)进行DFT得到的S（k）经过IDFT运算后所得s1(n)与序列s(n)一样，所以所编写的程序正确。（1） 只检测到了5ms的信号，通过DFT操作是否可以将四个频率分量全部检测到i. 代码 clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf决定t的取值范围，可更改st=sin(2*pi*1050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %连续信号s(t)

27、subplot(221)plot(t,st);title('连续信号s(t)');xlabel('t');ylabel('s(t)'); %绘制连续信号图像Ts=1; %采样周期，可更改n=0:1:tf/Ts; %采样点n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+sin(2*pi*1500*n/20000); %采样信号s(n)subplot(222)stem(n,sn);title('采样信号s(n)');xlabel(&

28、#39;n');ylabel('s(n)'); %绘制采样信号图像Sk=dft(sn); %s(n)做DFT得到S(k)ii. 图像图2.3 只检测到5ms信号的s(t)的DFTiii. 结论：不能够检测到，需进行第二问，（2） a）补零；b）改变采样率；c）增加被分析信号的有效长度。通过实验说明上述三种措施。a) 补零：clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf决定t的取值范围，可更改st=sin(2*pi*1050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*p

29、i*1200*t/20000)+sin(2*pi*1500*t/20000); %连续信号s(t)subplot(221)plot(t,st);title('连续信号s(t)');xlabel('t');ylabel('s(t)'); %绘制连续信号图像Ts=1; %采样周期，可更改n=0:1:tf/Ts; %采样点n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+sin(2*pi*1500*n/20000); %采样信号s(n)n1=0:1:t

30、f/Ts+1000;sn1=sn,zeros(1,1000);subplot(222)stem(n1,sn1);title('采样信号s(n)');xlabel('n');ylabel('s(n)'); %绘制采样信号图像Sk=dft(sn1); %s(n)做DFT得到S(k)title('补零所得图像')图2.4 补零所得图像b) 改变采样频率：clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf决定t的取值范围，可更改st=sin(2*pi*1050*t/2000

31、0)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %连续信号s(t)subplot(221)plot(t,st);title('连续信号s(t)');xlabel('t');ylabel('s(t)'); %绘制连续信号图像Ts=0.00001; %采样周期，可更改n=0:1:tf/Ts; %采样点n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000)+si

32、n(2*pi*1500*n/20000); %采样信号s(n)subplot(222)stem(n,sn);title('采样信号s(n)');xlabel('n');ylabel('s(n)'); %绘制采样信号图像Sk=dft(sn); %s(n)做DFT得到S(k)title('改变采样频率所得图像')图2.5 改变采样频率所得图像c) 增加被分析信号有效长度：clear all;close all;clc;dt=0.00001;tf=0.005;t=0:dt:tf; %tf决定t的取值范围，可更改st=sin(2*pi*1

33、050*t/20000)+sin(2*pi*1100*t/20000)+sin(2*pi*1200*t/20000)+sin(2*pi*1500*t/20000); %连续信号s(t)subplot(221)plot(t,st);title('连续信号s(t)');xlabel('t');ylabel('s(t)'); %绘制连续信号图像Ts=1; %采样周期，可更改n=0:1:tf/Ts; %采样点n的向量sn=sin(2*pi*1050*n/20000)+sin(2*pi*1100*n/20000)+sin(2*pi*1200*n/20000

34、)+sin(2*pi*1500*n/20000); %采样信号s(n)subplot(222)stem(n,sn);title('采样信号s(n)');xlabel('n');ylabel('s(n)'); %绘制采样信号图像N=length(sn)+1000;n=0:1:N-1; %n的行向量k=0:1:N-1; %k的行向量W=exp(-j*2*pi/N); %Wn因子nk=n'*k; %产生一个含nk值的N乘N维矩阵Sk=sn*(W.nk); %DFT系数的行向量disp('输出Xk');disp(Sk); %Sk的

35、输出subplot(223)stem(k,Sk);title('s(n)做DFT所得S(k)');xlabel('k');ylabel('S(k)');%DFT的图像绘制title('增加被分析信号有效长度所得图像') 图2.6 增加被分析信号有效长度所得图像（1） 编写一个DIT-基2FFT函数和DIF-基2FFT函数，对自定义的一个N=8的序列x(n)进行FFT操作，并对比两个函数的输出结果，最后通过matlab自带的FFT函数对自己所编写的两个函数进行验证。i. DIT-基2FFT函数a) 代码：clc;close all;

36、clear;format compact; xn=1 3 5 9 6 7 0 3; %初始序列xnM=nextpow2(length(xn), N=2M, for m=0:N/2-1; %旋转因子指数范围 WN(m+1)=exp(-j*2*pi/N)m; %计算旋转因子 end A=xn,zeros(1,N-length(xn); %数据输入，长度不够补零disp('输入到各存储单元的数据:'),disp(A); J=N/2; %给倒序数赋初值 for I=1:N/2-1; %按序交换数据和算倒序数 if I<J; %条件判断及数据交换 T=A(I+1);A(I+1)=A

37、(J+1);A(J+1)=T; end K=N/2; while J>=K; J=J-K;K=K/2; end J=J+K; end %disp('倒序后各存储单元的数据:'),disp(A); for L=1:M; %分级计算 % disp('运算级次:'),disp(L); B=2(L-1); for R=0:B-1; %各级按序蝶算 P=2(M-L)*R; for K=R:2L:N-2; %每序依次计算 T=A(K+1)+A(K+B+1)*WN(P+1); A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1); A(K+1)=T; end

38、 end %disp('本级运算后各存储单元的数据:'),disp(A); end %disp('输出各存储单元的数据:')Xk1=A;subplot(211)stem(Xk1); title('8点DITFFT计算X(k)')%disp('调用fft函数运算的结果:'),Xk2=fft(xn,N);subplot(212)stem(Xk2);title('对结果进行FFT验证X(k)')b) 图像结果：图2.7 DIT-基2FFTc) 结论：两者计算结果一致，故程序正确。ii. DIF-基2FFT函数a) 代码：

39、clear allclose allclcxn=1 3 5 9 6 7 0 3; %初始序列xnN=8;n=0:N-1;m=log2(N); %蝶形运算的次数d=bin2dec(fliplr(dec2bin(0:N-1,m)+1; %将其倒置排序Xk = xn(d); %将到序排列的x作为Xk初始值for M=1:m % 将DFT作m次基2分解,从左到右，对每次分解作DFT运算 Nm=2M;u=1; % 旋转因子u初始化为WN0=1 WN=exp(-i*2*pi/Nm); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nmr) for p=1:Nm/2 % 本次跨越间隔内的各次蝶形

40、运算 for k=p:Nm:N % 本次蝶形运算的跨越间隔为Nmr=2mm kp=k+Nm/2; % 确定蝶形运算的对应单元下标 a=Xk(kp)*u; % 蝶形运算的乘积项 Xk(kp)=Xk(k)-a; % 蝶形运算 Xk(k)=Xk(k)+a; % 蝶形运算 end u=u*WN; % 修改旋转因子,多乘一个基本DFT因子WN end end y=Xk;ditfft=abs(y);subplot(211);stem(n,ditfft,'b'); %绘制傅里叶变换xlabel('k');ylabel('X(k)');title('8点

41、DIFFFT计算X(k)')grid on;check=fft(xn,N);fft_abs=abs(check);subplot(212)stem(n,fft_abs,'b'); %进行快速傅里叶验证xlabel('k');ylabel('X(k)');title('对结果进行FFT验证X(k)')grid on;b) 图像结果：图2.8 DIF-基2FFTc) 结论：两者计算结果一致，故程序正确。（2）基于（1）中的两个FFT函数，分别通过两个方法实现IFFT操作，并通过与原x(n)的对比验证IFFT操作的正确性。i.

42、旋转因子指数变极性法a) 代码：clc;close all;clear;format compact; Xk=1 3 5 9 6 7 0 3; %初始序列XkM=nextpow2(length(Xk), N=2M, for m=0:N/2-1; %旋转因子指数范围 WN(m+1)=exp(j*2*pi/N)m; %计算旋转因子(改变极性) end A=Xk,zeros(1,N-length(Xk); %数据输入，长度不够补零disp('输入到各存储单元的数据:'),disp(A); J=N/2; %给倒序数赋初值 for I=1:N/2-1; %按序交换数据和算倒序数 if I

43、<J; %条件判断及数据交换 T=A(I+1);A(I+1)=A(J+1);A(J+1)=T; end K=N/2; while J>=K; J=J-K;K=K/2; end J=J+K; end %disp('倒序后各存储单元的数据:'),disp(A); for L=1:M; %分级计算 % disp('运算级次:'),disp(L); B=2(L-1); for R=0:B-1; %各级按序蝶算 P=2(M-L)*R; for K=R:2L:N-2; %每序依次计算 T=A(K+1)+A(K+B+1)*WN(P+1); A(K+B+1)=A(K

44、+1)-A(K+B+1)*WN(P+1); A(K+1)=T; end end %disp('本级运算后各存储单元的数据:'),disp(A); end %disp('输出各存储单元的数据:')xn1=A*1/N; %结果乘以1/Nsubplot(211)stem(xn1); title('旋转因子指数变极性法实现IFFT计算xn')%disp('调用fft函数运算的结果:'),xn2=ifft(Xk,N);subplot(212)stem(xn2);title('对结果进行IFFT验证xn')b) 图像结果：图2

45、.9 旋转因子指数变极性法c) 结论：通过图像可以看出程序编写无误。ii. 直接调用FFT子程序法：a) 代码：clc;close all;clear;format compact; Xk=1 3 5 9 6 7 0 3; %初始序列XkXk1=conj(Xk); %取初始序列的共轭%调用FFT子程序计算M=nextpow2(length(Xk1), N=2M, for m=0:N/2-1; %旋转因子指数范围 WN(m+1)=exp(-j*2*pi/N)m; %计算旋转因子 end A=Xk1,zeros(1,N-length(Xk1); %数据输入，长度不够补零disp('输入到各

46、存储单元的数据:'),disp(A); J=N/2; %给倒序数赋初值 for I=1:N/2-1; %按序交换数据和算倒序数 if I<J; %条件判断及数据交换 T=A(I+1);A(I+1)=A(J+1);A(J+1)=T; end K=N/2; while J>=K; J=J-K;K=K/2; end J=J+K; end %disp('倒序后各存储单元的数据:'),disp(A); for L=1:M; %分级计算 % disp('运算级次:'),disp(L); B=2(L-1); for R=0:B-1; %各级按序蝶算 P=2

47、(M-L)*R; for K=R:2L:N-2; %每序依次计算 T=A(K+1)+A(K+B+1)*WN(P+1); A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1); A(K+1)=T; end end %disp('本级运算后各存储单元的数据:'),disp(A); end %disp('输出各存储单元的数据:')xn1=conj(A)*1/N; %共轭并乘以1/Nsubplot(211)stem(xn1); title('直接调用FFT子程序法实现IFFT计算xn')%disp('调用fft函数运算的结果:'

48、;),xn2=ifft(Xk,N);subplot(212)stem(xn2);title('对结果进行IFFT验证xn')b) 图像结果： 图2.10 直接调用FFT子程序法c) 结论：由MATLAB图像对比可知，程序编写无误。实验三实验目的： 加深对FIR滤波器和IIR滤波器的认识与理解。实验原理与方法 数字滤波器基于单位脉冲响应的特性分为无限长单位脉冲响应数字滤波器（IIR）和有限长单位脉冲响应数字滤波器（FIR）。 IIR滤波器可以利用成熟的模拟滤波器进行设计，但是是非线性相位，常用方法有脉冲响应不变法和双极性不变法；FIR滤波器可严格线性相位并能够任意幅度特性，且为因

49、果稳定系统，可用FFT计算，但是阶次比IIR高的多，常用方法窗函数法。实验内容及步骤1、利用设计IIR低通、带通、高通三个滤波器，对音频信号中的低频部分、中频部分、高频部分进行提取；2、利用设计FIR低通、带通、高通三个滤波器，对音频信号中的低频部分、中频部分、高频部分进行提取；3、对FIR滤波器和IIR滤波器的处理效果进行对比，包括：感性对比和理性分析。4、试将音频中所提取的各部分的能量进行改变，再将各部分重组，比较一下重组后的音频与原始音频的变化。实验结果分析及结论总结1. 利用设计IIR低通、带通、高通三个滤波器，对音频信号中的低频部分、中频部分、高频部分进行提取；i. 代码：a) 主函

50、数：y,Fs=audioread('E晓霞m文件实验3大王叫我来巡山.wav');%读取音频，获得采样频率N=length(y); %采样点数T=N/Fs; %采样时间turn=(0:N-1)/Fs; %定义时间轴furn=(-N/2:N/2-1)/N*Fs; %定义频率轴% 未处理前信号波形figure(1)subplot(2,1,1)plot(turn,y(:,1) %画出时域波形图xlabel('t'),ylabel('s')title('音频信号时域波形图')yf=fftshift(fft(y(:

51、,1); %对信号做快速傅里叶变换subplot(2,1,2)plot(furn,abs(yf) %画出频域波形图xlabel('f'),ylabel('s')title('音频信号频域波形图')ylow=IIR_LP(y,Fs,furn); %低频部分ybp=IIR_BP(y,Fs,furn); %中频部分yhigh=IIR_HP(y,Fs,furn); %高频部分b) 低通：function ylow=IIR_LP(y,Fs,furn)%椭圆型低通滤波器设计Wp=2*1000/Fs; %通带截止频率Ws=2*1110/Fs; %阻带截止频率Rp=3; %通带内最大衰减系数Rs=60; %阻带内最大衰减系数No,Wc=ellipord(Wp,Ws,Rp,Rs); %设计滤波器的阶数以及3dB截止频率Bz,Az=ellip(No,Rp,Rs,Wc,'low'); %计算滤波器系统函数分子分母多项式 ylow=filter(Bz,Az,y); % 进行低通滤波ylowf1=fftshift(fft(ylow(:,1); % 对信号做FFT变换 figure(2)plot(furn,abs(ylowf1)xlabel('f'),ylabel('s')title('信号低