现代数字信号处理课后作业【第三章】及实践-DIT2-FFT编程&谱分析

---

Author: Peter
Date:2020-11-6
Location:FZU

---

现代数字信号处理课后作业【第三章】

3-3 若周期实序列$x_p(n)$是n的偶函数，则$X_p(k)$也是实序列且为$k$的偶函数，试证明之。

解：依题目要求，即证明$X_p(k)=X_p(-k)=X_p^{\ast }(k)$

• 方法一：

  $X_p(k)=\sum _{n=0}^{N-1}{x}_p(n)W_N^{kn}$

  $X_p^{\ast }(k)=\sum _{n=0}^{N-1}[x_p(n)W_N^{kn}{]}^{\ast }=\sum _{n=0}^{N-1}{x}_p^{\ast }(n)W_N^{-kn}$

∵$x_p(n)是实序列，所以x_p^{\ast }(n)=x_p(n)$

  $X_p^{\ast }(k)=\sum _{n=0}^{N-1}{x}_p(n)W_N^{-kn}=X_p(-k)$

∵$x_p(n)是n的偶函数，所以x_p(n)=x_p(-n)$

  $X_p(-k)=\sum _{n=0}^{N-1}{x}_p(n)W_N^{-kn}=\sum _{n=0}^{N-1}{x}_p(-n)W_N^{-kn}=\sum _{n=1-N}^0{x}_p(u)W_N^{ku}$

     $=\sum _{n=0}^{N-1}{x}_p(u)W_N^{ku}=X_p(k)$

  $综上所述X_p(k)=X_p(-k)=X_p^{\ast }(k)，当x_p(n)为实偶函数时，X_p(k)也为实偶函数$

注: 上式中$\sum _{n=1-N}^0{x}_p(u)W_N^{ku}=\sum _{n=0}^{N-1}{x}_p(u)W_N^{ku}利用了周期性质$

• 方法2:

  $X_p(k)=\sum _{n=0}^{N-1}{x}_p(n)e^{-\frac{j2\pi kn}{N}}=\sum _{n=0}^{N-1}{x}_p(n)[cos(\frac{2\pi kn}{N})-jsin(\frac{2\pi kn}{N})]$

  $因为x_p(n)是偶函数，即x(n)=x(N-n)，x(n)关于n=\frac{N}{2}偶对称$

  $而sin(\frac{2\pi kn}{N})关于n=\frac{N}{2}奇对称，即sin(\pi )=0点$

  $因此X_p(k)=\sum _{n=0}^{N-1}{x}_p(n)cos(\frac{2\pi kn}{N})=X_p(-k)=X_p^{\ast }(k)$

  $综上所述X_p(k)=X_p(-k)=X_p^{\ast }(k)，当x_p(n)为实偶函数时，X_p(k)也为实偶函数$

3-4 设$x(n)=\{\begin{array}{ll}1& 0⩽n⩽3\\ 0& 其它\end{array}$  $h(n)=\{\begin{array}{ll}1& 4⩽n⩽6\\ 0& 其它\end{array}$${x}_p(n)=\sum _{r=-\infty }^{+\infty }x(n+7r)$  $h_p(n)=\sum _{r=-\infty }^{+\infty }h(n+7r)$
求$x_p(n)$与$h_p(n)$的周期卷积。

解：$两个序列周期相等，T=7，且大于x(n)、h(n)长度$

  $y_p(n)=\sum _{r=-\infty }^{+\infty }y(n+7r)$

  $y(n)=x(n)\ast h(n)，即x_p(n)主值区间与y_p(n)主值区间进行线性卷积$

  $其中x(n)=x_p(n)G_7(n)，h(n)=h_p(n)G_7(n)$

   根 据 不 进 位 乘 法 ， 易 得 y ( n ) = { 0 , 0 , 0 , 0 , 1 , 2 , 3 , 3 , 2 , 1 , 0 , 0 , 0 } 根据不进位乘法，易得y(n)=\{0,0,0,0,1,2,3,3,2,1,0,0,0\} 根据不进位乘法，易得y(n)={0,0,0,0,1,2,3,3,2,1,0,0,0}，注:两个长度为7的序列线性卷积后，长度为7+7-1=13。

  $y_p(n)为y(n)的周期延拓，如下图所示$

3-10 如果有限长序列，$x(n)=G_N(n)$
(1) 求$ZT[x(n)]$，并画出零极点图。
(2) 求频谱$X(e^{jw})$，并画出幅度特性和相位特性。
(3)求$DFT[x(n)]$，并画出幅度特性和相位特性 。

解：

• (1)
    $x(n)=u(n)-u(n-N)$

    $ZT[x(n)]=\sum _{n=-\infty }^{\infty }x(n)z^{-n}=\sum _{n=0}^{N-1}{z}^{-n}=\frac{1-z^{-N}}{1-z^{-1}}=1+z^{-1}+z^{-2}+...+z^{-(N-1)}$

    $由1-z^{-N}=0推出z^N=1$

     z = e j 2 π k N ， { k = 0 , 1 , . . . N − 1 } ， 即 单 位 圆 上 的 N 等 分 点 z=e^{\frac{j2\pi k}{N}}，\{k=0,1,...N-1\}，即单位圆上的N等分点 z=eNj2πk​，{k=0,1,...N−1}，即单位圆上的N等分点

    $由于k=0时，z=e^{\frac{j2\pi k}{N}}=1,分子分母零极对消$

     因 此 零 点 为 ： z = e j 2 π k N ， { k = 1 , . . . N − 1 } 因此零点为：z=e^{\frac{j2\pi k}{N}}，\{k=1,...N-1\} 因此零点为：z=eNj2πk​，{k=1,...N−1}

    $极点为z=0（N-1阶极点）$

   零极点图:

注： 如果有疑问为什么z=1即不是零点也不是极点，可参考频谱图-幅频特性$X(e^{jw}){|}_{w=0}$时的幅值即不为零也不为无穷大。

• (2)
    $X(e^{jw})=X(z){|}_{z=e^{jw}}=\frac{e^{jw}(1-e^{-jwN})}{e^{jw}-1}$

  $=\frac{sin\frac{wN}{2}}{sin\frac{w}{2}}\cdot e^{\frac{jw(1-N)}{2}}$

  $|X(e^{jw})|=\left|\frac{sin\frac{wN}{2}}{sin\frac{w}{2}}\right|其中|X(e^{j0})|=\frac{sa(\frac{wN}{2})}{sa(\frac{w}{2})}\cdot \frac{\frac{wN}{2}}{\frac{w}{2}}{|}_{w=0}=N$

  $\varphi (w)=\frac{w(1-N)}{2}当w=\frac{2k\pi }{N}，(k\ne N的整数倍时）X(e^{jw})=0，此时\varphi (w)=0$

   频谱图:

• Matlab代码块:

B=[1,1,1,1,1,1,1,1];
A=1;
subplot(2,2,1);
zplane(B,A);title('N=8时零极点图');
[H,w]=freqz(B,A,'whole');
subplot(2,2,3);
plot(w/pi,abs(H));
xlabel('\omega/\pi');ylabel('|X(e^j^\omega)|');grid on;title('N=8时，幅频特性');
subplot(2,2,4);
plot(w/pi,angle(H),'-');
xlabel('\omega/\pi');ylabel('\phi(\omega)');grid on;axis([0,2,-5,5]);title('N=8时，相频特性');

• (3)
    $X(k)=\sum _{n=0}^{N-1}x(n)e^{-\frac{j2\pi kn}{N}}，X(0)=\sum _{n=0}^{N-1}x(n)=N$

  $$\begin{gathered} 当k=1,2,....N-1时X(k)=DFT[X(n)]=X(e^{jw}){|}_{w=\frac{2\pi k}{N}}=\frac{e^{\frac{j2\pi k}{N}}[1-e^{-j2\pi k}]}{e^{\frac{j2\pi k}{N}}-1} \\ 因为e^{-j2\pi k}=1，所以X(k)=0k=1,2,...N-1 \end{gathered}$$

$$X(k)=\{\begin{array}{ll}N& k=0\\ 0& 1⩽k⩽N-1\end{array}$$

   频谱图:

• Matlab代码块:

xn=[1,1,1,1,1,1,1,1];N=8;
Xk=fft(xn,N);
figure;
x_k=0:N-1;
subplot(2,2,1);stem(x_k,abs(Xk),'.');
axis([-1,9,0,9]);grid on;title('|X_8(k)|');
subplot(2,2,2);stem(x_k,angle(Xk),'.');
grid on;title('arg[X_8(k)]');

3-12 已知$x(n)$是长为N的有限长序列并且$X(k)=DFT[x(n)]$
设$y(n)=\sum _{r=0}^{N-1}x(n+rN)$，求$DFT[y(n)]$与$X[k]$之间关系。

解：$X(k)=\sum _{n=0}^{N-1}x(n)W_N^{kn}$

  $y(n)长度为N^2，对y(n)进行N^2点DFT:$

  $DFT[y(n)]=\sum _{n=0}^{N-1}\sum _{r=0}^{N-1}x(n+rN)W_{N^2}^{k(n+rN)}=\sum _{r=0}^{N-1}[\sum _{n=0}^{N-1}x(n)e^{\frac{-j2\pi kn}{N^2}}]e^{\frac{-j2\pi kr}{N}}=\sum _{r=0}^{N-1}[\sum _{n=0}^{N-1}x(n)W_N^{\frac{kn}{N}}]e^{\frac{-j2\pi kr}{N}}$

  $=\sum _{r=0}^{N-1}X(\frac{k}{N})e^{\frac{-j2\pi kr}{N}}$

  $=X(\frac{k}{N})\sum _{r=0}^{N-1}{e}^{\frac{-j2\pi kr}{N}}$

$$DFT[y(n)]=X(\frac{k}{N})\sum _{r=0}^{N-1}{e}^{\frac{-j2\pi kr}{N}}=\{\begin{array}{ll}N\cdot X(\frac{k}{N})& \frac{k}{N}为整数\\ & \\ X(\frac{k}{N})\cdot \frac{1-e^{-j2\pi k}}{1-e^{\frac{-j2\pi k}{N}}}=0& \frac{k}{N}不为整数\end{array}$$
  由此可见，截取长度$M$等于$x_p(n)$的整数个周期，即$M=mN$，m为正整数时，$X_M(k)$也能表示$x_p(n)$的频谱结构，此题$m=N，M=N^2。$

4-3 某台计算机平均一次复乘需要100$us$，一次复加需要20$us$，今用来进行N=1024点DFT运算，问直接算法需要多少时间？时间抽奇偶分解FFT算法需要多少时间？

解：$直接法：X(k)=DFT[x(n)]=\sum _{n=0}^{N-1}x(n)W_N^{kn}$

  $从上式可以看出每计算X(k)的一个值，需要N次复数乘法和(N-1)复数加法$

  $所以T=N[100N+20(N-1)]，N=1024$

  $T_1=125.80864s$

  $采用时间抽奇偶分解FFT算法：复乘次数为\frac{N}{2}lo{g}_{2}N，复加次数为Nlo{g}_{2}N$

  $T_2=100\cdot \frac{N}{2}lo{g}_{2}N+20\cdot Nlo{g}_{2}N=0.7168s$

注: $时域基2FFT，N点FFT可以分为M=lo{g}_{2}N级，每级有\frac{N}{2}个蝶形$

  $一个蝶形需要1次复乘，2次复加，所以N点FFT总共有\frac{N}{2}M复乘，\frac{N}{2}\cdot 2M复加$

实验题

1、写一个基2时抽FFT程序，不使用系统自带的fft函数(Matlab)

• 编程思想:

 1.进行$N=2^M$点FFT，若输入序列长度不足N，则尾部补0直到总长度为N。

 2.输入序列的索引先转成二进制再经过倒序后，序列按照新的索引值进行重新排列，matlab自带函数fliplr可实现该功能。

 3.总共分成$M=lo{g}_{2}N$级进行蝶形运算，$L=1,...M$为当前级数,$B=2^{L-1}$ 为蝶形运算两个输入数据的距离即$X_L(k)=X_{L-1}(k)+X_{L-1}(k+B)W_N^P$，每级含有$B=2^{L-1}$个旋转因子$W_N^P=W_{2^L}^J，J=0,...B-1$，
所以$P=J\cdot 2^{M-L}$

 4.在当前级中寻找拥有相同旋转因子的蝶形进行运算，同一旋转因子对应着间隔为$2^L$点的$2^{M-L}$个蝶形，由于每个蝶形的两个输入数据只对计算本蝶形有用，意味着计算完一个蝶形后所有数据可立即存入原输入数据所占用的存储单元。

 5.循环计算当前级下一个旋转因子蝶形运算，存入下一级，直到M级输出结果

• 程序框图：

• DIT2-FFT函数代码块:

%----------------函数调用方法：Myfft(xn,N)------------------
%----------------xn为要进行FFT计算的序列，N为FFT点数=2^n------------
function return_num=Myfft(xn,N)   

M=log2(N);          %计算级数
len=length(xn);     %序列xn的长度

if(len<N)
    xn=[xn,zeros(1,N-len)];  %如果xn序列长度不足N点FFT，则补零到长度为N
end

%--------------------倒序部分调用自带函数fliplr------------------------
binary_index=dec2bin((0:N-1),M);  %将输入的n=0,1,...N-1转换成二进制
reverse_binary_index=fliplr(binary_index);  %将二进制序列逆序
reverse_decimal_index=bin2dec(reverse_binary_index)+1; %逆序完将索引值转为十进制，记得加1，matlab数组索引从1开始到N结束
xn=xn(reverse_decimal_index);   %按照新的索引值重新排序，完成逆序输入

%--------------------分级进行蝶形运算---------------------------------
for L=1:M                               %DIT-FFT变换，M级蝶形变换
    B=2^(L-1);                          %两个输入数据相距B个点
    for J=0:B-1;                        %旋转因子处理
        P=2^(M-L)*J;                    %旋转因子上标
        W=exp(-1i*2*pi*P/N);            %对应旋转因子
        for k=(J+1):2^L:N;              %每级相同旋转因子间隔，数组索引从1开始到N结束，每次循环间隔2^L点          
            T=xn(k)+xn(k+B)*W;          %进行蝶形运算,先暂存，此时的xn(k)还要参与下一个xn(k+B)计算，不能更新存储值
            xn(k+B)=xn(k)-xn(k+B)*W;    %等相减部分计算完存储后
            xn(k)=T;                    %再存回
        end
    end
end
return_num=xn;

• 倒序编程思想:

 1.从二进制倒序表可以看出规律，顺序数I加1，是在二进制数最低位加1，逢2向高位进位。而倒序数J则是在M位二进制数最高位加1，逢2向次高位(低位)进位。故可以从任一倒序值求得下一个倒序值。

 2.用order_i表示当前顺序数十进制值，back_j表示当前倒序数十进制值，$N=2^M$，M位二进制数最高位的十进制权值为$\frac{N}{2}$，且从左向右二进制位的权值依次为$\frac{N}{4}，\frac{N}{8}...$。因此最高位加1相当于十进制运算back_j+$\frac{N}{2}$，如果最高位为0，即back_j<$\frac{N}{2}$，则直接由back_j+$\frac{N}{2}$得出下一个倒序值。如果最高位为1，即back_j$⩽\frac{N}{2}$，则先将最高位变成0，即back_j=back_j-$\frac{N}{2}$,然后次高位加1：back_j+$\frac{N}{4}$，但次高位加1时还需和次高位权值$\frac{N}{4}$进行比较，循环判断，直到完成二进制最高位加1，逢2向右进位的运算。

 3.形成倒序数back_j后，将原数组$x(n)$中存放的输入序列重新按倒序排列，假设N=8时，输入序列是$x(0),x(1),x(2),....x(7)$，由上图可知，倒序后第一个序列值$x(0)$和$x(N-1)$不需要重排，所以顺序数order_i初值为1，终值为N-2，倒序数back_j初值为$\frac{N}{2}$。每计算出一个倒序值back_j，便与循环语句自动生成的顺序order_i比较，当order_i=back_j时，不需要交换，当二者不等时交换，为了避免重复调换只对order_i<back_j进行交换。

• 倒序程序框图：

• 含具体倒序方法的DIT2-FFT函数代码:

%----------------函数调用方法：Myfft2(xn,N)------------------
%----------------xn为要进行FFT计算的序列，N为FFT点数=2^n------------
function return_num=Myfft2(xn,N)

M=log2(N);          %计算级数
len=length(xn);     

if(len<N)
    xn=[xn,zeros(1,N-len)];  %如果xn序列长度不足N点FFT，则补零到长度为N
end

%------------------------自己编写的倒序法------------------------
back_j=N/2;                        %back_j为倒序初始值N/2，注意由于matlab数组索引从1开始到N结束，故xn(back_j+1)，xn(order_i+1)
for order_i=1:N-2                  %order_i为顺序计数，由于x(0)和x(N-1)不需要重排，故从x(1)到x(N-2)进行倒序
    if order_i<back_j              %如果顺序数小于倒序数则进行交换，二者相等时不用交换
        temp=xn(back_j+1);
        xn(back_j+1)=xn(order_i+1);
        xn(order_i+1)=temp;
    end                             
                                    %用当前倒序值计算下一个倒序值
    high_bit=N/2;                   %二进制最高位十进制权值为N/2
    while back_j>=high_bit          %当前倒序值如果大于最高位值
        back_j=back_j-high_bit;     %则要去掉最高位，判断次高位
        high_bit=high_bit/2;        %二进制次高位十进制权值为N/2/2，依次类推每到下一位权值为上一位权值一半
    end
    back_j=back_j+high_bit;         %直到"无进位时"，二进制最高位加1,完成下一位倒序计算
end

%--------------------分级进行蝶形运算---------------------------------
for L=1:M                               %DIT-FFT变换，M级蝶形变换
    B=2^(L-1);                          %两个输入数据相距B个点
    for J = 0:B-1;                      %旋转因子处理
        P=2^(M-L)*J;
        W=exp(-1i*2*pi*P/N);            %对应旋转因子
        for k=(J+1):2^L:N;              %每级相同旋转因子间隔，数组索引从1开始到N结束，每次循环间隔2^L点           
            T=xn(k)+xn(k+B)*W;          %进行蝶形运算,先暂存，此时的xn(k)还要参与下一个xn(k+B)计算，不能更新存储值
            xn(k+B)=xn(k)-xn(k+B)*W;    %等相减部分计算完存储后
            xn(k)=T;                    %再存回
        end
    end
end
return_num=xn;

• 测试代码块:

xn=[1 2 3 4 5 6 9 10];      %要转换的序列
N=32;                        %FFT点数 N=2^n
system_fft=fft(xn,N);       %matlab内置fft变换，用于对比
n=0:N-1;                    %x轴索引
subplot(3,1,1);stem(n,abs(system_fft),'r.');title(['系统自带函数',num2str(N),'点FFT']);
grid on; %添加网格
My_fft=Myfft(xn,N);           
subplot(3,1,2);stem(n,abs(My_fft),'.') ;title(['自己编写函数',num2str(N),'点FFT'])          
grid on; %添加网格
My_fft2=Myfft2(xn,N);
subplot(3,1,3);stem(n,abs(My_fft2),'.') ;title(['自己编写倒序函数',num2str(N),'点FFT'])          
grid on; %添加网格

• 运行效果:

2、研究如何利用FFT程序分析确定性时间连续信号，设输入信号为$x(t)=cos(2\pi Ft)$。
①F=50，N=32，T=0.000625
②F=50，N=32，T=0.005
③F=50，N=32，T=0.0046875
④F=50，N=32，T=0.004
⑤F=50，N=64，T=0.000625

• Matlab代码块:

F=50;		%输入信号频率
N=32;		%N点FFT
T=[0.000625,0.005,0.0046875,0.004];		%抽样间隔
n=0:63;%对连续时间信号x(t)采样64个点   %采样点数应足够大
figure('name','谱分析');
for i=1:5
        if i<5
            T1=T(i); 
        else 
            T1=T(1);
            N=64;
        end
        xn=cos(2*pi*F*n*T1);%将x(t)抽样成离散x(nT)
        Xk=Myfft2(xn,N);%对抽样得到的xn进行FFT变换，此处采用自己编写的FFT函数，可换成系统自带函数Xk=fft(xn,N);
        abs_Xk=abs(Xk);  %取模值
        xlabel_f=(0:N-1)/(N*T1);%谱分析的横轴坐标为f/Hz   谱分辨率=1/(N*T1)
        subplot(3,2,i);
        stem(xlabel_f(1:N/2),abs_Xk(1:N/2),'.');  %由于实序列的dft共轭对称，即模值关于k=N/2偶对称，所以只需查看前N/2个点即可  
        xlabel('频率/Hz');ylabel('振幅');grid on;title(['F=50','  ','N=',num2str(N),'  ','T=',num2str(T1)]);
end

• 实验效果：

• 总结:

 由于输入信号频率为$50Hz$，5种条件下的抽样频率分别为$1600Hz，200Hz，213.33Hz，250Hz，1600Hz$都满足奈奎斯特抽样速率，避免了频谱混叠。
 ①FFT点数N=32，抽样频率1600Hz，谱分辨率50Hz
 ②FFT点数N=32，抽样频率200Hz，谱分辨率6.25Hz
 ③FFT点数N=32，抽样频率213.33Hz，谱分辨率6.67Hz
 ④FFT点数N=32，抽样频率250Hz，谱分辨率7.81Hz
 ⑤FFT点数N=64，抽样频率1600Hz，谱分辨率25Hz
 从以上数据和实验图可以看出，当输入信号频率能被谱分辨率整除时，频谱图中的峰值处的频率便是输入的连续时间信号的频率。另从④③②三组数据可看出，当N为固定值时，降低采样频率可以提高谱分辨率（谱分辨率值越小说明分辨能力越强）。