数字信号的 FFT 分析

DSP 实验报告

实验名称：实验二数字信号的 FFT 分析

1．实验要求

数字信号处理的一个重要分支就是信号分析，而信号分析的基本工具是离散傅立叶变换。利用傅立叶变换和级数所形成的频谱分析技术作为处理连续信号的重要工具已经应用得很久了,1956年库力（Cooley ）和图基（Tukey ）所发展的近似频谱的快速算法为频谱分析的数字信号的谱分析铺平了道路。因此，DFT （FFT ）得到广泛应用。本次实验设计了两个内容：

(1) 已知信号

0n N-1()0 n 0, n N

n Q x n ?≤≤=?<≥? 这里，N=25，Q= 0.9+j0.3。可以推导出，

11,011)()()(k k 10nk 10-=--===∑∑-=-=N k QW Q QW W n x k X N N n

N N n N

N n ，首先根据这个式子计算X(k) 的理论值，然后计算输入序列x(n) 的32个值，再利用基2时间抽选的FFT 算法，计算x(n) 的DFT X(k)，与X(k) 的理论值比较（要求计算结果最少6位有效数字）。

(2) 假设信号 x(n) 由下述信号组成：

()0.001*c o s (0.45)s i n (0.3)c o s (0.302)

4x n n n n ππππ=+-- 这个信号有两根主谱线 0.3pi 和 0.302pi 靠的非常近，而另一根谱线 0.45pi 的幅度很小，请选择合适的长度 N 和窗函数，用 DFT 分析其频谱，得到清楚的三根谱线。通过本次实验，应该掌握：

(a) 用傅立叶变换进行信号分析时基本参数的选择。

(b) 经过离散时间傅立叶变换（DTFT ）和有限长度离散傅立叶变换（DFT ）后信号频谱上的区别，前者 DTFT 时间域是离散信号，频率域还是连续的，而 DFT 在两个域中都是离散的。

(d) 获得一个高密度频谱和高分辨率频谱的概念和方法，建立频率分辨率和时间分辨率的概念，为将来进一步进行时频分析（例如小波）的学习和研究打下基础。

(e) 建立 DFT 从整体上可看成是由窄带相邻滤波器组成的滤波器组的概念，此概念的一个典型应用是数字音频压缩中的分析滤波器，例如 DVD AC3 和MPEG Audio 。

2. 程序分析

第一部分：

1.25点DFT变换

Q= 0.9+j*0.3;

N=25;

n=0:1:N-1;

k=0:1:N-1;

WN=exp(-j*2*pi/N);

X=(1-Q.^N)./(1-Q*WN.^k);

X=vpa(X,6)

magX=abs(X(1:1:25)); %abs为幅值

k=0:1:24;

w=2*pi*k/25;

stem(w/pi,magX); %stem函数作图

title('25点DFT');

xlabel('2π内幅度特性');

得到X(k)的理论值为：

X =[ 1.83992 + 2.88851*i, 10.0793 + 7.6341*i, 0.751948 - 5.86055*i, 0.231478 - 2.56282*i, 0.28018 - 1.62973*i, 0.333453 - 1.18569*i, 0.374736 - 0.920169*i, 0.407012 - 0.73897*i, 0.433275 - 0.603825*i, 0.455542 - 0.496182*i, 0.475155 - 0.405862*i, 0.493047 - 0.326701*i, 0.50991 - 0.254627*i, 0.526307 - 0.18669*i, 0.542748 - 0.120528*i, 0.559746 - 0.0540151*i, 0.577892 + 0.0150219*i, 0.597934 + 0.0890957*i, 0.620917 + 0.171456*i, 0.64843 + 0.266736*i,

0.683099 + 0.382148*i, 0.729676 + 0.529983*i, 0.797838 + 0.733482*i, 0.910733 + 1.04306*i,

1.13973 + 1.59277*i]

2.计算输入序列x(n)的32个值

n=0:1:N-1;

x=[Q.^n,zeros(1,32-N)]; %zeros为零矩阵

x=vpa(x,6)

序列的后几位数添为0，得到结果为

x =[ 1.0, 0.9 + 0.3*i, 0.72 + 0.54*i, 0.486 + 0.702*i, 0.2268 + 0.7776*i, - 0.02916 + 0.76788*i, - 0.256608 + 0.682344*i, - 0.43565 + 0.537127*i, - 0.553224 + 0.352719*i, - 0.603717 + 0.15148*i, - 0.588789 -

0.0447828*i, - 0.516476 - 0.216941*i, - 0.399746 - 0.35019*i, - 0.254714 - 0.435095*i, - 0.0987144 -

0.467999*i, 0.0515569 - 0.450814*i, 0.181645 - 0.390265*i, 0.28056 - 0.296745*i, 0.341528 - 0.182903*i, 0.362246 - 0.0621539*i, 0.344667 + 0.0527352*i, 0.29438 + 0.150862*i, 0.219684 + 0.22409*i, 0.130488 + 0.267586*i, 0.0371637 + 0.279974*i, 0, 0, 0, 0, 0, 0, 0]

3.利用基2时间抽选的FFT算法，计算x(n) 的DFT X(k)

代码为n=0:1:N-1;

x=[Q.^n,zeros(1,32-N)]

X1=fft(x);

X1=vpa(X1,6)

得到结果为

X1 =[ 1.83992 + 2.88851*i, 4.26186 + 8.57514*i, 10.0249 - 7.66712*i, 1.28721 - 2.56899*i, 0.1705 - 1.94185*i, 0.46568 - 1.84223*i, 0.865851 - 1.25033*i, 0.70937 - 0.63457*i, 0.361783 - 0.53736*i, 0.377467 - 0.707116*i, 0.637133 - 0.615516*i, 0.669139 - 0.292187*i, 0.449665 - 0.162353*i, 0.371812 - 0.295177*i, 0.549183 -

0.329189*i, 0.655053 - 0.113394*i, 0.508892 + 0.0581363*i, 0.375767 - 0.022061*i, 0.490126 - 0.117904*i, 0.647145 + 0.0302847*i, 0.564986 + 0.249659*i, 0.385201 + 0.237077*i, 0.434079 + 0.100955*i, 0.64188 + 0.191054*i, 0.638627 + 0.481008*i, 0.406323 + 0.58139*i, 0.360727 + 0.423857*i, 0.641415 + 0.453716*i, 0.790304 + 0.903673*i, 0.484692 + 1.3182*i, 0.219607 + 1.27159*i, 0.713671 + 1.3331*i]

3.32点DFT变换

n=0:1:N-1;

x=[Q.^n,zeros(1,32-N)]

X1=fft(x);

X1=vpa(X1,6)

magX1=abs(X1(1:1:32));

k=0:1:31;

w=2*pi*k/32;

stem(w/pi,magX1);

title('32点DFT');

xlabel('2π内幅度特性');

3、运行结果

第二部分：

1、DFT

function [Xk] = dft(xn,N)

% Computes Discrete Fourier Transform Coefficients % [Xk] = dft(xn,N)

% Xk = DFT coeff. array over 0 <= k <= N-1

% xn = input signal

% N = length of DFT

n = [0:1:N-1];

k = [0:1:N-1];

WN = exp(-j*2*pi/N);

nk = n'*k;

WNnk = WN .^ nk;

Xk = xn * WNnk;

2.主函数

n=0:1:999;

x=0.001*cos(0.45*n*pi)+sin(0.3*n*pi)-cos(0.302*n*pi-0.25*pi);

stem(n,x);

title('x(n)序列图');

xlabel('n(0<=n<=999)');

X=dft(x,1000); % 计算1000点DFT

magX=abs(X(1:1:501)); % 周期序列镜像对称

k=0:1:500;

w=2*pi*k/1000;

stem(w/pi,magX);

title('π内幅频特性');

xlabel('1000点DFT');

axis([0.29,0.31,0,500]); %axis函数为取横纵坐标的一段长度

xlabel('frequency between 0.29pi and 0.31pi'); %取x坐标的0.29到0.31的范围，来观察0.3π与030π的幅度

stem(w/pi,magX);

axis([0.44,0.46,0,0.5]);

xlabel('frequency between 0.44pi and 0.46pi');

运行结果

3. 总结

通过本次实验，我们使用matlab软件对数字信号处理课中的内容进行了实验仿真，得到了与理论课相近的结果；Matlab是一款非常好用，易上手的数学编程软件。Matlab软件拥有强大的函数库，通过调用丰富的函数，可以实现强大数字处理的功能。

通过本次实验，我们使用了Matlab实现了利用基2时间抽选的FFT算法，计算x(n) 的DFT X(k)，与X(k) 的理论值比较，并且通过结果发现，与实际要求相差不大。

由于Matlab软件功能太多，我们没有必要掌握所有的函数。通过本次课程，我们锻炼了自学的本领，能够在短时间内完成一个没怎么接触过的任务，锻炼了快速学习法，按照自己需要来学习，这也是我们大学生非常需要的一项本领。