频谱分析实验报告.docx

基于 Matlab 的频谱分析一、实验目的1．频谱分析的基本概念及应用；2．常用信号的频谱：余弦信号、方波信号、三角波信号3 .应用DFT进行频谱分析，MATLAB编程仿真(1) 整周期采样、非整周期采样一频谱泄漏；(2) 填零运算一高密度频谱；( 3)非整周期采样时，加窗改进；4. 应用频谱分析仪测量常用信号的频谱二、实验原理1、利用 DFT 计算连续信号的频谱分析过程X a( j°) X (小X (ej 3)NX (e j 3 )N 3 = 2 兀 k / N从图中看出，这是一次近似过程，首先，用离散采样信号的 DTFTX (ej3)来近似连续信号x (t)的傅立叶变换X (j° )，其次，将x(n)截短，这一过程 aa等效于用一矩形序列RN (n)与x(n)相乘，其DTFT为X ( e j3 ) = X ( e j3 ) * R (e j3 )NN最后，再对截短的信号作 DFT2、频谱分析的近似过程出现问题及解决办法( 1 )混叠①混叠现象对连续信号X(t)进行数字处理前，要进行采样x (nT ) = x(t) Y 6 (t 一 n T )an = —g采样序列的频谱是连续信号频谱的周期延拓，周期为fs,如采样率过低，不满足 采样定理，fsv2fh，则导致频谱混迭，使一个周期内的谱对原信号谱产生失真， 无法恢复原信号，进一步的数字处理失去依据。

另外，干扰信号的频率高于折叠 频率 f / 2 的也会产生频谱混叠s②对于干扰信号的频率高于 f /2产生频谱混叠的解决办法：s在采样前加一抗混叠滤波器，滤除高于折叠频率 f / 2 的频率成分s(2) 泄漏① 产生原因处理实际信号序列x (n)时，一般总要将它截断为一有限长序列，长为N 点，相当于乘以一个矩形窗w(n)=RN(n)矩形窗函数，其频谱有主瓣，也有许 多副瓣，窗口越大，主瓣越窄，当窗口趋于无穷大时，就是一个冲击函数我们知道，时域的乘积对应频域的卷积，所以，加窗后的频谱实际是原信号 频谱与矩形窗函数频谱的卷积，卷积的结果使频谱延伸到了主瓣以外，且一直延 伸到无穷当窗口无穷大时，与冲击函数的卷积才是其本身，这时无畸变，否则 就有畸变例如，信号为e^『T，是一单线谱，但当加窗后，线谱与抽样函数进行卷 积，原来在Q0处的一根谱线变成了以Q0为中心的，形状为抽样函数的谱线序 列，原来在一个周期(Qs)内只有一个频率上有非零值，而现在一个周期内几 乎所有频率上都有非零值，即x (eQ丁)的频率成份从Q 0处“泄漏”到其它频率 处去了考虑各采样频率周期间频谱“泄漏”后的互相串漏，卷积后还有频谱混迭现 象产生。

② 解决办法：增加截短的长度N，以减小泄漏 3 )栅栏效应①栅栏效应的产生N点DFT是在频率区间［0，2n ］上对信号频谱进行N点等间隔采样，得 到的是若干个离散的频谱点X (k)，且它们限制在基频的整数倍上，这就好像 在栅栏的一边通过缝隙看另一边的景象一样，只能在离散点处看到真实的景象， 其余部分频谱成分被遮挡， 所以称之为栅栏效应②减小栅栏效应方法：尾部补零，使谱线变密，增加频域采样点数，原来漏 掉的某些频谱分量就可能被检测出来4) DFT 的分辨率填补零值可以改变对DTFT的采样密度，人们常常有一种误解，认为补零可 以提高DFT的频率分辨率事实上我们通常规定DFT的频率分辨率为f / n， s 这里的N是指信号x(n)的有效长度，而不是补零的长度不同长度的x(n)其DTFT 的结果是不同的；而相同长度的x(n)尽管补零的长度不同其DTFT的结果应是相 同的，他们的DFT只是反映了对相同的DTFT采用了不同的采样密度参数选择的一般原则：(2)根据频率分辩率af(1)若已知信号的最高频率f ，为防止混叠，选定采样频率f > 2 f ; m ax s m ax，确定所需DFT的长度n = f / af ； s(3)f和N确定以后，即可确定相应模拟信号的时间长度T = N / f = NT (这ss里 T 是采样周期)。

三、实验内容及结果1、余弦波频谱分析对余弦函数x(t)=cos(2兀*50t)进行频谱分析，采样频率为10000Hz，对其进 行整周期采样，非整周期采样，频谱泄露及其改善方法其程序代码和运行结果 如下程序代码：fs = 10000;t = 0:1/fs:1.5;%xx = cos(2*pi*50*t);subplot(321),plot(t,x), axis([0 0.2 -1.2 1.2]); title('x(t)=cos(2*pi*50*t)');grid on;% x(n) nx=[0:99];x0=cos(2*pi*1/50*nx); subplot(322),stem(nx,x0);title('x(n)');grid on;% 整周期采样N1=100;x1=x0(1:N1);X1=fft(x1,N1);magX1=abs(X1(1:N1/2+1))*2/N1; w1=[0:N1/2]*2*pi/N1;subplot(323);stem(wl/pi,magXl); title('整周期采样');grid on;% 265点DFT, high resolution spectrumN2=256;X=fft(x0,N2);magX=abs(X(1:N2/2+1))*2/N2;w=[0:N2/2]*2*pi/N2;subplot(324);stem(w/pi,magX) ;title('256点DFT, high resolution spectrum');grid on;%加长矩形窗长度mx=[0:149];xm=cos(2*pi*1/50*mx);Nm=256;Y=fft(xm,Nm);magY=abs(Y(1:Nm/2+1))*2/Nm;wm=[0:Nm/2]*2*pi/Nm;subplot(325);stem(wm/pi,magY);title('矩形窗加长到 150' );grid on;%改窗N=256;wn = hamming(length(x0));hn = x0.*wn';H = fft(hn,N);magH=abs(H(1:N/2+1))*2/N;wH=[0:N/2]*2*pi/N;subplot(326);stem(wH/pi,magH) ;title('力加hamming窗')；grid on;运行结果：2、方波信号频谱对方波函数x(t)= square (2 *50t)进行频谱分析，采样频率为10000Hz，对 其进行整周期采样，非整周期采样，频谱泄露及其改善方法。

其程序代码和运行 结果如下程序代码：fs = 10000;t = 0:1/fs:1.5;%xx = square(2*pi*50*t);subplot(321),plot(t,x), axis([0 0.2 -1.2 1.2]); title('x(t)=square(2*pi*50*t)');grid on;% x(n) nx=[0:99];x0=square(2*pi*1/50*nx); subplot(322),stem(nx,x0);title('x(n)');grid on;% 整周期采样N1=100;x1=x0(1:N1);X1=fft(x1,N1);magX1=abs(X1(1:N1/2+1))*2/N1; w1=[0:N1/2]*2*pi/N1;subplot(323);stem(w1/pi,magX1);title('整周期采样')；grid% 265点DFT, high resolution spectrumN2=256;X=fft(x0,N2);magX=abs(X(1:N2/2+1))*2/N2;w=[0:N2/2]*2*pi/N2;subplot(324);stem(w/pi,magX);title('256点DFT, highresolution spectrum');grid on;% 加长矩形窗长度mx=[0:149];xm=square(2*pi*1/50*mx);Nm=256;Y=fft(xm,Nm);magY=abs(Y(1:Nm/2+1))*2/Nm;wm=[0:Nm/2]*2*pi/Nm;subplot(325);stem(wm/pi,magY);title('矩形窗加长到 150' );grid on;% 改窗N=256;wn = hamming(length(x0));hn = x0.*wn';H = fft(hn,N);magH=abs(H(1:N/2+1))*2/N;wH=[0:N/2]*2*pi/N;subplot(32 6);stem(wH/pi,magH);title('力口hamming窗');gridon;on;运行结果：3、锯齿波信号频谱对方波函数X(t)= sawtooth (2 *50t)进行频谱分析，采样频率为10000Hz, 对其进行整周期采样，非整周期采样，频谱泄露及其改善方法。

其程序代码和运 行结果如下程序代码：fs = 10000;t = 0:1/fs:1.5;%xx = sawtooth(2*pi*50*t);subplot(321),plot(t,x), axis([0 0.2 -1.2 1.2]); title('x(t)=sawtooth(2*pi*50*t)');grid on;% x(n) nx=[0:99];x0=sawtooth(2*pi*1/50*nx); subplot(322),stem(nx,x0);title('x(n)');grid on;% 整周期采样N1=100;x1=x0(1:N1);X1=fft(x1,N1);magX1=abs(X1(1:N1/2+1))*2/N1; w1=[0:N1/2]*2*pi/N1;subplot(323);stem(w1/pi,magX1);title('整周期采样')；grid% 265点DFT, high resolution spectrumN2=256;X=fft(x0,N2);magX=abs(X(1:N2/2+1))*2/N2;w=[0:N2/2]*2*pi/N2;subplot(324);stem(w/pi,magX);title('256点DFT, highresolution spectrum');grid on;% 加长矩形窗长度mx=[0:149];xm=sawtooth(2*pi*1/50*mx);Nm=256;Y=fft(xm,Nm);magY=abs(Y(1:Nm/2+1))*2/Nm;wm=[0:Nm/2]*2*pi/Nm;subplot(325);stem(wm/pi,magY);title('矩形窗加长到 150' );grid on;N=256; 改%窗wn = hamming(length(x0));hn = x0.*wn';H = fft(hn,N);magH=abs(H(1:N/2+1))*2/N;wH=[0:N/2]*2*pi/N;on;on;subplot(32 6);stem(wH/pi,magH);title('力口hamming窗');grid 运行结果：对方波函数 x(t)= 。