用MATLAB理解正弦波序列的DFT

前言

在穿戴式嵌入式模拟前端（AFE-Analog Front End）设计中，常常要测量人体组织的阻抗。通常的做法是向人体组织发出一段某一频率的正弦波，测量组织两端的电压与流过组织的电流，绕后分别计算电压与电流的离散傅里叶变换，得到二者的频域特性。最后将电压幅值除以电流幅值的到组织阻抗幅值，电压相位减去电流相位得到组织的阻抗相位。

本文分析了正弦波、余弦波被抽样后的得到的序列的DFT数据特征，以及相位不为零的情况下的DFT数据特征，并用MATLAB工具验证了数据特征，这对阻抗测量的嵌入式软件设计具有一定的指导作用。

已知单位矩形序列的DFT为单位冲击序列的N倍
$$X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(R(n))=N\delta (k)$$

以及时域内的调制对应频域内的搬移：
$$\begin{gathered} 如果：X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(x(n)) \\ 和：y(n)=x(n)W_N^{-ln} \\ 则有：Y(k)=\mathrm{D}\mathrm{F}\mathrm{T}(y(n))=X((k-l){)}_N \end{gathered}$$
其中：$W_N=\mathrm{e}\mathrm{x}\mathrm{p}(\frac{-2\pi i}{N})$

现在想要推出正弦波序列和余弦波序列的离散傅里叶变换

欧拉公式

正弦信号、余弦信号与指数信号的关系可用欧拉公式表示：
$$\mathrm{e}\mathrm{x}\mathrm{p}(ix)=\cos (x)+i\sin (x)$$

$$\mathrm{e}\mathrm{x}\mathrm{p}(\frac{2\pi ni}{N})=\cos (\frac{2\pi n}{N})+i\sin (\frac{2\pi n}{N})$$

正余弦函数的DFT

1. 包含一个正余弦周期、相位为0的序列的DFT

长度为$N$的指数序列的DFT:
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\mathrm{e}\mathrm{x}\mathrm{p}(\frac{2\pi ni}{N}))=N\delta ((k-1){)}_N \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
长度为$N$的余弦序列的DFT：
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\cos (\frac{2\pi n}{N}))=\frac{N}{2}(\delta ((k-1){)}_N+\delta ((k-N+1){)}_N) \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
长度为$N$的正弦序列的DFT
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\sin (\frac{2\pi n}{N}))=\frac{Ni}{2}(\delta ((k-1){)}_N-\delta ((k-N+1){)}_N) \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
以下MATLAB上运行的用m语言写的8个数据、一个周期、相位为0的余弦波序列 x 及其用FFT算法计算的离散傅里叶变换 X，以及它们的波形图 。余弦波序列由8个数据组成, 三张图从左到右依次为 序列x, 傅里叶变换 X 的实部和 X 的虚部。

从图上可以看出，8个数据、一个周期、相位为0的余弦波序列的DFT的特征为：

实部：$X(2)=X(8)=4$，其余为零；
虚部：全为0.
注解：MATLAB的序列号是从 1 开始计算，而不是从 0 开始计算的。

N=8; 
n = 0:1:(N - 1); 
x=cos(2*pi*(n-0)/N); 
figure(1); 
plot(x);
X=fft(x); 
figure(2); 
stem(real(X),'filled'); 
figure(3); 
stem(imag(X),'filled'); 

2. 包含一个正余弦周期、相位不为0的序列的DFT

长度为$N$，相位右移$m$的指数序列的DFT:
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\mathrm{e}\mathrm{x}\mathrm{p}(\frac{2\pi (n-m)i}{N}))=N\delta ((k-1){)}_N\mathrm{e}\mathrm{x}\mathrm{p}(\frac{-2\pi mk}{N}) \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
长度为$N$，相位右移$m$的余弦序列的DFT：
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\cos (\frac{2\pi (n-m)}{N}))=\frac{N}{2}(\delta ((k-1){)}_N+\delta ((k-N+1){)}_N)\cos (\frac{-2\pi mk}{N}) \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
长度为$N$，相位右移$m$的正弦序列的DFT
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\sin (\frac{2\pi (n-m)}{N}))=\frac{Ni}{2}(\delta ((k-1){)}_N-\delta ((k-N+1){)}_N)\sin (\frac{-2\pi mk}{N}) \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
以下MATLAB上运行的用m语言写的8个数据、一个周期、相位右移 1 的余弦波序列 x 及其用FFT算法计算的离散傅里叶变换 X，以及它们的波形图 。余弦波序列由8个数据组成, 三张图从左到右依次为 序列x, 傅里叶变换 X 的实部和 X 的虚部。

从图上可以看出，8个数据、一个周期、相位右移 1 的余弦波序列的DFT的特征为：

实部：$X(2)=X(8)=\frac{8}{2}\cos (-\frac{2\pi }{8})=2.828$，其余为零；
虚部：$X(2)=X(8)=\frac{8}{2}\sin (-\frac{2\pi }{8})=-2.828$，其余为零；
注解：MATLAB的序列号是从 1 开始计算，而不是从 0 开始计算的。

N=8; 
n = 0:1:(N - 1); 
x=cos(2*pi*(n-1)/N); 
figure(1); 
plot(x);
X=fft(x); 
figure(2); 
stem(real(X),'filled'); 
figure(3); 
stem(imag(X),'filled'); 

3. 包含多个正余弦周期、相位为0的序列的DFT

包含 $l$ 个周期、每个周期$N$个数据、长度为$Nl$的指数序列的DFT:
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\mathrm{e}\mathrm{x}\mathrm{p}(\frac{2\pi nli}{N}))=Nl\delta ((k-l){)}_N, \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
包含 $l$ 个周期、每个周期$N$个数据、长度为$Nl$的余弦序列的DFT：
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\cos (\frac{2\pi nl}{N}))=\frac{Nl}{2}(\delta ((k-l){)}_N+\delta ((k-N+l){)}_N) \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
包含 $l$ 个周期、每个周期$N$个数据、长度为$Nl$的正弦序列的DFT
$$\begin{gathered} X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\sin (\frac{2\pi nl}{N}))=\frac{Ni}{2}(\delta ((k-l){)}_N-\delta ((k-N+l){)}_N) \\ n,k=0,1,2,...,Nl-1 \end{gathered}$$
以下MATLAB上运行的用m语言写的16个数据、2个周期、相位为0的余弦波序列 x 及其用FFT算法计算的离散傅里叶变换 X，以及它们的波形图 。每个周期的余弦波由8个数据组成, 三张图从左到右依次为 序列x, 傅里叶变换 X 的实部和 X 的虚部。

从图上可以看出，16个数据、2个周期、相位为0的余弦波序列的DFT的特征为：

实部：$X(3)=X(15)=8$，其余为零；
虚部：全为0.
注解：MATLAB的序列号是从 1 开始计算，而不是从 0 开始计算的。

N=8; 
n = 0:1:(2*N - 1); 
x=cos(2*pi*(n-0)/N); 
figure(1); 
plot(x);
X=fft(x); 
figure(2); 
stem(real(X),'filled'); 
figure(3); 
stem(imag(X),'filled'); 

4. 包含多个正余弦周期、相位不为0的序列的DFT

包含 $l$ 个周期、每个周期$N$个数据、长度为$Nl$、相位右移 $m$的指数序列的DFT:
$$X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\mathrm{e}\mathrm{x}\mathrm{p}(\frac{2\pi (n-m)li}{N}))=Nl\delta ((k-l){)}_N\mathrm{e}\mathrm{x}\mathrm{p}(\frac{-2\pi mk}{N})$$
包含 $l$ 个周期、每个周期$N$个数据、长度为$Nl$、相位右移 $m$的余弦序列的DFT：
$$X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\cos (\frac{2\pi (n-m)l}{N}))=\frac{Nl}{2}(\delta ((k-l){)}_N+\delta ((k-N+l){)}_N)\cos (\frac{-2\pi mk}{N})$$
包含 $l$ 个周期、每个周期$N$个数据、长度为$Nl$、相位右移 $m$的正弦序列的DFT
$$X(k)=\mathrm{D}\mathrm{F}\mathrm{T}(\sin (\frac{2\pi (n-m)l}{N}))=\frac{Nl}{2}(\delta ((k-l){)}_N-\delta ((k-N+l){)}_N)\sin (\frac{-2\pi mk}{N})$$
以下MATLAB上运行的用m语言写的16个数据、2个周期、相位右移 1 的余弦波序列 x 及其用FFT算法计算的离散傅里叶变换 X，以及它们的波形图 。每个周期的余弦波由8个数据组成, 三张图从左到右依次为 序列x, 傅里叶变换 X 的实部和 X 的虚部。

从图上可以看出，16个数据、2个周期、相位右移 1 的余弦波序列的DFT的特征为：

实部：$X(3)=X(15)=\frac{8\ast 2}{2}\cos (-\frac{2\pi }{8})=5.656$，其余为零；
虚部：$X(3)=X(15)=\frac{8\ast 2}{2}\sin (-\frac{2\pi }{8})=-5.656$，其余为零；
注解：MATLAB的序列号是从 1 开始计算，而不是从 0 开始计算的。

N=8; 
n = 0:1:(2*N - 1); 
x=cos(2*pi*(n-1)/N); 
figure(1); 
plot(x);
X=fft(x); 
figure(2); 
stem(real(X),'filled'); 
figure(3); 
stem(imag(X),'filled'); 

结论

正弦波序列的DFT变换后的频率域序列$X(k)$只在某两项上不为零，其余项都为0. 不为零的两项数值相等，而且反映了正弦波序列的全部信息：其幅度代表正弦波的幅度，其相位代表了正弦波的相位。因此，
将采集到的正弦波序列通过DFT变换为频率域的序列，只需要计算非0项的幅度和相位，便可获得正弦波序列的幅度和相位信息。

正弦波的相位信息的精度取决于采样率。在正弦波的一个周期中如果有$N$个采样点，则相位的精度为$\frac{2\pi }{N}$.