《2022年数字信号处理zang实验》由会员分享，可在线阅读，更多相关《2022年数字信号处理zang实验（10页珍藏版）》请在金锄头文库上搜索。
1、实验名称：数字信号处理实验实验项目：用 FFT作谱分析 指导老师：李秋菊班级 ： 10电科姓名：张卫娟学号： 201000804084 成绩：一、实验目的：1、在理论学习的基础上， 通过本实验， 加深对 FFT的理解， 熟悉 MATLAB 中的有关函数。2、熟悉应用 FFT对典型信号进行频谱分析的方法。熟悉FFT算法原理和 FFT子程序的应用。3、学习用 FFT对连续信号和时域离散信号进行谱分析的方法。了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT 。二、实验原理：（一）、在各种信号序列中，有限长序列信号处理占有很重要地位，对有限长序列，我们可以使用离散傅里叶变
2、换(DFT)。这一变换不但可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x(n) 的长度为 N 时，它的 DFT定义为：反变换为：有限长序列的DFT是其 Z 变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。在信号处理中， DFT 的计算具有举足轻重的地位， ，信号的相关、滤波、谱估计等都要通过 DFT 来实现。然而，当 N 很大的时候， 求一个 N 点的 DFT 要完成NN次复 数乘 法和) 1(NN次复 数加 法， 其计 算量 相当 大 。 1965 年J.W.Cooley和 J.W.Tukey巧妙地利用NW因子的周期性和
3、对称性， 构造了一个 DFT快速算法，即快速傅立叶变换(FFT)。 FFT 并不是与 DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小点数的组合， 从而减少精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 10 页运算量。常用的 FFT是以 2 为基数的， 其长度。它的效率高， 程序简单，使用非常方便， 当要变换的序列长度不等于2 的整数次方时， 为了使用以 2 为基数的 FFT ，可以用末位补零的方法，使其长度延长至2 的整数次方。（二）、在运用 DFT进行频谱分析的过程中可能
4、的产生混叠误差序列的频谱是被采样信号频谱的周期延拓，当采样速率不满足Nyquist定理时，就会发生频谱混叠， 使得采样后的信号序列频谱不能真实的反映原信号的频谱。避免混叠现象的唯一方法是保证采样速率足够高，使频谱混叠现象不致出现，即在确定采样频率之前，必须对频谱的性质有所了解，在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前，先用低通模拟滤波器对信号进行滤波。（三）、 matlab 函数应用：MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有 Fft、Ifft 、Fft2 、Ifft2, Fftn、ifftn 和 Fftshift、Ifftshift 等。
5、当所处理的数据的长度为 2 的幂次时，采用基 -2 算法进行计算， 计算速度会显著增加。 所以， 要尽可能使所要处理的数据长度为2 的幂次或者用添零的方式来添补数据使之成为 2 的幂次。1、fft 和 Ifft 函数（1）调用方式： Y=fft(X) 参数说明如果 X 是向量，则采用傅立叶变换来求解X 的离散傅立叶变换；如果 X 是矩阵，则计算该矩阵每一列的离散傅立叶变换；如果 X 是（N)D维数组，则是对第一个非单元素的维进行离散傅立叶变换；（2） Yfft(X,N) 参数说明精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 10 页N
6、 是进行离散傅立叶变换的X 的数据长度， 可以通过对 X 进行补零或截取来实现。（3）Yfft(X,dim) 或 Yfft(X,N,dim) 参数说明在参数 dim 指定的维上进行离散傅立叶变换；当 X 为矩阵时，dim 用来指定变换的实施方向： dim=1,表明变换按列进行；dim=2 表明变换按行进行。函数 Ifft 的参数应用与函数Fft 完全相同。应用说明【实例 1】fft 的应用X=2 1 2 8; Y=fft(X ，4) 运行结果Y13.0000 0+7.0000i -5.0000 0-7.0000 【实例 2】fft(X,N,dim)的应用A=2 5 7 8; 1 4 0 5;
7、3 8 5 1; 9 1 2 7; Z=fft(A,1) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 10 页2、 Fftshift 和 Ifftshift函数调用方式Z=fftshift(Y) 此函数可用于将傅立叶变换结果Y（频域数据） 中的直流成分 （即频率为 0处得值）移到频谱的中间位置。【实例 3】 fftshift 的应用X=rand(5,4); y=fft(X); z=fftshift(y);% 只将傅立叶变换结果y 中的直流成分移到频谱的中间位置. 运行结果：y= 3.2250 2.5277 1.4820 1.6314
8、 0.3294+0.2368i 0.0768+0.3092i 0.6453+0.4519i -0.7240-0.4116i -0.2867-0.6435i 0.5657+0.4661i -0.5515+0.2297i -0.0573-0.0881i -0.2867+0.6435i 0.5657-0.4661i -0.5515-0.2297i -0.0573+0.0881i 0.3294-0.2368i 0.0768-0.3092i 0.6453-0.4519i -0.7240+0.4116i Z= -0.5515-0.2297i -0.0573+0.0881i -0.2867+0.6435i
9、0.5657-0.4661i 0.6453-0.4519i -0.7240+0.4116i 0.3294-0.2368i 0.0768-0.3092i 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 10 页1.4820 1.6314 3.2250 2.5277 0.6453+0.4519i -0.7240-0.4116i 0.3294+0.2368i 0.0768+0.3092i -0.5515+0.2297i -0.0573-0.0881i -0.2867-0.6435i 0.5657+0.4661i 【实例 4】fft 在信号分析
10、中的应用使用频率分析方法对模拟信号x(t) x=sin(2*pi*100*t) 进行频谱分析。采样频率fs=1000，采样点数 N=512。并画出信号的时域波形及FFT 变换后的幅频响应与相频响应。程序：fs=512;% 采样频率N=512;% 数据点数n=0:N-1; t=0:1/fs:(N-1)/fs;% 采样时间序列f0=100;% 信号频率x=sin(2*pi*f0*t); subplot(3,1,1); plot(t,x); xlabel(t); ylabel(sin(2*pi*100*t); title( 时域信号 ); Y=fft(x,N);% 对信号进行FFT 变换magY=a
11、bs(Y);% 求得 FFT 变换后的幅度angY=angle(Y)*180/pi;% 求得 FFT 变换后的相位f=n*fs/N;% 频率序列subplot(3,1,2); plot(f,magY);% 画出幅频响应xlabel(f); ylabel(幅度 ); title(N=512); grid; subplot(3,1,3); plot(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f); ylabel(相位 ) title(N=512); grid; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 10
12、 页00.10.20.30.40.50.60.7-101tsin(2*pi*100*t)时 域 信 号010020030040050060070080090010000200400f幅度N=51201002003004005006007008009001000-2000200f相位N=512结果分析：假设采样频率为Fs，信号频率 F，采样点数为N。那么 FFT 之后结果就是一个为 N 点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。假设原始信号的峰值为A，那么 FFT 的结果的每个点 （除了第一个点直流分量之外）的模值就是 A 的 N/2 倍。（例如 f=100H
13、z 点处，A=1*N/2=1*512/2=256, 而第一个点就是直流分量， 它的模值就是直流分量的N 倍)。magY=abs(Y);% 求得 FFT 变换后的幅度如果要换算为实际的幅度需要改为：magY=magY/(N/2); magY(1)=magY(1)/2; 每个点的相位 ,就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz） ， 而最后一个点 N 的再下一个点（假设的第 N+1 个点）则表示采样频率Fs，这中间被 N-1 个点平均分成 N 等份，每个点的频率依次增加。第n 个点所表示的频率为：Fn=(n-1)*Fs/N。 由上面的公式可以看出， Fn所能分辨到频率为为Fs/N
14、，如果采样频率 Fs 为 1024Hz，采样点数为 1024 点，则可以分辨到1Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。由于 FFT 结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。前面的程序由： plot(f(1:N/2),magY(1:N/2);% 画出幅频响应三、实验内容：1、对被白噪声污染的信号2/150cos5 .16/100cos32Stt进行频谱分析，从中鉴别出有用的信号。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 10 页要求：将信号的幅度换算成实际的幅度，信号的频率换算成
15、实际的频率程序：clear; fs=512;% 采样频率N=512;% 数据点数n=0:N-1; t=0:1/fs:(N-1)/fs;% 采样时间序列f0=100;% 信号频率x=2+3*cos(100*pi*t-pi/6)+1.5*cos(150*pi*t+pi/2)+randn(1,N); subplot(3,1,1); plot(t,x); xlabel(t); ylabel(sin(2*pi*100*t); title( 时域信号 ); Y=fft(x,N);% 对信号进行FFT 变换magY=abs(Y)/N*2;% 求得FFT 变换后的幅度magY(1)=magY(1)/2; an
16、gY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;% 频率序列subplot(3,1,2); plot(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f); ylabel(幅度 ); title(N=512); grid; subplot(3,1,3); plot(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f); ylabel(相位 ) title(N=512); grid; 图形：2、对连续的单一频率周期信号 , t2sinS , 信号频率 f=1Hz 按采样频率 fs=8Hz采样，截取长度N分别选 N=2
17、0和 N=16，观察其 DFT结果的幅度谱。幅频N=20 程序：clear; fs=8;%采样频率N=20;% 数据点数n=0:N-1; t=0:1/fs:(N-1)/fs;% 采样时间序列f0=100;% 信号频率x=sin(2*pi*t); subplot(3,1,1); stem(t,x); xlabel(t); ylabel(sin(2*pi*100*t); title(时域信号 ); Y=fft(x,N);% 对信号进行FFT 变换magY=abs(Y)/(N/2);% 求得 FFT 变换后的幅度magY(1)=magY(1)/2; 精选学习资料 - - - - - - - - -
18、名师归纳总结 - - - - - - -第 7 页，共 10 页angY=angle(Y)*180/pi;% 求得 FFT 变换后的相位f=n*fs/N;% 频率序列subplot(3,1,2); stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f); ylabel(幅度 ); title(N=20); grid; subplot(3,1,3); stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f); ylabel(相位 ) title(N=20); grid; N=16 程序：clear; fs=16;%采样频率N=16;% 数
19、据点数n=0:N-1; x=1,2,3,4,4,3,2,1,zeros(1,8); subplot(3,1,1); stem(x); xlabel(n); ylabel(sin(2*pi*100*t); title( 时域信号 ); Y=fft(x);% 对信号进行FFT 变换magY=abs(Y)/N*2;% 求得FFT 变换后的幅度magY(1)=magY(1)/2; angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;% 频率序列subplot(3,1,2); stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f); yl
20、abel(幅度 ); title(N=16); grid; subplot(3,1,3); stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f); ylabel(相位 ) title(N=16); grid; 图形：3、对三角波序列74, n-830, 1)(nnnnxc，elsenxc,0)(进行频谱分析设 n=16。程序：clear; fs=16;%采样频率N=16;% 数据点数n=0:N-1; x=1,2,3,4,4,3,2,1,zeros(1,8); 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共
21、 10 页subplot(3,1,1); stem(x); xlabel(n); ylabel(sin(2*pi*100*t); title( 时域信号 ); Y=fft(x);% 对信号进行FFT 变换magY=abs(Y)/N*2;% 求得FFT 变换后的幅度magY(1)=magY(1)/2; angY=angle(Y)*180/pi;% 求得 FFT 变换后的相位f=n*fs/N;% 频率序列subplot(3,1,2); stem(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f); ylabel(幅度 ); title(N=16); grid; subplo
22、t(3,1,3); stem(f(1:N/2),angY(1:N/2);%画出相频响应xlabel(f); ylabel(相位 ) title(N=16); grid; 图形：4、对矩形波进行频谱分析。矩形波为脉冲宽度为2s，持续时间为 -55s，采用矩形脉冲函数 rectpuls(t,2); 程序：clear; fs=64;%采样频率N=64;% 数据点数n=0:N-1; t=-5:1/fs:5; x=rectpuls(t,2); subplot(3,1,1); plot(t,x); xlabel(n); ylabel(sin(2*pi*100*t); title( 时域信号 ); Y=ff
23、t(x);% 对信号进行FFT 变换magY=abs(Y)/(N/2);% 求得 FFT 变换后的幅度magY(1)=magY(1)/2; angY=angle(Y)*180/pi;%求得 FFT 变换后的相位f=n*fs/N;% 频率序列subplot(3,1,2); plot(f(1:N/2),magY(1:N/2);%画出幅频响应xlabel(f); ylabel(幅度 ); title(N=64); grid; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 10 页图形：四 实验结论通过此次试验，加深对FFT的理解，在各种信号
24、序列中，对有限长序列，我们可以使用离散傅里叶变换(DFT)。这一变换不但可以很好的反映序列的频谱特性，当序列 x(n) 的长度为 N时，它的 DFT定义为：，反变换为：在信号处理中， 信号的相关、 滤波、谱估计等都要通过DFT 来实现。但是，当 N很大的时候，计算很困难，我们可以巧妙地利用NW因子的周期性和对称性，构造了一个 DFT 快速算法，即快速傅立叶变换(FFT)来计算。当然，在进行频谱分析的过程中， 当采样速率不满足Nyquist 定理时，就会发生频谱混叠， 使得采样后的信号序列频谱不能真实的反映原信号的频谱。为避免混叠现象， 在一般情况下，为了保证高于折叠频率的分量不会出现，在采样前， 先用低通模拟滤波器对信号进行滤波。精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 10 页
