本节主要介绍了数字信号处理中的一些频域上的处理，包括低通带通采样定理，DTFT、DFS、DFT。

带通采样定理

采样是频谱的左右搬移（以采样频率$f_s$左右搬移）
那么对于一个$f_L\sim f_H$的窄带信号采样，有：
- 只有当频谱不发生交叠时，才能完全恢复出时域信号
- 那么需满足如下条件：
  $$
  f_L\ge -f_L+(k-1)f_s
  $$
  
  $$
  f_H\le -f_H+kf_s
  $$
- 整合上述两式可得到：
  $$
  \frac{2f_H}{k}\le f_s \le \frac{2f_L}{k-1}
  $$
  
  $$
  1\le k\le int[\frac{f_H}{B}],B=f_L-f_H
  $$

用一个matlab示例程序来演示带通采样频率的取值

clc;                      %清屏
clear all;                %清变量

fL=60;
fH=69;
B=fH-fL;
k=1:(fH/B);
fs_min=2*fH./k
fs_max=2*fL./(k-1)

结果如下：

奈奎斯特采样定理

采样频率$f_s\ge 2f_H$（信号最高频率为$f_H$）

用一个matlab程序示例来理解

clc;                      %清屏
clear all;                %清变量

N=70;                     %数据长度
fs=64000;                 %采样频率
f0=67000;
f1=3000;

t=0:(N-1);
t=t/fs;                   %时间序列

s0=sin(2*pi*f0*t);
s1=sin(2*pi*f1*t);

%绘图
subplot(211);plot(t,s0); title('f0=67kHz');
subplot(212);plot(t,s1); title('f1=3kHz');

结果如下：
可见采样后67KHz的频率和3KHz的频率的时域信号一样
其原因在于，采样是频谱以$f_s$的线性搬移，搬移后的频率会通过一个低通滤波器，也就是只留下频谱在$\pm \frac{f_s}2$之间的频率，即：
$$
f_{out}=\pm k\times f_s \pm f_i\
f_{out}\le \frac{f_s}2
$$
所以67KHz会有一个搬移到3KHz（67KHz-64KHz），3KHz在32KHz内，故采样滤波后会输出一个3KHz的正弦波

时域离散傅里叶变换

离散时间傅里叶变换的公式：
$$
x(e^{jw})=\sum_{n=-\infty}^{+\infty}x(n)e^{-jwn}
$$

$$
逆变换:x(n)=\frac1{2\pi}\int_{2\pi}x(e^{jw})e^{jwn}dw
$$
$x(e^{jw})$是角频率以$2\pi$为周期的函数
理想低通滤波器的iDTFT
- 逆变换的公式推导：
  $$
  \begin{aligned}
  h(n)&=\frac1{2\pi}\int_{2\pi}x(e^{jw})e^{jwn}dw\
  &=\frac1{2\pi}\int_{w_c}^{w_c}e^{jwn}dw\
  &=\frac1{2\pi jn}[e^{jw_cn}-e^{-jw_cn}]\
  &=\frac{sin(w_cn)}{\pi n}
  \end{aligned}
  $$
- 由Matlab可以画出理想低通滤波器的时域波形
  1
  2
  3
  4
  5
  6
  7
  8
  9
  clc; %清屏
  clear all; %清变量
  
  wc=0.5;
  n=-200:200;
  hn=sin(wc*n)./(pi*n);
  stem(n,hn)
  xlabel('n');
  ylabel('hn');
- 结果如下：

时域离散傅里叶级数

离散傅里叶级数的公式：
$$
x(k)=\frac1N\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}Nkn}\quad k=0\sim N-1
$$

$$
x(n)=\sum_{k=0}^{N-1}x(k)e^{j\frac{2\pi}{N}kn}
$$

更详细的解释见我的另一篇文章：信号处理中的小知识点

离散傅里叶变换

离散傅里叶变换的公式：
$$
x(k)=\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi k}Nn} \quad k=0\sim N-1
$$

$$
逆变换:x(n)=\frac1N\sum_{k=0}^{N-1}x(k)e^{j\frac{2\pi k}Nn}
$$
DFT一定是关于$\frac N2$共轭对称的，也就是$x(k_0)=x(N-k_0)$共轭

1.MATLAB辅助理解DFT

通过Matlab来理解DFT的计算过程：

clc;                      %清屏
clear all;                %清变量

N=16;                     %矩形序列的周期
L=4;                      %矩形序列的长度
rn=ones(1,L);             %矩形的长度
xn=[rn,zeros(1,N-L)];
n=0:N-1;

xk=fft(rn,N);             %离散谱

w=0:0.1:2*pi;
xw=(1-exp(-j*w*L))./(1-exp(-j*w)); %连续谱

subplot(311)
stem(n,xn);title('x(n)');grid on;
subplot(312);
stem(n,abs(xk));title('X(K)'); grid on;

subplot(313)
plot(w,abs(xw));title('X(jw)');grid on;

结果如下：

借助matlab讨论N点不同时对幅频特性的影响

不同序列信号，参与DFT计算的N不同时，DFT可以一样，也可以不一样

clc;clear all;            %清屏清变量

n1=0:3; n2=4:7;
n=[n1,n2];

x2n1=n1+1;
x2n2=8-n2;
x2n=[x2n1,x2n2];

x3n1=4-n1;
x3n2=n2-3;
x3n=[x3n1,x3n2];

subplot(321);stem(x2n);title('xn2');
subplot(322);stem(x3n);title('xn3');
subplot(323);stem(abs(fft(x2n,8)));title('DFT(x2n,8)');
subplot(324);stem(abs(fft(x3n,8)));title('DFT(x3n,8)');
subplot(325);stem(abs(fft(x2n,16)));title('DFT(x2n,16)');
subplot(326);stem(abs(fft(x3n,16)));title('DFT(x3n,16)');

figure(2)
subplot(221);stem([x2n x2n x2n x2n]);title('xn2\_8\_repeat4');
subplot(222);stem([x3n x3n x3n x3n]);title('xn3\_8\_repeat4');
x2n=[x2n,zeros(1,8)];
x3n=[x3n,zeros(1,8)];
subplot(223);stem([x2n x2n x2n x2n]);title('xn2\_16\_repeat4');
subplot(224);stem([x3n x3n x3n x3n]);title('xn3\_16\_repeat4');

结果如下：

DFT_diff_N_8_16_answer

观察到，输入有限序列xn2和xn3不同，但当做8点ffts时，二者结果相同，做16点fft时，二者结果不同
究其原因，首先画出xn2和xn3分别进行8和16点周期延拓的序列图像
从图中可以看到，进行8点周期延拓后，xn2和xn3是相同的周期信号，只是相位不同，但并不会影响DFT的幅频特性
进行16点周期延拓后，xn2和xn3显然不是同一个周期信号，故DFT的结果也会不一致
从上述分析得出结论，DFT是将M个有限序列进行N点（N>M时补零）进行周期延拓，然后选取一个周期内的序列，进行DFT计算

对同一序列，进行不同N的DFT计算，结果也会不同

clc;clear all;            %清屏清变量

fs=100;                   %采样频率
f=10;                     %信号频率
t=0:1/fs:10;              %时间序列

x=sin(2*pi*f*t);          %正弦信号

%绘制幅频响应
figure(1);
subplot(211);stem(abs(fft(x,40)));
subplot(212);stem(abs(fft(x,45)));

figure(2)
x40=[x(1:40),x(1:40),x(1:40),x(1:40)];
x45=[x(1:45),x(1:45),x(1:45),x(1:45)];
subplot(211);plot(x40); title('x40');
subplot(212);plot(x45); title('x45');

结果如下：

DFT_sin_diff_N_answer

从结果图中观察到对sin函数进行40点DFT和45点DFT是不一样
究其原因，首先可以画出在时域上分别采样40点和45点的时域波形
观察上图，因为在时域上取45点序列进行周期延拓后形成的波形，并不是一个正常的正弦波，故最后的DFT结果当然跟理想的sin频谱不一样，其实这里就是出现频谱泄露的原因，也就是后面为什么会有窗函数这个玩意
关于为什么40点是正常的正弦波而45点不是正常的正弦波的解释：
- 信号频率为10Hz，采样频率为100Hz，也就是一个周期采10个点，40个点也就是采了4个周期
- 而45并不是10的倍数，采样了4个多一点的周期，所以周期延拓后并不是一个正常的正弦波

借sin函数讨论不同N对幅频特性中幅度的影响

clc;clear all;            %清屏清变量

fs=100;                   %采样频率
f=10;                     %信号频率
t=0:1/fs:10;              %时间序列

x=sin(2*pi*f*t);          %正弦信号

%绘制幅频响应
figure(1);
subplot(411);stem(abs(fft(x,20)));title('DFT\_x20');
subplot(412);stem(abs(fft(x,80)));title('DFT\_x80');
subplot(413);stem(abs(fft(x,160)));title('DFT\_x160');
subplot(414);stem(abs(fft(x,1000)));title('DFT\_x1000');

结果如下：
这张图完美诠释了，N点double，FFT幅度会double
但前提是N<M，且N需要是$\frac{f_s}{f_0}$的倍数

2.循环卷积原理

公式如下：
$$
x_1(n)\circledast x_2(n)\rightarrow x_1(k)\cdot x_2(k)
$$

$$
x_1(n)\cdot x_2(n)\rightarrow \frac1N[x_1(k)\circledast x_2(k)]
$$
循环卷积长度：若$x(n):L,h(n):M$，那么循环长度$N$为：
$$
N=max(L,M)
$$
- 是对$L$或者$M$的长度补零，补到比$L+M-1$长，那么循环卷积和线性卷积相等
当循环卷积长度$\ge$线性卷积长度$(L+M-1)$时，二者计算得到的结果相等

循环卷积的计算示例：

			$a(n)$	1	2	3	4
$b(n)$翻转后的周期延拓	0	7	6	5	0	7	6	5	0	7	$c(0)=50$
翻转延拓后向右移动1	5	0	7	6	5	0	7	6	5	0	$c(1)=44$
翻转延拓后向右移动2	6	5	0	7	6	5	0	7	6	5	$c(2)=34$
翻转延拓后向右移动3	7	6	5	0	7	6	5	0	7	6	$c(3)=52$

MATLAB的验证

clc;
an=[1 2 3 4];
bn=[5 6 7];

%4点DFT验证循环卷积的计算
an_fft_4=fft(an,4);
bn_fft_4=fft(bn,4);
cn_circle_4=ifft((an_fft_4.*bn_fft_4),4)

%6点DFT验证循环卷积与线性卷积相等
an_fft_6=fft(an,6);
bn_fft_6=fft(bn,6);
cn_linear=conv(an,bn)
cn_circle_6=ifft((an_fft_6.*bn_fft_6),6)

%大于6点DFT时循环卷积的结果
an_fft_8=fft(an,8);
bn_fft_8=fft(bn,8);
cn_circle_8=ifft((an_fft_8.*bn_fft_8),8)

结果如下：

Reference

学以致用系列课程之数字信号处理]_哔哩哔哩_bilibili

ssy的小天地

数字信号处理之频域处理