【 MATLAB 】用 MATLAB 实现离散时间傅里叶变换（DTFT）的两个案例分析

先给出离散时间傅里叶变换的简单介绍：

如果 x(n) 是绝对可加的，即

$\sum_{- \infty}^{+ \infty} \left | x(n) \right | < \infty$

那么它的离散时间傅里叶变换给出为：

$X(e^{jw})= F[x(n)]=\sum_{n = - \infty}^{+ \infty}x(n)e^{-jwn}$

w 称为数字频率，单位是每样本 rad（弧度）或（弧度/样本）（rad/sample）

案例1：

求 $x(n) = (0.5)^nu(n)$ 的离散时间傅里叶变换 $X(e^{jw})$ ，并用MATLAB将 $X(e^{jw})$ 在 $[0,\pi]$ 之间的501个等分点上求值，并画出它的幅度、相位、实部和虚部。

题解：

由于x(n)是无限长的序列，所以不能直接用MATLAB直接从x(n)得到 $X(e^{jw})$ 。然而，我们可以用它对表达式 $X(e^{jw})$ 在 $[0,\pi]$ 频率点上求值，然后画出它的幅度和相位（实部或虚部）。

$X(e^{jw})=\sum_{-\infty}^{\infty}x(n)e^{-jwn}=\sum_{0}^{\infty}(0.5)^ne^{-jwn} =\sum_{0}^{\infty}(0.5e^{-jw})^n$

$=\frac{1}{1-0.5e^{-jw}}= \frac{e^{jw}} {e^{jw}-0.5}$

脚本：

clc
clear
close allw = [0:500]*pi/500; %[0,pi] axis divided into 501 points
X = exp(j*w)./( exp(j*w) - 0.5*ones(1,501) );
magX = abs(X);
angX = angle(X);
realX = real(X);
imagX = imag(X);subplot(2,2,1)
plot(w/pi,magX);
grid;
title('Magnitude Part');
xlabel('frequency in pi units');ylabel('Magnitude');subplot(2,2,2)
plot(w/pi,angX);
grid;
title('Angle Part')
xlabel('frequency in pi units');ylabel('Radians');subplot(2,2,3)
plot(w/pi,realX);
grid
title('Real Part');
xlabel('frequency in pi units');ylabel('Real');subplot(2,2,4);
plot(w/pi,imagX);
grid;
title('Imaginary Part');
xlabel('frequency in pi units');ylabel('Imaginary');

案例2：

求下面有限长序列的离散时间傅里叶变换：

$x(n)= \left \{ 1,2,3,4,5\right \}, -1 \leq n \leq 3$

在[0,pi]之间的501个等分频率上进行数值求值。

题解：

$X(e^{jw})=\sum_{-\infty}^{+\infty}x(n)e^{-jwn}=e^{jw}+2+3e^{-jw}+4e^{-j2w}+5e^{-j3w}$

我们可以直接对上式进行MATLAB编程，但是这种方法在有限长序列的DTFT中不是太方便，我们可以直接由 x(n) 来求它的DTFT。

我们使用向量化的编程方法，最后得到一个通用的公式。推导如下：

用MATLAB实现如下：

k = [0:M];
n = [n1:n2];
X = x * (exp(-j * pi/M)).^(n'*k);

给出MATLAB脚本语言如下：

clc
clear
close alln = -1:3;
x = 1:5;
k = 0:500;
w = (pi/500)*k;
X = x * (exp(-j * pi/500)).^(n' * k);magX = abs(X);
angX = angle(X);
realX = real(X);
imagX = imag(X);subplot(2,2,1);
plot(w/pi,magX);
title('Magnitude Part');
xlabel('w/pi');ylabel('Magnitude');subplot(2,2,2);
plot(w/pi,angX);
title('Angle Part');
xlabel('w/pi');ylabel('Radians');subplot(2,2,3);
plot(w/pi,realX);
title('Real part');
xlabel('w/pi');ylabel('Real');subplot(2,2,4);
plot(w/pi,imagX);
title('Imaginary Part');
xlabel('w/pi');ylabel('Imaginary');