离散傅里叶变换

预备知识·

傅里叶级数

傅里叶变换是 以时间 $t$ 为自变量的时域信号 和 以频率 $f$ 为自变量的频域函数 之间的变换关系。

对于周期为 $T_p$ 的连续时间信号 $\tilde{x}(t)$，其频谱 $X(jk\Omega_p)$ 是离散非周期函数

$$\begin{aligned} X(jk\Omega_p) &= \dfrac{1}{T_p}\int_{\tfrac{1}{T_p}}^{\tfrac{1}{T_p}} x(t)e^{-jk\Omega_p t}dt\\ \tilde{x}(t) &= \sum_{k =-\infty}^{+\infty} X(jk\Omega_p)e^{jk\Omega_p t} \end{aligned}$$

其中 $\Omega_p=\dfrac{2\pi}{T_p}$ 为离散频谱两谱线间的角频率间隔，$k$ 为谐波序号。

对于连续的非周期时间信号 $x(t)$，其频谱 $X(j\Omega)$ 是一个连续的非周期函数。满足

$$\begin{aligned} X(j\Omega) &= \int_{-\infty}^{+\infty} x(t)e^{-j\omega t}\\ x(t) &= \dfrac{1}{2\pi}\int_{-\infty}^{+\infty} X(j\Omega)e^{j\Omega t}d \Omega \end{aligned}$$

对于离散非周期信号，其频谱 $\tilde{X}(e^{j\omega})$ 是连续周期信号

$$\begin{aligned} \tilde{X}(e^{j\omega}) &= \sum_{n =-\infty}^{+\infty} x [n] e^{-j\omega n}\\ x [n] &= \dfrac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})e^{j\omega n}d\omega \end{aligned}$$

其中 $\omega$ 是数字角频率，满足 $\omega=\Omega_p T_s$。

可以看出，时间域的周期造成频谱的离散，时间域的非周期造成频谱的连续。

不过，上面的三种傅里叶变换总有一个域是连续的，这不能利用计算机辅助计算。

周期序列的离散傅里叶级数（DFS）·

离散非周期信号的频谱是连续的，将连续的傅里叶变换 $N$ 点采样后，频域离散，时域周期延拓。因此，离散傅里叶级数对周期序列才存在。

设 $\tilde{x}[n]$ 是周期为 $N$ 的一个离散周期序列

$$\tilde{x}[n] = \tilde{x}[n+kN], k\in N$$

由离散时间傅里叶变换（DTFT）可知，该周期序列可以表示为

$$\tilde{x}[n] = \dfrac{1}{N}\sum_{k = 0}^{N-1}\tilde{X}(k)e^{j\tfrac{2\pi}{N}kn}, n = 0,\pm 1, \cdots$$

其中 $\tilde{X}(k)$ 是离散傅里叶级数的系数，

$$\tilde{X}(k) = \sum_{n = 0}^{N-1}\tilde{x}[n] e^{-j\tfrac{2\pi}{N}nk}, k = 0,\pm1,\cdots$$

并且

$$\tilde{X}(k+mN) = \sum_{n = 0}^{N-1}\tilde{x}[n] e^{-j\tfrac{2\pi}{N}n(k+mN)}= \sum_{n = 0}^{N-1}\tilde{x}[n] e^{-j\tfrac{2\pi}{N}nk}=\tilde{X}(k), k = 0,\pm1,\cdots$$

也是周期序列。因此，时域的离散周期序列的离散傅里叶级数（DFS）在频域也是离散的周期序列。

令 $W_N=e^{-j\tfrac{2\pi}{N}}$，离散傅里叶级数可表示为

$$\begin{aligned} \tilde{X}(k) &= DFS [\tilde{x}[n]] = \sum_{n = 0}^{N-1}\tilde{x}[n] e^{-j\tfrac{2\pi}{N}kn} = \sum_{n = 0}^{N-1}\tilde{x}[n] W_{N}^{kn}, k = 0,\pm1,\cdots\\ \tilde{x}[n] &= IDFS [\tilde{X}(k)] = \dfrac{1}{N}\sum_{k = 0}^{N-1}\tilde{X}(k)e^{j\tfrac{2\pi}{N}nk} = \dfrac{1}{N}\sum_{k = 0}^{N-1}\tilde{X}(k)W_{N}^{nk}, n = 0,\pm 1, \cdots \end{aligned}$$

连续周期信号与离散周期序列的对比

分类	基频序列	周期	基频	$k$ 次谐波序列
连续周期	$e^{j\Omega_p t}=e^{j\tfrac{2\pi}{T_p}t}$	$T_p$	$\Omega_p=\dfrac{2\pi}{T_p}$	$e^{jk\tfrac{2\pi}{T_p}t}$
离散周期	$e^{j\omega n}=e^{j\tfrac{2\pi}{N}n}$	$N$	$\omega =\dfrac{2\pi}{N}$	$e^{jk\tfrac{2\pi}{N}n}$

由于 $e^{j(k+\lambda N)\tfrac{2\pi}{N}n}=e^{jk\tfrac{2\pi}{N}n}$，离散傅里叶级数的所有谐波成分中只有 $N$ 个是独立的，因此在展开成离散傅里叶级数时，只能取 $N$ 个独立的谐波分量。

有限长序列的离散傅里叶变换（DFT）·

将长度为 $N$ 的有限长序列 $x[n]$ 看作是周期为 $N$ 的周期序列的一个周期，计算该周期序列的离散傅里叶级数，即可得到有限长序列的离散傅里叶变换。

$$x [n] = \left\{ \begin{aligned} &\tilde{x}[n]&& 0\leq n\leq N-1\\ &0&& else\\ \end{aligned} \right.$$

也可以表示为周期延拓或时窗运算

$$\tilde{x}[n] = x [[n]]_N \Leftrightarrow x [n] = \tilde{x}[n] R_N [n]$$

可以得到有限长序列的离散傅里叶变换为

$$\begin{aligned} X(k) &= DFT [x[n]] = \sum_{n = 0}^{N-1}x [n] W_N^{kn}, 0\leq k\leq N-1\\ x [n] &= IDFT [X(k)] = \frac{1}{N} \sum_{n = 0}^{N-1}X(k)W_{N}^{-nk}, 0\leq n\leq N-1 \end{aligned}$$

有限长序列 $x[n]$ 的 DFT $X(k)$，就是 $x(n)$ 的周期延拓序列 $\tilde{x}[n]$ 的 DFS 系数 $\tilde{X}(k)$ 的主值序列。DFT 的时域与频域都是有限长序列，方便用计算机表示。

频域采样（频域采样点数 $N$ 的选取）

设时域周期序列 $\tilde{x}[n]$ 的周期为 $M$，

$$\tilde{x}[n] = \dfrac{1}{N}\sum_{k = 0}^{N-1}X(k)e^{j\tfrac{2\pi}{N}nk}$$

将 $X(k) = \displaystyle\sum_{m=0}^{M-1}x[m]e^{-j\tfrac{2\pi}{N}km}$ 带入上式，得到

$$\begin{aligned} \tilde{x}[n] &= \dfrac{1}{N}\sum_{k = 0}^{N-1}\left(\sum_{m = 0}^{M-1}x [m] e^{-j\tfrac{2\pi}{N}km}\right)e^{j\tfrac{2\pi}{N}nk}\\ &= \sum_{m = 0}^{M-1}x [m]\left( \dfrac{1}{N}\sum_{k = 0}^{N-1} e^{j\tfrac{2\pi}{N}(n-m)k}\right) \end{aligned}$$

其中，当 $n-m=\lambda N$ 时

$$\dfrac{1}{N}\sum_{k = 0}^{N-1} e^{j\tfrac{2\pi}{N}(n-m)k}= \left\{ \begin{aligned} 1, n-m =\lambda N\\ 0, n-m\not =\lambda N \end{aligned} \right.$$

则

$$\begin{aligned} \tilde{x}[n] &= \sum_{m = 0}^{M-1}x [m]\sum_{\lambda =-\infty}^{+\infty}\delta(n-m-\lambda N)\\ &= \sum_{\lambda =-\infty}^{+\infty}\sum_{m =-\infty}^{+\infty}x [m]\delta(n-m-\lambda N)\\ &= \sum_{\lambda =-\infty}^{+\infty}x [n-\lambda N] \end{aligned}$$

说明，频域 $N$ 点采样等价于时域上以 $N$ 为周期延拓，因此 $N\geq M$ 时，时域不发生混叠。

• 循环移位性质

$$\begin{aligned} x [[n+m]]_NR_N [n] &\leftrightarrow X(k)W_N^{-km}\\ X((k+l))_NR_N(k) &\leftrightarrow x [n] W_N^{lk} \end{aligned}$$

频域移位点数 与对应的 时域频偏频率 的关系：

假设捕获模块的输入信号时经过 $m$ 倍内插的信号，即 $f_s=mR_c$，则在频域进行 $l$ 点循环移位对应时域生成的载波为

$$e^{-j\tfrac{2\pi}{N}nl} = e^{-j\tfrac{2\pi}{m\cdot SSR}(f_st)l} = e^{-j\tfrac{2\pi R_c}{SSR}tl} = e^{-j2\pi R_s tl}$$

其中 $SSR$ 表示扩频比，$m$ 为内插系数，$f_s$ 为采样频率，$R_c$ 为码片速率，$R_s$ 为符号速率。

可以得到结论：频域循环移位 $l$ 点对应时域乘以频率为为 $l\cdot R_s$ 的载波.

• 共轭对称性（类似于偶对称和奇对称）

在 DFT 中，涉及的序列 $x[n]$ 和 $X(k)$ 均为有限长序列，主值区间为 $[0,N-1]$，这里的对称性指的是关于 $\dfrac{N}{2}$ 的圆周对称，长度为 $N$ 的有限长序列 $x[n]$ 的圆周共轭对称分量 $x_{ep}[n]$ 和圆周共轭反对称分量 $x_{op}[n]$ 分别定义为

$$\begin{aligned} x_{ep}[n] &= \dfrac{1}{2}(x [n] + x^*[N-n])\\ x_{op}[n] &= \dfrac{1}{2}(x [n]-x^*[N-n])\\ \end{aligned}$$

满足

$$\begin{aligned} x_{ep}[n] &= x^*_{ep}[N-n]\\ x_{op}[n] &= -x^*_{op}[N-n]\\ \end{aligned}$$

因此，任何有限长序列都可分解为圆周共轭对称分量 $x_{ep}[n]$ 和圆周共轭反对称分量 $x_{op}[n]$ 的和，

$$x [n] = x_{ep}[n] + x_{op}[n], 0\leq n\leq N-1$$

其中

$$\begin{aligned} DFT [x_{ep}[n]] &= \dfrac{1}{2}[X(k) + X^*(k)] = \Re[X(k)]\\ DFT [x_{op}[n]] &= \dfrac{1}{2}[X(k) - X^*(k)] = \Im[X(k)]\\ \end{aligned}$$

$$DFT[x^*[n]] = X^*(N-k)$$

则

$$\begin{aligned} X(k) &= DFT [x_{ep}[n]] + DFT [x_{op}[n]]\\ &=\Re [X(k)] + \Im [X(k)]\\ &= X_R(k) + jX_I(k) \end{aligned}$$

同理，将 $x[n]$ 表示为 $x_r[n] + jx_i[n]$，其中

$$\begin{aligned} x_r [n] &= \dfrac{1}{2}(x [n] + x^*[n])\\ jx_i [n] &= \dfrac{1}{2}(x [n] - x^*[n]) \end{aligned}$$

则

$$\begin{aligned} DFT [x_r[n]] &= \dfrac{1}{2}[X(k) + X^*(N-k)] = X_{ep}(k)\\ DFT [jx_i[n]] &= \dfrac{1}{2}[X(k) - X^*(N-k)] = X_{op}(k) \end{aligned}$$

$$DFT [x[n]] = X_{ep}(k) + X_{op}(k)$$

DFT 与 DTFT、ZT 的关系·

已知 $x[n]$ 为 $N$ 点有限长序列，则 $x[n]$ 的 $z$ 变换为

$$X(z) = \sum_{n = 0}^{N-1} x [n] z^{-n}$$

DTFT 为

$$X(e^{j\omega}) = \sum_{n = 0}^{N-1}x [n] e^{-j\omega n}|_{z = e^{j\omega}}$$

DFT 为

$$\begin{aligned} X(k) &= \sum_{n = 0}^{N-1}x [n] W_{N}^{kn} = \sum_{n = 0}^{N-1}x [n] e^{-j\tfrac{2\pi}{N}kn}\\ &= X(z)|_{z = e^{j\frac{2\pi}{N}k}} = X(e^{j\omega})|_{\omega = \frac{2\pi}{N}k} \end{aligned}$$

有限长序列 $x[n]$ 的 DTFT 是单位圆上的 $z$ 变换，$\omega$ 是连续的，体现出 DTFT 是以 $2\pi$ 为周期的连续谱。

有限长序列 $x[n]$ 的 DFT 是对 $X(z)$ 在单位圆上进行 $N$ 点的等间隔采样，也是 DTFT 一个周期 $[0,2\pi)$ 上的 $N$ 点等间隔采样。

如果频域采样满足频域采样定理（频域采样点数 $N$ 不小于时域序列长度 $M$），就可以用频域采样值恢复序列的 $z$ 变换，

$$\begin{aligned} X(z) &=\sum_{n = 0}^{N-1}x [n] z^{-n}\\ &=\sum_{n = 0}^{N-1}\left(\dfrac{1}{N}\sum_{k = 0}^{N-1}X(k)W_N^{-nk}\right)z^{-n}\\ &=\dfrac{1}{N}\sum_{k = 0}^{N-1}X(k)\left(\sum_{n = 0}^{N-1}W_{N}^{-nk}z^{-n}\right)\\ &=\dfrac{1}{N}\sum_{k = 0}^{N-1}X(k)\dfrac{1-W^{-k}z^{-1}}{1-W_N^{-k}z^{-N}}\\ &= \dfrac{1}{N}\sum_{k = 0}^{N-1}X(k)\dfrac{1-z^{-N}}{1-W_{N}^{-k}z^{-1}}\\ \end{aligned}$$

令 $\varPhi_k(z)=\dfrac{1}{N}\dfrac{1-z^{-N}}{1-W_{N}^{-k}z^{-1}}$，$X(z)$ 可以表示为

$$X(z) = \sum_{k = 0}^{N-1}X(k)\varPhi_k(z)$$

函数 $\varPhi_k(z)$ 有一个极点 $z=W_N^{-k}$，$N$ 个零点 $z=e^{j\tfrac{2\pi}{N}r},r=0,1,\cdots,N-1, r\not=k$，因此内插函数 $\varPhi_k(z)$ 在除了 $k$ 的其它 $N-1$ 个采样点均为 $0$，

同理，

$$\begin{aligned} X(e^{j\omega}) &=\dfrac{1}{N}\sum_{k = 0}^{N-1}X(k)\dfrac{1-e^{-j\omega N}}{1-e^{j\tfrac{2\pi}{N}k}e^{-j\omega}}\\ &= \sum_{k = 0}^{N-1}X(k)\varPhi(\omega - \dfrac{2\pi}{N}k)\\ \end{aligned}$$

其中 $\varPhi(\omega)=\dfrac{1}{N}\dfrac{\sin\tfrac{\omega N}{2}}{\sin\tfrac{\omega}{2}}e^{-j\tfrac{N-1}{2}\omega}$ 为内插函数，满足

$$\varPhi(\omega-\dfrac{2\pi}{N}k) = \left\{ \begin{aligned} &1,&\omega =\dfrac{2\pi}{N}k\\ &0,&\omega =\dfrac{2\pi}{N}i, i\not = k \end{aligned} \right.$$

也就是每个采样点上的 $X(e^{j\omega})$ 的值等于 $X(k)$，即 $X(e^{j\omega})|_{\omega=\tfrac{2\pi}{N}k}=X(k),k=0,1,\cdots,N-1$，各个采样点之间的值由各采样点的加权内插函数 $X(k)\varPhi(\omega - \dfrac{2\pi}{N}k)$ 在所求 $\omega$ 点上的值叠加得到。

由内插函数求频率响应

频域内插（与 $sinc$ 卷积），等价于时域截断（加矩形窗）

线性卷积与线性相关·

周期卷积

设 $\tilde{x}[n]$ 与 $\tilde{h}[n]$ 都是以 $N$ 为周期的序列，DFS 分别为 $\tilde{X}(k)$ 和 $\tilde{H}(k)$，若

$$\tilde{Y}(k) = \tilde{X}(k)\tilde{H}(k)$$

则

$$\tilde{y}[n] = IDFS [\tilde{Y}(k)] =\sum_{m = 0}^{N-1}\tilde{x}[m]\tilde{h}[n-m]$$

圆周卷积（循环卷积）

设当 $x[n]$ 与 $h[n]$ 都是长度为 $N$ 的有限长序列，DFT 分别为 $X(k)$ 和 $H(k)$，若

$$Y(k) = X(k)H(k)$$

则

$$\begin{aligned} y [n] = IDFT [Y(k)] &= \sum_{m = 0}^{N-1}x [m] h [[n-m]]_NR_N [n]\\ &= \sum_{m = 0}^{N-1}h [m] x [[n-m]]_NR_N [n] \end{aligned}$$

记作 $y[n]=x[n]\circledast h[n]$.

当信号通过线性时不变系统时，系统的输出 $y[n]$ 时输入 $x[n]$ 与单位脉冲响应 $h[n]$ 的线性卷积，即 $y[n]=x[n]*h[n]$，当 $x[n]$ 和 $h[n]$ 均为有限长序列时，可以考虑用圆周卷积代替线性卷积，下面分析代替的条件：

设序列 $x_1[n]$ 和 $x_2[n]$ 分别是长度为 $N_1$ 和 $N_2$ 的有限长序列，设 $y_l[n]$ 是两者的线性卷积

$$y_l [n] = \sum_{m =-\infty}^{+\infty} x_1 [m] x_2 [n-m] =\sum_{m = 0}^{N_1-1} x_1 [m] x_2 [n-m]$$

则 $y[n]$ 是长度为 $N_1+N_2-1$ 的有限长序列。

设 $y_c[n]$ 是两者的 $L$ 点圆周卷积

$$x_1 [n] = \left\{ \begin{aligned} &x_1 [n], &0\leq n \leq N_1-1\\ &0, &N_1\leq n\leq L-1 \end{aligned} \right.$$

$$x_2 [n] = \left\{ \begin{aligned} &x_2 [n], &0\leq n \leq N_2-1\\ &0, &N_2\leq n\leq L-1 \end{aligned} \right.$$

$$\begin{aligned} y_c [n] &= \left(\sum_{m = 0}^{L-1}x_1 [m] x_2 [[n-m]]_L\right)R_L [n]\\ &= \left(\sum_{m = 0}^{L-1}x_1 [m]\sum_{r =-\infty}^{+\infty}x_2 [n+rL-m]\right)R_L [n]\\ &=\left(\sum_{r =-\infty}^{+\infty}y_l [n+rL]\right)R_L [n] \end{aligned}$$

因此，$L$ 点圆周卷积 $y_c[n]$ 是线性卷积 $y_l[n]$ 以 $L$ 为周期延拓序列的主值序列，且不失真的条件是

$$L > N_1 + N_2 + 1$$

圆周卷积代替线性卷积

谱分析·

对于时域连续的非周期信号 $x(t)$，其 CTFT 为

$$X(jf) = \int_{-\infty}^{+\infty}x(t)e^{-j2\pi f t}dt$$

频谱泄漏·

理论分析·

实际要把观测的信号 $x[n]$ 限制在一定的时间间隔内，时域截断数学上表示为无限长时间序列乘以窗函数，频域上是两者频谱的卷积，由于窗函数不能无限宽，频谱不是冲激函数，因此，时域截断必然会造成频谱展宽（拖尾），造成频谱泄漏。

以正弦信号 $x(t)=e^{j\Omega t}$ 为例，先对其采样，采样频率为 $f_s$，

$$x [n] = e^{j\Omega n T_s} = e^{j\omega n}, \omega = \Omega T_s$$

做 DFT，得到

$$X(k) = \sum_{n = 0}^{N-1} x [n] W_N^{kn} = \sum_{n = 0}^{N-1} e^{j(\omega-\tfrac{2\pi}{N}k)n}= \left\{ \begin{aligned} N, \omega = \dfrac{2\pi}{N}k\\ 0, \omega \not = \dfrac{2\pi}{N}k\\ \end{aligned} \right.$$

是 $x(t)$ 的频谱以 $\dfrac{2\pi}{N}$ 为周期进行延拓的结果。

$$\omega = \Omega T_{smp} = \dfrac{2\pi}{N}k$$

则 $\Omega = \dfrac{2\pi k}{NT_{smp}}$，其中 $NT_{smp}$ 是样本的长度，并且，对于正弦信号，$\Omega=\dfrac{2\pi}{T_S}$，则有

$$NT_{smp} = k T_S$$

因此，当样本的长度为信号周期的整数倍时，DFT 能正确分析频谱。否则，会因为对信号的不当截断造成频谱泄漏。

MATLAB 仿真·

close all; clear; clc;

tic;

% FIXME System Parameters
fs = 160e3;
t_total = 1;

% Taking single tone signal as an example
fc = 20e3;
t = 0 : 1/fs : t_total - 1/fs;
s = cos(2 *pi* fc * t);

% FIXME Hamming Window Function (or Bartlett, Hanning, Blackman et al.)
N = 2048;  % FIXME The number of FFT points
hamming_window = hamming(length(s));
s_truncated = s .* hamming_window';

f = (-fs/2 : fs/N : fs/2 - fs/N);
S = 10 * log10(abs(fftshift(fft(s, N))) /N);
Hamming_window = 10 * log10(abs(fftshift(fft(hamming_window, N))));
S_truncated  = 10 * log10(abs(fftshift(fft(s_truncated , N))));


figure;
subplot(3, 2, 1);
plot(t, s);
xlabel("Time (s)"); ylabel("Amplitude");
title("The Original Signal"); grid on;

subplot(3, 2, 2);
plot(f, S);
xlabel("Frequency (Hz)"); ylabel("Amplitude (dB)");
title("Spectrum of the Original Signal"); grid on;

subplot(3, 2, 3);
plot(t, hamming_window);
xlabel("Time (s)"); ylabel("Amplitude");
title("Hamming Window"); grid on;

subplot(3, 2, 4);
plot(f, Hamming_window);
xlabel("Frequency (Hz)"); ylabel("Amplitude (dB)");
title("Spectrum of Hamming Window"); grid on;


subplot(3, 2, 5);
plot(t, s_truncated);
xlabel("Time (s)"); ylabel("Amplitude");
title("The Truncated Signal"); grid on;

subplot(3, 2, 6);
plot(f, S_truncated);
xlabel("Frequency (Hz)"); ylabel("Amplitude (dB)");
title("Spectrum of the Truncated Signal"); grid on;

sgtitle("Simulation of Spectrum Leakage");

toc;

栅栏效应·

理论分析·

$N$ 点 DFT 是将 $2\pi$ 周期分成 $N$ 份，在 $\dfrac{2k\pi}{N},k=0,1,\cdots,N-1$ 这几个离散的 栅栏点 上观察一个周期。

考虑一个 $N_1$ 点的时域序列 $\{x[0],x[1],\cdots,x[N_1-1]\}$，其 DFT 为

$$X_1 [k] = \sum_{n = 0}^{N_1-1}x [n] W_{N_1}^{nk}$$

其中 $W=e^{-j2\pi}$，对该时域序列的末尾补 $N_2$ 个 $0$，得到新的时域序列

$$\{x [0], x [1],\cdots, x [N_1+N_2-1]\}$$

其中 $x[N_1]=x[N_1+1]=\cdots=x[N_1+N_2-1]=0$，补零后的 DFT 为

$$X_2 [k] = \sum_{n = 0}^{N_1+N_2-1}x [n] W_{N_1+N_2}^{nk} =\sum_{n = 0}^{N_1-1}x [n] W_{N_1+N_2}^{nk}$$

可以发现 $X_1[k]$ 与 $X_2[k]$ 仅是点数不同，即频率的分辨率不同。

当 $N_2=mN_1$ 时，有

$$X_2 [(m+1)k] = \sum_{n = 0}^{N_1-1}x [n] W_{(m+1)N_1}^{n(m+1)k} =\sum_{n = 0}^{N_1-1}x [n] W_{N_1}^{nk}= X_1 [k]$$

表征补零前后的频谱谱线有相同的频点，在该频点处原信号的幅值被保留，如果补零前后没有相同的频点，则原频点只能由补零后的其它频点合成，能量也就被分散到那些频点上，发生 频谱泄漏。

末尾补零本质上是在增加频率的分辨率，当用 FFT 观察一个离散信号的频谱时，实际是在观察一个连续谱的一些离散点，就像是透过 栅栏 观察频谱。

MATLAB 仿真·

close all; clear; clc;

tic;

% FIXME System Parameters
fs = 1e6;
t_total = 1;

% Taking the dual-tone signal as an example
fc1 = 5e3; fc2 = 6e3;  % FIXME
t = 0 : 1/fs : t_total - 1/fs;
s = cos(2 * pi * fc1 * t) + cos(2 * pi * fc2 * t);
f = (-fs/2 : fs/length(s) : fs/2 - fs/length(s));
S = 10 * log10(abs(fftshift(fft(s))));

figure;
subplot(2, 1, 1);
plot(t, s);
xlabel("Time (s)"); ylabel("Amplitude");
title("The Original Signal"); grid on;

subplot(2, 1, 2);
plot(f, S);
xlabel("Frequency (Hz)"); ylabel("Amplitude (dB)");
title("Spectrum of the Original Signal"); grid on;


N = 2048;  % FIXME The number of FFT points
f_2048 = (-fs/2 : fs/N : fs/2 - fs/N);
S_2048 = 10 * log10(abs(fftshift(fft(s, N))));

N = 1024;
f_1024 = (-fs/2 : fs/N : fs/2 - fs/N);
S_1024 = 10 * log10(abs(fftshift(fft(s, N))));

% After padding zero at the end
N = 2048;
f_zero_padding = (-fs/2 : fs/N : fs/2 - fs/N);
s = [s(1:1024) zeros(1, 1024)];
S_zero_padding = 10 * log10(abs(fftshift(fft(s, N))));

figure;
subplot(3, 1, 1);
plot(f_2048, S_2048);
xlabel("Frequency (Hz)"); ylabel("Amplitude (dB)");
title("Spectrum of the Original Signal ((FFT-2048)"); grid on;

subplot(3, 1, 2);
plot(f_1024, S_1024);
xlabel("Frequency (Hz)"); ylabel("Amplitude (dB)");
title("Spectrum of the Original Signal (FFT-1024)"); grid on;

subplot(3, 1, 3);
plot(f_zero_padding, S_zero_padding);
xlabel("Frequency (Hz)"); ylabel("Amplitude (dB)");
title("Spectrum of the Zero-padding Signal (FFT-2048)"); grid on;

sgtitle("Simulation of Fence Effect");

toc;

内插与抽取·

在原有的离散数据间插入 $U-1$ 个数据，实现对数据的上采样，具体实现流程为 插值 和 滤波。

假设原数据为 $x[n] = \{x[0], x[1], \cdots, x[n-1], x[n]\}$，经过 $U(U=4)$ 倍上采样后得到 $x'[m]=\{x[0], 0, 0, 0, x[1], 0, 0, 0, x[2], \cdots, x[n-1],0,0,0, x[n]\}$.

原数据的 DFT 为

$$X(k) = \sum_{n = 0}^{N-1} x [n] W_{N}^{nk}, 0\leq k\leq N-1$$

经过上采样的数据 $x'[m]$ 的 DFT 为

$$\begin{aligned} X'(k') &= \sum_{n = 0}^{UN-1} x'[Un] W_{UN}^{Unk'}, 0\leq k'\leq UN-1\\ &= \sum_{n = 0}^{UN-1} x'[Un] W_{N}^{nk'}, 0\leq k'\leq UN-1\\ &= \sum_{n = 0}^{N-1} x [n] W_{N}^{nk'}, 0\leq k'\leq UN-1\\ \end{aligned}$$

由此可知，$X'(k')$ 是由原频谱 $X(k)$ 复制 $U$ 次的结果。