《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 页