离散信号的傅里叶分析

基础知识

定义

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

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

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

\[ \begin{equation} x(n) = e^{j \omega_0 n} \tag{1} \end{equation} \]

它可以看成是连续信号在整数点上的取值,和连续信号类似,但又有区别。可以使用如下的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]')

振荡速率

在连续信号情况下,每个\(\omega_0\)对应的是不同的信号,但是离散情况下就不同了,由于

\[ \begin{equation} e^{j (\omega_0+2\pi) n} = e^{j \omega_0 n} \cdot e^{j 2\pi n} = e^{j \omega_0 n} \tag{2} \end{equation} \]

所以,离散周期复指数信号在频率\(\omega_0+2\pi\)与频率\(\omega_0\)时是完全一样的。

所以,考虑这种离散时间周期复指数信号时,仅仅需要在某个\(2\pi\)间隔内选择\(\omega_0\)即可,一般取\(0 \le \omega_0 < 2\pi\)\(-\pi \le \omega_0 < \pi\)。由上面的讨论可知,\(\omega_0\)数值增加时,\(e^{j \omega_0 n}\)的振荡速率并不会不断增加。\(\omega_0\)\(0\)增加到\(\pi\)时,振荡速率越来越快,从\(\pi\)变化到\(2\pi\)时,振荡速率会下降。下图是\(\omega_0\)取不同值时对应正弦信号的图像:

图1 不同角频率下的图像

所以,离散复指数信号的低频部分位于\(\omega_0\)\(\pi\)的偶数倍值附近,高频部分则位于\(\omega_0\)\(\pi\)的奇数倍值附近。

周期公式

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

\[ \begin{equation} e^{j \omega_0 (n + N)} = e^{j \omega_0 N} \cdot e^{j \omega_0 n} = e^{j \omega_0 n} \tag{3} \end{equation} \]

所以\(e^{j \omega_0 N} = 1\),必须有个整数\(m\),满足

\[ \begin{equation} \omega_0 N = 2 \pi m \tag{4} \end{equation} \]

\[ \begin{equation} \frac{\omega_0}{2 \pi} = \frac{m}{N} \tag{5} \end{equation} \]

所以,只有\(\omega_0 / 2 \pi\)为有理数时,\(e^{j \omega_0 n}\)才是周期信号,否则就不是周期的。这是由于离散时间信号仅能在自变量的整数值上有定义造成的。

基波频率

\[ \begin{equation} \frac{2\pi}{N} = \frac{\omega_0}{m} \tag{6} \end{equation} \]

那么基波周期\(N\)也可以写成

\[ \begin{equation} N = m \left( \frac{2 \pi}{\omega_0} \right) \tag{7} \end{equation} \]

基于上面的原因,通常给定周期复指数信号的基波周期\(N\),然后求得频率\(\omega_0=2\pi/N\)(此时\(m=1\)),得到它的表达式\(e^{j \omega_0 n}\)

谐波信号

再讨论成谐波关系的周期离散时间复指数信号。下面的信号具有公共周期\(N\),频率都是基波周期\(2\pi/N\)的整数倍:

\[ \begin{equation} \phi_k(n) = e^{j k (2\pi/N) n}, \, k = 0, \pm 1, \cdots \tag{8} \end{equation} \]

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

\[ \begin{equation} \phi_{k+N}(n) = e^{j (k+N) (2\pi/N) n} = e^{j k (2\pi/N) n} \cdot e^{j 2\pi n} = \phi_k(n) \tag{9} \end{equation} \]

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

\[ \begin{equation} \phi_k(n) = \phi_{k+rN}(n) \tag{10} \end{equation} \]

离散傅里叶级数(DFS)

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

\[ \begin{equation} \tilde{x}(n) = \frac{1}{N} \sum_{k=\left \langle N \right \rangle} \tilde{X}(k) e^{j k \omega_0 n} = \frac{1}{N} \sum_{k=\left \langle N \right \rangle} \tilde{X}(k) e^{j k (2\pi/N) n} \tag{11} \end{equation} \]

求和限\(k=\left \langle N \right \rangle\)表示\(k\)\(N\)个相继整数的区间上变化,下同。系数\(\tilde{X}(k)\)的计算公式为

\[ \begin{equation} \tilde{X}(k) = \sum_{n=\left \langle N \right \rangle} \tilde{x}(n) e^{-j k \omega_0 n} = \sum_{n=\left \langle N \right \rangle} \tilde{x}(n) e^{-j k (2\pi/N) n} \tag{12} \end{equation} \]

很多教材求和限取\(0\)\(N-1\)区间,此时\((11)\)\((12)\)可写成

\[ \begin{equation} \tilde{x}(n) = \frac{1}{N} \sum_{k=0}^{N-1} \tilde{X}(k) e^{j k \omega_0 n} = \frac{1}{N} \sum_{k=0}^{N-1} \tilde{X}(k) e^{j k (2\pi/N) n} \tag{13} \end{equation} \]

\[ \begin{equation} \tilde{X}(k) = \sum_{n=0}^{N-1} \tilde{x}(n) e^{-j k \omega_0 n} = \sum_{n=0}^{N-1} \tilde{x}(n) e^{-j k (2\pi/N) n} \tag{14} \end{equation} \]

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

如果令\(W_N=e^{-j\frac{2\pi}{N}}\),则DFS公式对\((11)\)\((12)\)还可以写成如下形式

\[ \begin{equation} \tilde{X}(k) = \sum_{n=\left \langle N \right \rangle} \tilde{x}(n) W_N^{nk} \tag{15} \end{equation} \]

\[ \begin{equation} \tilde{x}(n) = \frac{1}{N} \sum_{k=\left \langle N \right \rangle} \tilde{X}(k) W_N^{-nk} \tag{16} \end{equation} \]

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

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

来看一个例子。现有一个离散周期信号\(\tilde{x}(n)\),其周期为\(N\),某个周期片段如下:从\(-N_1\)\(N_1\)的值为\(1\),其余为\(0\),即一个周期\(N\)内,有\(2N_1+1\)个时间点的值为\(1\)。如下的代码取周期信号\(\tilde{x}(n)\)在一个周期的片段\(x(n)\),使用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 周期离散信号及其频谱

读者可以尝试增大\(N\)的值,看频谱图如何变化。

离散时间傅里叶变换(DTFS)

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

\[ \begin{equation} x(n) = \frac{1}{2\pi} \int_{2\pi} X(e^{j\omega}) e^{j \omega n} \mathrm{d}\omega \tag{17} \end{equation} \]

\[ \begin{equation} X(e^{j\omega}) = \sum_{n=-\infty}^{+\infty} x(n) e^{-j \omega n} \tag{18} \end{equation} \]

从DTFT得到的频谱\(X(e^{j\omega})\)是连续函数,由于\(\omega_0\)相差\(2\pi\)的两个复指数信号是一样的,所以\(X(e^{j\omega})\)是周期函数,周期为\(2\pi\)。也就是说,离散时间非周期信号的频谱是连续的、周期的。频谱中靠近偶数倍\(\pi\)\(\omega\),对应的是低频率的信号;而靠近奇数倍\(\pi\)\(\omega\),对应的则是高频率的信号。

举个例子,取上图中周期信号的一个片段:\(n\)\(-2\)取到\(17\),一共取\(N=20\)个点。那么该信号的DTFT计算如下:

\[ X(e^{j\omega}) = e^{j 2\omega} + e^{-j 2\omega} + e^{j \omega} + e^{-j \omega} + 1 = 2 \cos 2\omega + 2 \cos \omega + 1 \]

一般来说\(X(e^{j\omega})\)是复函数,但此例中得到的是实函数,可直接画出它的图像,注意横轴标度是以\(\pi\)为单位。

图3 DTFT是周期连续函数

DTFS与DFS的关系

与图2的频谱图作对比可发现,频谱图的包络线就是本例中的周期连续函数\(X(e^{j\omega})\)。实际上,可以证明,将任意非周期函数\(x(n)\)的离散时间傅里叶变换\(X(e^{j\omega})\)等间隔采样,可得到\(x(n)\)的周期延拓\(\tilde{x}(n)\)的傅里叶级数系数\(\tilde{X}(k)\),采样间隔为\(2\pi/N\),即

\[ \begin{equation} \tilde{X}(k) = X(e^{j\omega}) \mid_{\omega = \frac{2\pi}{N} k} \tag{19} \end{equation} \]

\(\omega_1=2\pi / N\)\(\omega_k = k \omega_1\),上述关系还可以记为

\[ \begin{equation} \tilde{X}(k) = X(e^{j \omega_k}) = X(e^{j k \omega_1}) \tag{20} \end{equation} \]

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

离散傅里叶变换(DFT)

理论

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

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

\[ \begin{equation} X(k) = \sum_{n=0}^{N-1} x(n) e^{-j k \omega_0 n} = \sum_{n=0}^{N-1} x(n) e^{-j k (2\pi/N) n} \tag{21} \end{equation} \]

可以发现,\((21)\)和DFS的计算公式\((14)\)是相同的。反变换IDFT的定义式为

\[ \begin{equation} x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) e^{j k \omega_0 n} = \frac{1}{N} \sum_{k=0}^{N-1} X(k) e^{j k (2\pi/N) n} \tag{22} \end{equation} \]

同样地,使用简化记号,令\(W_N=e^{-j\frac{2\pi}{N}}\),DFT公式对\((21)\)\((22)\)还可以写成

\[ \begin{equation} X(k) = \sum_{n=0}^{N-1} x(n) W_N^{nk} \tag{23} \end{equation} \]

\[ \begin{equation} x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) W_N^{-nk} \tag{24} \end{equation} \]

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

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

第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在\(0\)\(2\pi\)区间等间隔取\(N\)个点即可得到DFT。

图4 N=8的DFT图像

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

图5 N=64的DFT图像

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

实例

某序列\(x(n)=\{0,1,2,3\}\),序列长度\(N=4\),求和区间从\(0\)\(3\)\(W_4=e^{-j \frac{2\pi}{4}} = -j\),DFT计算如下:

\[ X(k)=\sum_{n=0}^{3} x(n)W_4^{nk}, \, k=0,1,2,3 \]

所以

\[ \begin{align*} X(0) & = \sum_{n=0}^{3} x(n)W_4^{0 \cdot n} = \sum_{n=0}^{3} x(n) = x(0) + x(1) + x(2) + x(3) = 6 \\\\ X(1) & = \sum_{n=0}^{3} x(n)W_4^{n} = \sum_{n=0}^{3} x(n)(-j)^n = -2+2j \\\\ X(2) & = \sum_{n=0}^{3} x(n)W_4^{2n} = \sum_{n=0}^{3} x(n)(-j)^{2n} = -2 \\\\ X(3) & = \sum_{n=0}^{3} x(n)W_4^{3n} = \sum_{n=0}^{3} x(n)(-j)^{3n} = -2-2j \\\\ \end{align*} \]

在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版)