实验二：数字信号的 FFT 分析
题目1
假设信号 x(n) 由下述信号组成：
?n(0?.45 x(n)?0.001*cos)?ns?in(0.3??n)?c os(0.3024)这个信号有两根主谱线 0.3pi 和 0.302pi 靠的非常近，而另一根谱线 0.45pi 的幅度很小，请选择合适的长度 N 和窗函数，用 DFT 分析其频谱，
得到清楚的三根谱线。 步骤：
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]; % row vector for n k = [0:1:N-1]; % row vecor for k WN = exp(-j*2*pi/N); % Wn factor
nk = n'*k; % creates a N by N matrix of nk values WNnk = WN .^ nk; % DFT matrix
Xk = xn * WNnk; % row vector for DFT coefficients
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('signal x(n), 0<=n<=999');xlabel('n'); 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('DTFT Magnitude');xlabel('frequency in pi units'); axis([0.29,0.31,0,500]);
xlabel('frequency between 0.29pi and 0.31pi'); axis([0.44,0.46,0,0.5]);
xlabel('frequency between 0.44pi and 0.46pi');
3.图片：
4.分析：
x(n)由3个正弦函数叠加而成，周期分别是40, 20, 1000。这里序列长度选择n=1000，为周期的最小公倍数。从频谱分析中可以看出，此时的谱线中有3条清晰的谱线，所在位置分别是0.3pi, 0.302pi和0.45pi。这说明当取样点数为函数周期的整倍数时，其频谱中科出现3条清晰的谱线。
题目2
已知信号
?Qn 0?n?N-1 x(n)??
?0 n?0, n?N
这里，N=25，Q= 0.9+j0.3。可以推导出 ，
1?QNX(k)??x(n)WN??(QWN)?，k?0,1?N?1 k1?QWNn?0n?0N?1N?1 nk kn首先根据这个式子计算X(k) 的理论值，然后计算输入序列x(n) 的32个值，再利用基2时间抽选的FFT算法，计算x(n) 的DFT X(k)，与X(k) 的
理论值比较（要求计算结果最少6位有效数字）。 步骤：
1.计算X(k) 的理论值 写入代码：
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)
得到结果：
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)] x=vap(x,6)
得到32点的x(n),输出结果为：
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]
4. 25点x(n)的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); magX=abs(X(1:1:25)); k=0:1:24; w=2*pi*k/25; stem(w/pi,magX);
title('25 points DFT Magnitude'); xlabel('frequency in 2pi units');
图像：
5. 32点x(n)的DFT图像 写入代码：
N=25;
x=[Q.^n,zeros(1,32-N)];
X=fft(x);
magX=abs(X(1:1:32)); k=0:1:31; w=2*pi*k/32; plot(w/pi,magX);
title('32 points FFT Magnitude'); xlabel('frequency in 2pi units');
图像：
6.分析
经过计算25点的DFT以及32点的DFT（亦即基2时域的FFT算法）我们可以发现，在6位有效数字的情况下，两个序列的第一个点是完全相同的，这是因为第一个点分布在Z平面单位圆与x轴交点处，但是之后由于点数的不同，两个序列的点在单位圆上分布的位置也不相同，因此可以看出会有较大的偏差。因为通过数字的判断不太好说明，因此画出了他们的图像，我们可以比较清楚地看到，两个序列的趋势以及整体的外包络线应该是相同的，只不过因为32点的外包络更光滑一些。
总结：
通过本次实验，我明白了分析时域离散周期信号的频谱时应该取其周期的整倍数进行分析，这样可以清楚地看出信号的频率分量。同时周期信号的频谱在取样点数合适的情况下，2pi内是镜像对称的，而非周期信号的频谱则不是。
本实验并未取更多的点，来比较取点数的多少对频谱分辨率的影响。其实倘若题目1中并未取够1000点或者取点数不是1000的整倍数，会出现频谱模糊的情况，因为分立的谱线横坐标无法精确地取到信号的三个频率分量。取点数越多，频谱分辨率就越高。而题目2由于序列本身只有25个非零点，是非周期信号，当补零后取32点的情况下，其实无法提高频谱的分辨率，只能得到更高密度的频谱，无法获取更多的信息。
