离散信号的傅里叶分析

发表于 分类于 Waline： 本文字数： 3.9k 阅读时长 ≈ 14 分钟

基础知识

定义

本文讨论的离散信号指的是离散时间信号（discrete-time signal），信号只定义在离散的时刻点上，也就是自变量（时间）仅取在一组离散值上。例如，股票市场的指数就是离散时间信号。

在这里需要注意几个容易混淆的概念。如果信号的自变量和函数值都取连续值，称信号为模拟信号或连续时间信号；如果自变量取离散值，而函数值取连续值，则称信号为离散时间信号，这种信号通常来源于对模拟信号的采样；如果信号的自变量和函数值均取离散值，则称为数字信号，数字信号其实就是量化后的离散时间信号，是计算机中的表示形式。

离散信号的表示和连续信号类似，只不过把时间变量改成了。同样地，我们基于周期复指数信号讨论正弦信号的特性，它的表达式为

它可以看成是连续信号在整数点上的取值，和连续信号类似，但又有区别。可以使用如下的MATLAB代码绘制离散信号的图像：

1
2
3
4
n = -5:5;
x = sin(pi*n/5);
stem(n, x, '.'); line([-5, 5], [0, 0]);
axis([-5, 5, -1.2, 1.2]); xlabel('n'); ylabel('x[n]')

振荡速率

在连续信号情况下，每个对应的是不同的信号，但是离散情况下就不同了，由于

所以，离散周期复指数信号在频率与频率时是完全一样的。

所以，考虑这种离散时间周期复指数信号时，仅仅需要在某个间隔内选择即可，一般取或。由上面的讨论可知，数值增加时，的振荡速率并不会不断增加。从增加到时，振荡速率越来越快，从变化到时，振荡速率会下降。下图是取不同值时对应正弦信号的图像：

图1 不同角频率下的图像

所以，离散复指数信号的低频部分位于在的偶数倍值附近，高频部分则位于在的奇数倍值附近。

周期公式

现在讨论它的周期。根据周期信号的定义有

所以，必须有个整数，满足

或

所以，只有为有理数时，才是周期信号，否则就不是周期的。这是由于离散时间信号仅能在自变量的整数值上有定义造成的。

基波频率

那么基波周期也可以写成

基于上面的原因，通常给定周期复指数信号的基波周期，然后求得频率（此时），得到它的表达式。

谐波信号

再讨论成谐波关系的周期离散时间复指数信号。下面的信号具有公共周期，频率都是基波周期的整数倍：

在连续情况下，成谐波关系的信号都是不相同的。但是离散情况下则不一样，因为

所以，上述成谐波关系的信号中，仅有个互不相同的周期复指数信号，而且

离散傅里叶级数（DFS）

离散时间周期信号可以表示为成谐波关系的指数信号的加权和，即为离散傅里叶级数（discrete Fourier series, DFS）。由于周期为、基波频率、成谐波关系的周期复指数信号中只有个是不相同的，所以离散时间傅里叶级数的形式为有限个复指数信号的和

求和限表示在个相继整数的区间上变化，下同。系数的计算公式为

很多教材求和限取到区间，此时、可写成

由于成谐波关系的信号中仅有个互不相同的周期复指数信号，所以，离散时间周期信号的傅里叶级数一定以为周期往复变化的，也就是说，离散时间周期信号的频谱是离散的、周期的。和连续时间信号傅里叶级数不同的是，级数的项数是有限的，不存在收敛问题，而且能够完全精确地逼近。

如果令，则DFS公式对、还可以写成如下形式

下面的MATLAB代码实现了dfs和idfs函数，用于计算离散傅里叶级数及其逆运算，可以任意指定求和限：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
function [Xk, k] = dfs(xn, n1, n2)
% Compute Discrete Fourier Series Coefficients
% Xk = DFS coeff. array over n1 <= k <= n2
% xn = One period of periodic signal over n1 <= n <= n2
N = n2 - n1 + 1;       % Fundamental period of xn
n = n1:n2;             % row vector for n
k = n1:n2;             % row vector for k
WN = exp(-1j*2*pi/N);  % Wn factor
nk = n' * k;           % creates a N by N matrix of nk values
WNnk = WN .^ nk;       % DFS matrix
Xk = xn * WNnk;        % row vector for DFS coefficients
end

function [xn, n] = idfs(Xk, k)
% Compute Inverse Discrete Fourier Series
% xn = One period of periodic signal over n
% Xk = DFS coeff. array over k
n = k;
N = length(n);         % Fundamental period
WN = exp(-1j*2*pi/N);  % Wn factor
nk = n' * k;           % creates a N by N matrix of nk values
WNnk = WN .^ (-nk);    % DFS matrix
xn = (Xk * WNnk) / N;  % row vector for DFS coefficients
end

来看一个例子。现有一个离散周期信号，其周期为，某个周期片段如下：从到的值为，其余为，即一个周期内，有个时间点的值为。如下的代码取周期信号在一个周期的片段，使用dfs函数计算傅里叶级数系数，本例得到的频谱系数都是实数。最后绘制信号及其频谱，它们都是周期的，在代码中是通过周期延拓得到的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
N = 20;
N1 = 2; n1 = -N1; n2 = N - N1 - 1; n = n1:n2;

% xn为一个周期片段
xn = [ones(1, 2*N1+1), zeros(1, N-2*N1-1)];
[Xk, k] = dfs(xn, n1, n2);

% 延拓成3个周期
xnh = [xn, xn, xn];
nh = [n-length(n), n, n+length(n)];
Xkh = [Xk, Xk, Xk];
kh = [k-length(k), k, k+length(k)];

% 绘制周期信号的图像
subplot(2, 1, 1); stem(nh, xnh, '.'); line([min(nh), max(nh)], [0, 0]);
axis([min(nh), max(nh), min(xnh)-1, max(xn)+1]);
% 绘制周期信号的频谱图
subplot(2, 1, 2); stem(kh, real(Xkh), '.'); line([min(kh), max(kh)], [0, 0]);
axis([min(real(kh)), max(kh), min(real(Xkh))-1, max(real(Xkh))+1]);

信号及其频谱如下图所示

图2 周期离散信号及其频谱

读者可以尝试增大的值，看频谱图如何变化。

离散时间傅里叶变换（DTFS）

离散时间非周期信号可以表示为不同频率的指数信号的加权和，它的频谱是连续的，用离散时间傅里叶变换（discrete-time Fourier transform, DTFT）描述，后面简称DTFT。和连续信号的分析方法类似，不断增大周期信号的周期，频谱会越来越密集，这一点可以使用上面的代码验证。当时，频谱由离散变为连续，周期信号变为非周期信号，傅里叶级数转化为傅里叶变换。得到了如下的傅里叶变换公式对

从DTFT得到的频谱是连续函数，由于相差的两个复指数信号是一样的，所以是周期函数，周期为。也就是说，离散时间非周期信号的频谱是连续的、周期的。频谱中靠近偶数倍的，对应的是低频率的信号；而靠近奇数倍的，对应的则是高频率的信号。

举个例子，取上图中周期信号的一个片段：从取到，一共取个点。那么该信号的DTFT计算如下：

一般来说是复函数，但此例中得到的是实函数，可直接画出它的图像，注意横轴标度是以为单位。

图3 DTFT是周期连续函数

DTFS与DFS的关系

与图2的频谱图作对比可发现，频谱图的包络线就是本例中的周期连续函数。实际上，可以证明，将任意非周期函数的离散时间傅里叶变换作等间隔采样，可得到的周期延拓的傅里叶级数系数，采样间隔为，即

令，，上述关系还可以记为

DTFS与DFS的这种关系是后续引入离散傅里叶变换（discrete Fourier transform, DFT）的关键。

离散傅里叶变换（DFT）

理论

从上面的讨论我们知道，非周期离散信号的频谱是连续函数，也就是说频域是连续的，因此不便于计算机处理。如果在的周期内进行等间隔采样，采样点的个数和的长度相同，将得到的频域离散序列反变换到时域中（IDFS），便得到一个周期信号，它在一个周期内的值正好和相同，或者说是它的周期延拓。

由此定义离散傅里叶变换（discrete Fourier transform, DFT），它是DFS的一个主周期成分。DFT是可以在计算机中实现的傅里叶变换，可用于分析任意有限长序列的频谱。DFT的定义式如下，求和范围从到：

可以发现，和DFS的计算公式是相同的。反变换IDFT的定义式为

同样地，使用简化记号，令，DFT公式对、还可以写成

DFT的物理意义不明显，它实际上是信号的周期延拓的DFS序列，那我们为什么要使用它呢？原因如下：

  1. 的真实频谱是通过DTFT得到的，它是连续的，计算机不方便处理，而DFT得到的是离散的，能够在计算机上实现
  2. 可以证明，通过DFT得到的序列可以完全复原
  3. 的包络线就是在一个周期内的图像，如果的点足够多，包络线非常光滑，可以看成是的图像

第3点可以通过如下代码验证：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
% 信号的长度为N，开始L个值是1，其余为0
L = 4; N = 8;
x = zeros(1, N); x(1:L) = 1;

% 信号的DTFT表达式，L不变时表达式保持不变
w = 0:0.01:2*pi;
xejw = 1 + exp(-1j*w) + exp(-1j*2*w) + exp(-1j*3*w);

% 信号的DFT，magX为幅值，phaX为相角
X = fft(x, N);
magX = abs(X); phaX = rad2deg(angle(X));
n = 0:2/N:2-2/N;

% 实线为DTFT的幅值，圆圈为DFT的幅值
subplot(2, 1, 1); plot(w/pi, abs(xejw)); grid;
hold on; plot(n, magX, 'ro');
xlabel('Frequency in \pi units'); ylabel('Magnitude');
% 实线为DTFT的相角，圆圈为DFT的相角
subplot(2, 1, 2); plot(w/pi, rad2deg(angle(xejw))); grid;
hold on; plot(n, phaX, 'ro');
xlabel('Frequency in \pi units'); ylabel('Phase')

图4是N取8时代码的运行结果，可以看到，将DTFT在到区间等间隔取个点即可得到DFT。

图4 N=8的DFT图像

如果L是定值，那么DTFT的表达式保持不变。N取得比较大时（即在后面补零），DFT的点数更多，包络线会更加接近DTFT的图像，图5是N取64时的结果。

图5 N=64的DFT图像

DFT算法的时间复杂度为，难以用于大规模数据的计算。后来出现了优化算法，即快速傅里叶变换（fast Fourier transform, FFT），它的时间复杂度为，使得DFT成了实用的算法。

实例

某序列，序列长度，求和区间从到，，DFT计算如下：

所以

在MATLAB中，用fft命令计算DFT如下：

1
2
3
4
5
6
>> x = [0, 1, 2, 3]
x =
     0     1     2     3
>> y = fft(x)
y =
   6.0000 + 0.0000i  -2.0000 + 2.0000i  -2.0000 + 0.0000i  -2.0000 - 2.0000i

参考资料

  1. 信号与系统(奥本海姆)
  2. 数字信号处理（MATLAB版）