DSP第二次上机作业

DSP第二次上机作业

一、实验目的

  学习并掌握用MATLAB程序计算DFT的方法，深入理解DFT的概念和性质。

  数字信号处理中的五大变换：
  傅里叶级数（FS）：将周期性连续函数变换为离散频率点上的函数
  离散傅里叶级数（DFS）：将周期性离散函数变换为离散频率点上的周期性函数
  傅里叶变换（FT）：将连续函数变换为连续频率的函数
  离散时间傅里叶变换（DTFT）：将离散函数变换为连续频率的函数
  离散傅里叶变换（DFT）：将有限长离散函数变换为离散频率点上的函数

  快速傅里叶变换（FFT）：DFT的一种快速、高效的算法。

二、Task 1

2.1 题目要求

  计算序列$$x(n)=\{\begin{array}{rlr}0.1n+1,& & 0\le n\le 15\\ 0,& & 其他\end{array}$$的N=16及N=32点的DFT，即$X(k)=DFT[x(n)]$。分别输出$\left|X(k)\right|$，$\varphi (k)$曲线，并与FFT函数计算结果比较。

2.2 程序设计

clear;clc;
% x=input('从左向右输入序列x(n)的幅度值矩阵（默认x(n)起始位置为0）：'); 
x=zeros(1,16);
for i=0:15
    x(i+1)=0.1*i+1;
end
N=input('输入DFT的变换区间：');
N_x=length(x);%x(n)的长度

%画出x(n)图像
nx=0:N_x-1;
subplot(4,1,1);
stem(nx,x,'.');title('输入序列x(n)');xlabel('n');ylabel('x(n)')

%若x(n)的长度小于N，则对x(n)补零
if N_x<N
    x1=zeros(1,N-N_x);
    x=[x,x1];
end

%画出补零后x(n)图像
n=0:N-1;
subplot(4,1,2);
stem(n,x,'.');title('补零后的序列x(n)');xlabel('n');ylabel('x(n)')

%N点DFT
X=zeros(1,N);
Xm=zeros(1,N);  %幅值
fen=zeros(1,N); %相角
re=zeros(1,N);
im=zeros(1,N);
for k=n
    for j=n
        re(k+1)=re(k+1)+x(j+1)*cos(-2*pi/N*(k)*(j));%实部
        im(k+1)=im(k+1)+x(j+1)*sin(-2*pi/N*(k)*(j));%虚部
    end
    X(k+1)=re(k+1)+1i*im(k+1);
    Xm(k+1)=sqrt(re(k+1)^2+im(k+1)^2);
    fen(k+1)=atan(im(k+1)/re(k+1));
end

subplot(4,1,3);
stem(n,Xm,'.');title('幅频特性|X(k)|');xlabel('k');ylabel('Xm(k)')  %画出X(k)的幅度Xm(k)图像
subplot(4,1,4);
stem(n,fen,'.');title('相频特性fen(k)');xlabel('k');ylabel('fen(k)')  %画出X(k)的相角fen(k)图像

fprintf("X(k)的值为：");
X
fprintf('\n');

% 使用FFT函数
X1=fft(x,N);
fprintf("X1(k)的值为：");
X1
fprintf('\n');

2.3 程序结果

  N=16点DFT：
  与FFT函数的输出结果比较：


  N=32点DFT：

  与FFT函数的输出结果比较：

2.4 结果分析

  由输出结果图可知，作者根据DFT的定义式编写的DFT函数与MATLAB中的FFT函数输出结果相同，证明作者编写的DFT函数结果是正确的。其中，当输入序列x(n)的长度小于变换区间长度N时，需要在x(n)的末端补零，将x(n)的长度扩展为N后再做变换。
  但是，根据定义式直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时运算量过大，从而限制了DFT在信号频谱分析和实时信号处理中的应用。FFT是DFT的一种快速、高效率的算法，得到相同结果的同时大大减少了算法的运算量，例如DIT-FFT和DIF-FFT的复数乘法运算次数为$\frac{N}{2}lo{g}_{2}N$,复数加法运算次数为$Nlo{g}_{2}N$，N越大，FFT算法的效率越高。

三、Task 2

3.1 题目要求

  输出序列$\sin (0.5\pi n+0.2\pi )=x(n),0\le n\le N-1$ ，N自定。计算并输出x(n)的N点DFT：$X_1(k)$，以及x(n)的2N点DFT： $X_2(k)$。
  观察 $X_1(k)$ 和 $X_2(k)$，能得出什么结论？

3.2 程序设计

clear;clc;
% x=input('从左向右输入序列x(n)的幅度值矩阵（默认x(n)起始位置为0）：'); 
N=16; % 自定义：N=16
x=zeros(1,N);
for i=0:N-1
    x(i+1)=sin(0.5*pi*i+0.2*pi);
end

%N点DFT
n1=0:N-1;
x1=x;

X1=zeros(1,N);
Xm1=zeros(1,N);  %幅值
fen1=zeros(1,N); %相角
re1=zeros(1,N);
im1=zeros(1,N);
for k=n1
    for j=n1
        re1(k+1)=re1(k+1)+x1(j+1)*cos(-2*pi/N*(k)*(j));%实部
        im1(k+1)=im1(k+1)+x1(j+1)*sin(-2*pi/N*(k)*(j));%虚部
    end
    X1(k+1)=re1(k+1)+1i*im1(k+1);
    Xm1(k+1)=sqrt(re1(k+1)^2+im1(k+1)^2);
    fen1(k+1)=atan(im1(k+1)/re1(k+1));
end

subplot(3,2,1);
stem(n1,x1,'.');title('N点DFT输入');xlabel('n');ylabel('x1')  %画出x1(k)的图像
subplot(3,2,3);
stem(n1,Xm1,'.');title('幅频特性|X1(k)|');xlabel('k');ylabel('Xm1(k)')  %画出X1(k)的幅度Xm1(k)图像
subplot(3,2,5);
stem(n1,fen1,'.');title('相频特性fen1(k)');xlabel('k');ylabel('fen1(k)')  %画出X1(k)的相角fen1(k)图像

%2N点DFT
n2=0:2*N-1;
x0=zeros(1,N); % 对x(n)补零
x2=[x,x0];

X2=zeros(1,2*N);
Xm2=zeros(1,2*N);  %幅值
fen2=zeros(1,2*N); %相角
re2=zeros(1,2*N);
im2=zeros(1,2*N);
for k=n2
    for j=n2
        re2(k+1)=re2(k+1)+x2(j+1)*cos(-2*pi/(2*N)*(k)*(j));%实部
        im2(k+1)=im2(k+1)+x2(j+1)*sin(-2*pi/(2*N)*(k)*(j));%虚部
    end
    X2(k+1)=re2(k+1)+1i*im2(k+1);
    Xm2(k+1)=sqrt(re2(k+1)^2+im2(k+1)^2);
    fen2(k+1)=atan(im2(k+1)/re2(k+1));
end

subplot(3,2,2);
stem(n2,x2,'.');title('2N点DFT输入');xlabel('n');ylabel('x2')  %画出x2(k)的图像
subplot(3,2,4);
stem(n2,Xm2,'.');title('幅频特性|X2(k)|');xlabel('k');ylabel('Xm2(k)')  %画出X2(k)的幅度Xm2(k)图像
subplot(3,2,6);
stem(n2,fen2,'.');title('相频特性fen2(k)');xlabel('k');ylabel('fen2(k)')  %画出X2(k)的相角fen2(k)图像

3.3 程序结果

3.4 结果分析

  从输出结果图可以看出，变换区间N扩大2倍，采样频率不变时，频域的频谱间隔减小了一半，即频域的采样密度增加了一倍。因为$X(k)$是$x(n)$的离散傅里叶变换$X(e^{j\omega })$在$[0,2\pi )$区间上的N点等间隔采样，所以在原序列末尾补零相当于增加了频域对$X(e^{j\omega })$的采样点数，使真实频谱的包络能够更好地显现出来。

四、Task 3

4.1 题目要求

  用快速卷积算法计算下列两序列的线性卷积序列，并输出结果：
$$x_1(n)=0,2,2,1$$ $$x_2(n)=\{\begin{array}{rlr}(1.02{)}^n,& & 0\le n\le 15\\ (0.98{)}^n,& & 16\le n\le 28\end{array}$$   输出$x_1(k)$、$x_2(k)$及其FFT信号图形，输出卷积结果。（注意FFT点数应满足循环卷积与线性卷积相等条件）与第一次上机作业中时域线性卷积比较（指计算时间的比较，卷积结果应相同）。

4.2 程序设计

clear;clc;
x1=[0,2,2,1];
x2=zeros(1,29);
for i=0:28
    if i<=15
    x2(i+1)=1.02^i;
    else
    x2(i+1)=0.98^i; 
    end
end

%% 快速线性卷积算法
L=32; % L>=N1+N2-1=32，取L=32
x01=zeros(1,28);
x11=[x1,x01];% x1补28个零点
x02=zeros(1,3);
x22=[x2,x02];% x1补3个零点

X1=fft(x11,L);% L点DFT
X2=fft(x22,L);

yy=ifft(X1.*X2);
for i=0:31
    fprintf("%.2f ",yy(i+1));
end
fprintf('\n');

n=0:31;n1=0:3;n2=0:28;
subplot(3,2,1);
stem(n1,x1,'.');title('X1(k)图形');xlabel('k');ylabel('X1(k)');
subplot(3,2,3);
stem(n2,x2,'.');title('X2(k)图形');xlabel('k');ylabel('X2(k)');
subplot(3,2,5);
stem(n,yy,'.');title('y(n)图形');xlabel('n');ylabel('y(n)');

%% 第一次上机作业中的时域线性卷积函数
x=x1;
h=x2;
Start_x=0;Start_h=0;
N_x=length(x);%x(n)的长度
N_h=length(h);%h(n)的长度
nx=Start_x:Start_x+N_x-1;
nh=Start_h:Start_h+N_h-1;

subplot(3,2,2);
stem(nx,x,'.');title('x(n)图形');xlabel('n');ylabel('x(n)');    %画出x(n)图形
subplot(3,2,4);
stem(nh,h,'.');title('h(n)图形');xlabel('n');ylabel('h(n)');    %画出h(n)图形

Start=Start_x+Start_h;    %卷积后序列起始位置为x(n)的起始位置加h(n)的起始位置
N=N_x+N_h-1;    %卷积后序列的长度为x(n)的长度Nx与h(n)的长度Nh的和减一
n=Start:Start+N-1;
y=zeros(1,N);   %初始化y

for i=n    %线性卷积
    sum=0;
    for m=nx
        if i-m>=Start_h && i-m<=Start_h+N_h-1
            y(i-Start+1)=y(i-Start+1)+x(m-Start_x+1).*h(i-m-Start_h+1);
        end
    end
end

for i=n
    fprintf('%.2f ',y(i-Start+1));
end
fprintf('\n');

subplot(3,2,6);
stem(n,y,'.');title('y(n)图形');xlabel('n');ylabel('y(n)');

4.3 程序结果

4.4 结果分析

  当DFT的运算区间长度L满足$L\ge N_1+N_2-1$时，序列$x_1(n)$和$x_2(n)$的线性卷积结果和循环卷积结果相等。快速线性卷积算法，即在满足$L\ge N_1+N_2-1$的条件下，先在两输入序列后补零至长度为L，然后分别作L点DFT，得到结果后相乘，再根据DFT频域相乘，时域循环卷积的性质，将相乘后的结果作IDFT即得到输入两序列的线性卷积。从结果图中看出，根据快速线性卷积算法编写的程序与套用上机作业一中编写的线性卷积程序输出结果完全相同。

五、Task 4

5.1 题目要求

  若$x(n)=x_1(n)+j{x}_2(n)$，其中$x_1(n)=cos\frac{\pi }{4}n$，$x_2(n)=sin\frac{\pi }{8}n$。根据DFT的对称性，由$X(k)$求出$X_1(k)=DFT[x_1(n)]$和$X_2(k)=DFT[x_2(n)]$。（$X(k)=DFT[x(n)]$）

5.2 程序设计

clear;clc;
N=16; %设 x(n)=x1(n)+jx2(n) 长度为N
x1=zeros(1,N);
x2=zeros(1,N);
x=zeros(1,N);
for i=0:N-1
    x1(i+1)=cos(pi/4*i);
    x2(i+1)=sin(pi/8*i);
    x(i+1)=x1(i+1)+1i*x2(i+1);
end

X=fft(x,N);
XX=conj(X);   % X(k)共轭
XX=[XX,XX(1)];
Xep=zeros(1,N);
Xop=zeros(1,N);
for i=0:N-1
    Xep(i+1)=1/2*(X(i+1)+XX(N-i+1)); % 共轭对称分量
    Xop(i+1)=1/2*(X(i+1)-XX(N-i+1)); % 共轭反对称分量
end

X1=Xep;
X2=Xop/1i; 

% 利用fft函数计算x1(n)和x2(n)的DFT
X11=fft(x1,N);
X22=fft(x2,N);

%输出结果并比较
fprintf("X1(k)的值为：");
X1
fprintf("利用fft函数得到X1(k)的值为：");
X11
fprintf("X2(k)的值为：");
X2
fprintf("利用fft函数得到X2(k)的值为：");
X22

5.3 程序结果

5.4 结果分析

  根据DFT的对称性质，如果序列$x(n)$的DFT为$X(k)$，则$x(n)$的实部和虚部乘以$j$的DFT分别为$X(k)$的共轭对称分量和共轭反对称分量，而$x(n)$的共轭对称序列和共轭反对称序列的DFT分别为$X(k)$的实部和虚部乘以$j$。因此，在已知$X(k)$的情况下，要求出$x_1(n)$即$x(n)$实部的DFT，求出$X(k)$的共轭对称分量即可；要求出$x_2(n)$即$x(n)$虚部的DFT，需利用DFT的线性性质，先求出$X(k)$的共轭反对称分量，再除以$j$后得到。

5.4 结果分析

  根据DFT的对称性质，如果序列$x(n)$的DFT为$X(k)$，则$x(n)$的实部和虚部乘以$j$的DFT分别为$X(k)$的共轭对称分量和共轭反对称分量，而$x(n)$的共轭对称序列和共轭反对称序列的DFT分别为$X(k)$的实部和虚部乘以$j$。因此，在已知$X(k)$的情况下，要求出$x_1(n)$即$x(n)$实部的DFT，求出$X(k)$的共轭对称分量即可；要求出$x_2(n)$即$x(n)$虚部的DFT，需利用DFT的线性性质，先求出$X(k)$的共轭反对称分量，再除以$j$后得到。
  由程序结果图可以看出，根据上述思路编写程序得到的结果和用fft函数直接计算$x_1(n)$和$x_2(n)$的DFT得到的结果完全相同，说明该思路是正确的。