MATLAB实现FFT 及信号的谱分析

2022-07-20 14:23:36 浏览数 (2)

一、实验目的 １.通过实验加深对 FFT 的理解，熟悉 FFT 程序、结构及编程方法。

２.熟练应用 FFT 对典型信号进行谱分析的方法。

３.了解应用 FFT 进行信号频域分析可能出现的问题以便在实际中正确应用FFT。

 4. 理解 FFT 与 IFFT 的关系。

 5.. 熟悉应用 FFT 实现两个序列的线性卷积的方法。

二、实验原理及方法         在各种信号序列中，有限长序列信号处理占有很重要的位置，对有限长序列，我们可以使用离散 Fourier 变换（DFT）。这一变换不但可以很好地反映序列的频谱特性，而且已于永快速算法在计算机上实现，当序列 x(n) 的长度为 N 时，它的 DFT 定义为：

        有限长序列的 DFT 是其 Z 变换在单位圆上的等距采样，或者说是序列 Fourier 变换的等距采样，因此可以用于序列的谱分析。FFT 并不是与 DFT 不同的另一种变幻，而是为了减少 DFT运算次数的一种快速算法。它是对变换式进行一次次的分解，是其成为若干小点数的组合，从而减少运算量。常用的 FFT 是以 2 为基数的，其长度 N = 2L 。它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于 2 的整数次方时，为了使用以２为基数的 FFT，可以用末位补零的方法，是其长度延长至 2 的整数次方。

（一）、在利用 DFT 进行频谱分析时可能会出现三种误差。 (1) 混叠         为了利用计算一个连续信号的频谱，首先需要对这个连续信号进行取样，如果取样频率太低，也即抽样周期太大，在频域内将产生混叠现象，这样就不可能无失真的恢复原连续信 号。对带限信号，当所处理模拟信号最高频率 fh 与抽样频率 fs 满足fs  ≥ 2 fh时就不会出现频谱混叠现象。         然而，时域内有限长的信号，其频谱宽度是无限的，为了使有限长信号满足抽样定理，在进行抽样之前，可以先用低通模拟滤波器对信号进行滤波,从而保证高于折叠频率的分量不会出现。 （2）泄漏         实际信号序列往往很长，甚至是无限长序列。为了方便，我们往往用截短的序列来近似他们。这样可以使用较短的 DFT 来对信号进行频谱分析。对序列 x(n) 截短的过程就是将 原信号序列与矩形窗函数相乘的过程，在频域就是两者频谱的卷积。一般情况下这样都会造成由此得到的频谱不同于信号原来的频谱，这种现象叫做泄漏。在实际应用中，可以选用频 谱主瓣小、旁瓣小、尽量接近于(Ω) 的窗函数来减少泄漏。         泄漏不能与混叠完全分开，因为泄漏导致频谱的扩展，从而造成混叠。为了减少泄漏的影响，可以选择适当的窗函数是频谱的扩散减至最小。 （3）栅栏效应         DFT 是对单位圆上 z 变换的均匀取样，所以它不可能将频谱视为一个连续函数。这样就产生了栅栏效应。就一定的意义上看，用 DFT 来观看频谱就好像通过一个尖桩的栅栏来观看一个图景一样，只能在离散点上看到真实的频谱。这样就有可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所挡住，不能被我们观察到。减小栅栏效应的一个方法是借助于原序列的末端增加一些零值，从而变动 DFT 的点数。这一方法实际上是人为的改变了对真实频谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得原来看不到的频谱的峰点或谷点就有可能看到了。         IFFT 一般可以通过 FFT 程序来完成，只要对 X[k]取共轭，进行 FFT 运算，然后再取共轭，并乘以因子 1/N，就可以完成 IFFT。         实验中用到的信号序列：

 三、实验内容 1．用三种不同的 DFT 程序计算 x(n) = R₈ (n) 的傅里叶变换 运行时间。 （1）用 for loop  语句的M 函数文件dft1.m，用循环变量逐点计算 X (k ) ； （2）编写用 MATLAB 矩阵运算的 M 函数文件 dft2.m,  完成下列矩阵运算； （3）调用 FFT 库函数，直接计算 X (k ) ； （4）分别利用上述三种不同方式编写的 DFT 程序计算序列 x(n) 的傅立叶变换 X (e ʲʷ ) ，并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。 M 函数文件如下：

代码语言：javascript复制
dft1.m:
function[Am,pha]=dft1(x)
N=length(x);
w=exp(-j*2*pi/N);
for k=1:N
sum=0;
for n=1:N
sum=sum x(n)*w^((k-1)*(n-1));
end
Am(k)=abs(sum);
pha(k)=angle(sum);
end
dft2.m:
function[Am,pha]=dft2(x)
N=length(x);
n=[0:N-1];
k=[0:N-1];
w=exp(-j*2*pi/N);
nk=n'*k;
wnk=w.^(nk);
Xk=x*wnk;
Am=abs(Xk);
pha=angle(Xk);

四、思考题 FFT   在什么条件下也可以用来分析周期信号序列的频谱？ 如果正弦信号sin(2Πf₀k ), f₀ = 0.1Hz ，用 16 点来做 FFT 运算，得到的频谱是信号本身的真实谱吗？为什么？

五、实验报告要求 1．简述实验原理及目的。 2．给出所编制的实验主程序、实验信号序列的时域和频域图形并分析所得图形，说明参数 改变时对时域和频域信号波形的影响。 3．简要回答思考题。

代码语言：javascript复制
x=ones(1,8);
n=0:length(x)-1;
t0=clock;
[m1,a1]=dft1(x);
subplot(321),stem(n,m1,'.');
title('用dft1实现 mag');
subplot(322),stem(n,a1,'.');
title('angle');
t1=etime(clock,t0);
t0=clock;
[m2,a2]=dft2(x);
subplot(323),stem(n,m2,'.');
title('用dft2实现 mag');
subplot(324),stem(n,a2,'.');
title('angle');
t2=etime(clock,t0);
t0=clock;
xk=fft(x);
subplot(325),stem(n,abs(xk),'.');
title('用FFT实现 mag');
subplot(326),stem(n,angle(xk),'.');
title('angle');
t3=etime(clock,t0);
figure;
subplot(311),stem(t1,'.');ylabel('单位：s');title('用dft1实现');
subplot(312),stem(t2,'.');ylabel('单位：s');title('用dft2实现');
subplot(313),stem(t3,'.');ylabel('单位：s');title('用FFT实现');
调用的函数：
�t1.m
function[Am,pha]=dft1(x)
N=length(x);
w=exp(-j*2*pi/N);
for k=1:N
sum=0;
for n=1:N
sum=sum x(n)*w^((k-1)*(n-1));
end
Am(k)=abs(sum);
pha(k)=angle(sum);
end
�t2.m
function[a,p]=dft2(x)
N=length(x);
n=0:N-1;
k=n;
w=exp(-j*2*pi/N);
nk=n'*k;
wnk=w.^(nk);
xk=x*wnk;
a=abs(xk);
p=angle(xk);
代码语言：javascript复制
p=8;
n=0:15;
figure;
for k=1:3
q=2^k;
x1=exp(-((n-p).^2)/q);
subplot(3,2,2*k-1),stem(n,x1,'.');
title('q变化 x(n)');
[m,a,w]=dtft(x1);
subplot(3,2,2*k),plot(w/pi,m);
title('q变化 abs[X(w)]');
end
q=8;
figure;
for k=1:3
p=-2*(k^2) 11*k-1;
x1=exp(-((n-p).^2)/q);
subplot(3,2,2*k-1),stem(n,x1,'.');
title('p变化 x(n)');
[m,a,w]=dtft(x1);
subplot(3,2,2*k),plot(w/pi,m);
title('p变化 abs[X(w)]');
end
figure;
f1=[0.0625,0.4375,0.5625];
for k=1:3
f=f1(k);
a=0.1;
x2=exp(-a)*sin(2*pi*f*n);
subplot(3,2,2*k-1),stem(n,x2,'.');
title('f 变化 x(n)');
[m,a,w]=dtft(x2);
subplot(3,2,2*k),plot(w/pi,m);
title('f 变化 abs[X(w)]');
end
调用的函数：
%dtft.m
function[m,a,w]=dtft(x)
N=length(x);
n=0:N-1;
w=linspace(-2*pi,2*pi,500);
y=x*exp(-j*n'*w);
m=abs(y);
a=angle(y);
代码语言：javascript复制
3.程序：
f=64;
dt=1/f;
N=[16,32,64];
for k=1:3
n=0:N(k)-1;
nt=n*dt;
x=cos(8*pi*nt) cos(16*pi*nt) cos(20*pi*nt);
y=fft(x);
subplot(3,3,3*k-2),stem(n,x,'.');
title('x(n)');
subplot(3,3,3*k-1),stem(n,abs(y),'.');
title('abs(X(k))');
subplot(3,3,3*k),stem(n,angle(y),'.');
title('angle(X(k))');
end
代码语言：javascript复制
n=0:127;
x=cos(n*pi/36) cos(1.5*n*pi/36);
y=fft(x);
figure;
subplot(311),stem(n,x,'.');
title('x(n)');
subplot(312),stem(n,abs(y),'.');
title('X(k) 幅度谱');
subplot(313),stem(n,angle(y),'.');
title('X(k) 相位谱');
N=length(n);
n1=0:N/4-1;
x1=cos(4*n1*pi/36) cos(1.5*4*n1*pi/36);
x1=[x1,zeros(1,N*3/4)];
y1=fft(x1,128);
figure;
subplot(311),stem(n,x1,'.');
title('x(4n)');
subplot(312),stem(n,abs(y1),'.');
title('对应的X(k) 幅度谱');
subplot(313),stem(n,angle(y1),'.');
title('对应的X(k) 相位谱');
n2=0:(N-1)*4;
x2=zeros(1,N*4-3);
for k=0:127
x2(k*4 1)=cos(k*pi/36) cos(1.5*k*pi/36);
end
x3=zeros(1,N);
x3=x2(1:N);
y2=fft(x3,128);
figure;
subplot(311),stem(n,x3,'.');
title('x(n/4)');
subplot(312),stem(n,abs(y2),'.');
title('对应的X(k) 幅度谱');
subplot(313),stem(n,angle(y2),'.');
title('对应的X(k) 相位谱');
fft fs

0 人点赞