matlab数字信号处理实验——验证DFT变换的物理意义 超详细代码解析

  

 我们以这道题目为例，用matlab对离散序列进行DTFT变换以及DFT变换，绘制x(n)的连续频谱与10点DFT的离散频谱并做对比，目的是验证DFT变换的物理意义。

一、连续频谱函数的公式推导

   首先推导该序列的连续频谱函数，思路是由DTFT变换的定义式进行公式推导，需要用到等比数列求和的知识。

离散时间傅里叶变换（DTFT)，时域是非周期、离散信号，频域是连续、周期信号。DFT变换的变换公式如下:

  这是DTFT傅里叶正变换的公式，用于进行从时域到频域的转换，其中累加的范围是n从负无穷到正无穷，实际问题中n的范围可以更加给定时域序列的范围进行调整，以此来化简运算。

这是DTFT傅里叶负变换的公式，用于进行频域到时域的转换。

回归本题，由DTFT变换的定义式推导x(n)的连续频谱函数的过程如下：

 下面我们用matlab来求出10点DFT变换的离散频谱函数，并绘制DFT变换与DTFT变换得到的频谱函数图像。

二、本题的完整matlab代码与超详细注解

N=10; n=0:N-1;

%对xn的频谱函数采样2048个点可以近似看作是连续的频谱
%频谱函数的一个周期为2pi，对一个周期的频谱函数进行等间隔采样
%采样点数2048，产生具有2048个值的行向量，代表离散的w序列
%可以理解为由于采样间隔无限小而采样点数无限多，产生的频谱函数趋近于连续了
%因此Xw可认为近似等于xn经过DTFT序列傅里叶变换后得到的连续频谱函数
w=2*pi*(0:2047)/2048;   

%连续频谱函数通过DTFT变换的定义推导出，推导过程见附图
%由于w是一个具有2048个点的离散序列，这里要注意用点乘
%点乘表示矩阵（这里是行向量）中元素与元素对应相乘，结果是对应的一个序列
%如果写成乘法，matlab会执行矩阵与矩阵相乘，结果是一个值，就错了
Xw=...
   (1/2).*( (1-exp(-1i*(w+pi/(2*N))*N))./(1-exp(-1i*(w+pi/(2*N))))...
   +(1-exp(-1i*(w-pi/(2*N))*N))./(1-exp(-1i*(w-pi/(2*N)))));

xn=cos(n*pi/(2*N));   %产生xn，用于计算Xk
Xk=fft(xn,N);         %用fft快速傅里叶变换计算N点的DFT

%注意这里是画幅度对于w/pi的曲线，对w进行了去弧度，此时w的单位成为了1/s
subplot(2,2,1);plot(w/pi,abs(Xw));grid on;             %分区作图，在绘图时显示网格线
xlabel('w/pi');title('DTFT变换得到的连续频谱的幅度曲线');

%绘制连续频谱一个周期（0~2pi)的曲线，由于对w进行了去弧度，横坐标范围就是0~2
subplot(2,2,2);plot(w/pi,angle(Xw));grid on;
%在点[0,N]和[0,0]之间绘制一条紫色虚线。将 Color 和 LineStyle 属性设置为名称-值对组。
axis([0,2,-pi,pi]);line([0,N],[0,0],'Color','magenta','lineStyle','--');
xlabel('w/pi');title('DTFT变换得到的连续频谱的相位曲线');


k=0:N-1;
subplot(2,2,3);stem(k,abs(Xk),'filled');grid on;
xlabel('k(w=2pi*k/N)');ylabel('|X(k)|');  %DFT隐含周期性，是DTFT变换0~2pi周期内的N等分
%在DFT离散频谱上叠加连续频谱Xw的幅度曲线
hold on;
%两个函数绘制在一个图里时一定要注意单位变换，与横轴单位一致
%DFT隐含周期性，是DTFT变换0~2pi周期内的N等分,所以w=2pi*k/N
plot(N/2*w/pi,abs(Xw),'--'); 
title('在DFT频谱上叠加DTFT连续频谱的幅度曲线');

subplot(2,2,4);stem(k,angle(Xk),'filled');grid on;
axis([0,N,-pi,pi]);line([0,N],[0,0],'Color','magenta','lineStyle','--');
xlabel('k(w=2pi*k/N)');ylabel('Arg[X(k)]');
hold on;
%在DFT离散频谱上叠加连续频谱的相位曲线
plot(N/2*w/pi,angle(Xw),'--');
title('在DFT频谱上叠加DTFT连续频谱的相位曲线');

三、代码中值得注意的地方注解

• 在绘制DTFT变换的连续频谱时，选取“w/pi”而非“w”作为横坐标，这样绘制出来的图像是连续频谱函数对于w/pi的幅度曲线。由于w=2*pi*k/N,(其中N=2048为采样点数，k为0~2047), 单位为（rad/s）,所以w/pi的单位为（1/s），且范围为0~2。如下面这段代码所示：

  %注意这里是画幅度对于w/pi的曲线，对w进行了去弧度，此时w的单位成为了1/s
  subplot(2,2,1);plot(w/pi,abs(Xw));grid on;             %分区作图，在绘图时显示网格线
  xlabel('w/pi');title('DTFT变换得到的连续频谱的幅度曲线');
  
  %绘制连续频谱一个周期（0~2pi)的曲线，由于对w进行了去弧度，横坐标范围就是0~2
  subplot(2,2,2);plot(w/pi,angle(Xw));grid on;
  %在点[0,N]和[0,0]之间绘制一条紫色虚线。将 Color 和 LineStyle 属性设置为名称-值对组。
  axis([0,2,-pi,pi]);line([0,N],[0,0],'Color','magenta','lineStyle','--');
  xlabel('w/pi');title('DTFT变换得到的连续频谱的相位曲线');
  

  那如果我们不对w作去弧度化，就绘制连续频谱对于w的幅度曲线，则需要修改代码如下：

  %注意这里是画幅度对于w的曲线
  subplot(2,2,1);plot(w,abs(Xw));grid on;             %分区作图，在绘图时显示网格线
  xlabel('w');title('DTFT变换得到的连续频谱的幅度曲线');
  
  %绘制连续频谱一个周期（0~2pi)的曲线
  subplot(2,2,2);plot(w,angle(Xw));grid on;
  %在点[0,N]和[0,0]之间绘制一条紫色虚线。将 Color 和 LineStyle 属性设置为名称-值对组。
  axis([0,7,-pi,pi]);line([0,N],[0,0],'Color','magenta','lineStyle','--');
  xlabel('w');title('DTFT变换得到的连续频谱的相位曲线');

  做一个周期内的频谱函数曲线图，0~2pi约等于0~6.18，所以横坐标范围设置在0~7就足够且合适。我们将两种情况所绘制出来的图像作对比：

  因此，选取“w/pi”而非“w”作为横坐标显然能够帮助我们更清晰地观察图像，由于pi是一个值不确定的无限不循环小数，所以以“w”作横坐标时我们很难准确地定位图像对称轴，范围线的位置，而以w/pi作横坐标，由于除去了不确定的pi值，曲线对称轴的位置，即周期的一半对应“pi/pi”，也就是’1’,而一个周期结束的位置对应”2pi/pi”也就是“2 ”。 因此转换横坐标除pi的意义是让曲线图更具分析和观察的优势。

• 本题在绘制DTFT变换图像时，由于DTFT变换的频域是连续的信号，我们采用了无限密集的离散序列去逼近连续信号的思想来做图。

  w=2*pi*(0:2047)/2048; 

  对xn的频谱函数采样2048个点可以近似看作是连续的频谱，频谱函数的一个周期为2pi，对一个周期的频谱函数进行等间隔采样，采样点数为2048，产生了一个有2048个值的行向量，代表离散的w序列。可以理解为由于采样间隔无限小而采样点数无限多，产生的频谱函数趋近于连续了。因此Xw可认为近似等于xn经过DTFT序列傅里叶变换后得到的连续频谱函数。 这个思想和高数里面的“无限分割求积分”有点类似，前者是将连续的曲线分割为多个长度无限小的区间，而在足够足够小的距离内，一小段曲线被近似等于一小段直线，整个曲线与坐标轴围成的图形被分割为无数个无限小的矩形，而总的面积（即要求的积分的值）就可以转换为求一连串儿矩形的面积之和。 而这里我们在绘制DTFT变换得到的连续频谱函数图像时，而是对一个周期的函数进行疯狂采样（类似于分割），频域的采样间隔（类似于区间长度）无限小，产生一个无限密集的离散序列，它的包络就非常非常逼近连续曲线了。

• Xw=...
     (1/2).*( (1-exp(-1i*(w+pi/(2*N))*N))./(1-exp(-1i*(w+pi/(2*N))))...
     +(1-exp(-1i*(w-pi/(2*N))*N))./(1-exp(-1i*(w-pi/(2*N)))));

  这个式子,emmm，就是推导正确了也很难在代码中写对，很容易少敲一两个符号，出现“括号不匹配”的错误，一定要小心。   另外我们在用matlab进行数字信号处理实验的时候一定要注意“点乘”与“乘”的区别，很容易犯错！“乘号”在matlab中执行的是矩阵与矩阵相乘，而矩阵相乘的结果是一个数值。而matlab中“点乘”表示矩阵中元素与元素对应相乘，结果是对应的一个序列。  一个避免犯错的简单思想就是，“看你要的是一个序列还是一个数值”，比如这个题，我们意在产生一个点数无限多（有2048个点）的离散序列去逼近连续频谱，那肯定用点乘法啊。再者，由于w是一个具有2048个点的行向量，包括w序列乘除的地方都注意要用点乘。而像2*N ，N*pi这种单个数值之间的乘法，用乘号还是点乘号就都可以。（点除和除的区别类似，不赘述了）

• subplot(2,2,4);stem(k,angle(Xk),'filled');grid on;
  axis([0,N,-pi,pi]);line([0,N],[0,0],'Color','magenta','lineStyle','--');
  xlabel('k(w=2pi*k/N)');ylabel('Arg[X(k)]');
  hold on;
  %在DFT离散频谱上叠加连续频谱的相位曲线
  plot(N/2*w/pi,angle(Xw),'--');
  title('在DFT频谱上叠加DTFT连续频谱的相位曲线');

  这里用hold on；在DFT离散频谱上叠加了连续频谱的相位曲线，既然是在人家的图上面叠加，那一定要注意单位变换，要与人家的单位一致。这里横轴的单位是k ，一个0~2047范围内的离散序列，其实就是对连续频谱采样时的采样点们。

  w=2*pi*(0:2047)/2048; 

  DFT变换隐含周期性，是DTFT变换得到的连续频谱在0~2pi周期内的N点等间隔采样,所以w=2pi*k/N。N=2048, k的范围是0~2047，所以k=N/2*w/pi。至此，连续频谱函数Xw的单位由w被转换为k，与横轴单位一致了，它才能被正确地叠加在人家的图里啊。

四、实验结论

 这是运行代码后绘制成功的图像，可以看出，DFT相当于对DTFT的一个周期进行等间隔采样。本实验验证了离散傅里叶变换DFT的物理意义，即有限长序列的离散频域表示。需要注意的是，DFT变换隐含周期性，它所处理的有限长序列都是作为周期序列的一个周期来表示的，是DTFT变换得到的连续频谱在0~2pi周期内的N等分。