(1) 采样的说明
1. 为什么要取样
- 我们要使用计算机去处理信号。信号都是连续的，在任意一个非空区间内都有无限个值，但是计算机的内存是一个有限值，只能完成有限的数据存储和运算，所以要进行采样和量化等操作，把信号离散成有限个点就可以使用计算机处理了。当信号在计算机中处理完，按照一定的规则恢复成连续的状态就可以了。
- 注意采样定理和后面 DFT（discrete Fourier transform） 之间的联系。
2. 什么是取样
- 取样就是利用取样脉冲序列从连续信号中 抽取离散的样本值。
- 从时域上看
f
s
(
t
)
=
f
(
t
)
×
s
(
t
)
f_s(t) = f(t) \times s(t)
fs(t)=f(t)×s(t)
- 从频域上看
F
s
(
j
ω
)
=
1
2
π
F
(
j
ω
)
∗
S
(
j
ω
)
F_s(j\omega) = \displaystyle\frac{1}{2\pi}F(j\omega) * S(j\omega)
Fs(jω)=2π1F(jω)∗S(jω)
- 冲激取样（理想取样）的例子
- 可以看到在时域上是通过乘积完成对应离散点的选取，在频域上完成的是信号频谱的周期延拓。
(2) 采样定理
1. 为什么要有奈奎斯特频率
- 首先要明白一件事，使用采样点恢复出原来的信号需要的是进行一个低通滤波，把频域上的低频波形滤出来就可以了。
- 从上面的图可以看出来，取样在频域中相当于进行了频谱的周期延拓，所以就会出现一个问题，平移距离不够会引起频域波形的重叠，当发生重叠之后就无法完成滤波了，即无法完成信号的恢复了。因此需要对平移的距离有一个限制，也就是采样定理中的奈奎斯特频率。
2. 什么是采样定理
- 一个频谱在区间
(
−
ω
m
,
ω
m
)
(-\omega_m,\omega_m)
(−ωm,ωm) 以外为 0 的带限信号
f
(
t
)
f(t)
f(t)，可唯一的由其在均匀间隔
T
s
[
T
s
<
2
π
/
ω
m
]
T_s [T_s < 2\pi/\omega_m]
Ts[Ts<2π/ωm] 上的样值点
f
(
n
T
s
)
f(nT_s)
f(nTs) 确定。
- 注意必须是带限信号。像冲激函数这样的就无法取样，因为在频域上的无限意味着在时域信号的存在时间无限趋于 0。
- 取样频率不能太低，必须
f
s
>
2
f
m
f_s > 2f_m
fs>2fm 。最低取样频率
f
s
=
2
f
m
f_s = 2f_m
fs=2fm称为奈奎斯特频率。
(3) 信号的恢复
- 参量的说明
- 低通滤波器的截止角频率：
ω
c
\omega_c
ωc，从图上明显可以看出需要有
ω
m
<
ω
c
<
ω
s
−
ω
m
\omega_m<\omega_c<\omega_s-\omega_m
ωm<ωc<ωs−ωm，为方便取
ω
c
=
0.5
ω
s
\omega_c = 0.5\omega_s
ωc=0.5ωs。
- 采样角频率：
ω
s
\omega_s
ωs，注意根据采样定理
ω
s
>
2
ω
m
\omega_s > 2\omega_m
ωs>2ωm。
- 带限信号的最大角频率：
ω
m
\omega_m
ωm。
- 信号恢复的过程是信号采样过程的逆过程。实际上在进行频域乘积滤波的过程中，时域进行了卷积的平移，平移到不同位置的函数叠加就恢复出原始信号。
(4) Matlab的Sa函数取样仿真
1. 采样信号Sa函数的说明
- Matlab 中自带的函数是
sinc 函数，其形式是
s
i
n
(
π
t
)
π
t
\displaystyle\frac{sin(\pi t)}{\pi t}
πtsin(πt)，我们要在仿真中使用的是 Sa 函数，其形式是
s
i
n
(
t
)
t
\displaystyle\frac{sin(t)}{t}
tsin(t)，因此 sa = sinc(t/pi)。
- 代码：
%% 打印出来sa函数
t = -20:0.001:20;
L = length(t);
x = sinc(t / pi);
plot(t,x,'LineWidth',3);
xlabel('t');ylabel('Amplitude'); title('Sa(t)')
- 结果：
2. 进行参数的说明及相关计算
- 参数说明
-
s
a
(
t
)
sa(t)
sa(t) 的傅里叶变换结果是
π
g
2
(
ω
)
\pi g_2(\omega)
πg2(ω)，就是一个门宽为 2 的门函数。因此可以知道
ω
m
=
1
\omega_m = 1
ωm=1。
- 根据奈奎斯特采样定律，这里选取
ω
s
=
2
ω
m
\omega_s = 2\omega_m
ωs=2ωm，
ω
s
=
1.5
ω
m
\omega_s=1.5\omega_m
ωs=1.5ωm，
ω
s
=
4
ω
m
\omega_s=4\omega_m
ωs=4ωm。分别模拟临界采样，欠采样和过采样三种情况。相应的选取信号还原时低通滤波器的截止频率
ω
c
=
0.5
ω
s
\omega_c = 0.5\omega_s
ωc=0.5ωs。
- 这里选取时域的正半轴取样点一共 N 个，下面使用
∞
\infin
∞ 推公式，但是最后要用
N
N
N。
- 信号取样
- 冲激取样函数：
δ
T
s
(
t
)
=
∑
n
=
−
∞
∞
δ
(
t
−
n
T
s
)
\delta_{T_s}(t)=\displaystyle\sum_{n=-\infin}^{\infin}\delta(t-nT_s)
δTs(t)=n=−∞∑∞δ(t−nTs)。
- 通过采样的定义可知
f
s
(
t
)
=
f
(
t
)
×
s
a
(
t
)
f_s(t) = f(t) \times sa(t)
fs(t)=f(t)×sa(t)，在matlab中只需要
fs = sinc(t/pi)。
- 信号恢复
- 采样后的信号在时域上的表达式为
f
s
(
t
)
=
f
(
t
)
∑
n
=
−
∞
∞
δ
(
t
−
n
T
s
)
=
∑
n
=
−
∞
∞
δ
(
t
−
n
T
s
)
f
(
n
T
s
)
f_s(t)=f(t)\displaystyle\sum_{n=-\infin}^{\infin}\delta(t-nT_s)=\displaystyle\sum_{n=-\infin}^{\infin}\delta(t-nT_s)f(nT_s)
fs(t)=f(t)n=−∞∑∞δ(t−nTs)=n=−∞∑∞δ(t−nTs)f(nTs)
- 假设采样后的信号在频域上的表达式为
F
s
(
j
ω
)
F_s(j\omega)
Fs(jω)，并选取低通滤波器
H
(
ω
)
=
{
T
s
,
∣
ω
∣
≤
ω
c
0
,
∣
ω
∣
>
ω
c
H(\omega)=\begin{cases} T_s ,&|\omega|\leq \omega_c\\ 0, & |\omega|> \omega_c \end{cases}
H(ω)={Ts,0,∣ω∣≤ωc∣ω∣>ωc
可以算出
H
(
ω
)
H(\omega)
H(ω) 在时域上的表达式
h
(
t
)
=
T
s
ω
c
π
s
a
(
ω
c
t
)
h(t)=T_s\displaystyle\frac{\omega_c}{\pi}sa(\omega_ct)
h(t)=Tsπωcsa(ωct)。之所以选取
H
(
ω
)
H(\omega)
H(ω)的放大倍数为
T
s
T_s
Ts 是因为此时
h
(
t
)
h(t)
h(t) 的系数是 1（因为
ω
c
=
0.5
ω
s
\omega_c = 0.5\omega_s
ωc=0.5ωs）。
- 根据前面的讨论，让取样后的信号通过低通滤波器相当于频域相乘即
F
(
j
ω
)
=
F
s
(
j
ω
)
×
H
(
ω
)
F(j\omega) = F_s(j\omega)\times H(\omega)
F(jω)=Fs(jω)×H(ω)。同时根据时域和频域的关系，
f
(
t
)
=
f
s
(
t
)
∗
h
(
t
)
f(t) = f_s(t) * h(t)
f(t)=fs(t)∗h(t)。带入前面的结果可以得到
f
(
t
)
=
T
s
ω
c
π
∑
n
=
−
∞
∞
f
(
n
T
s
)
s
a
(
ω
c
(
t
−
n
T
s
)
)
f(t)=T_s\displaystyle\frac{\omega_c}{\pi}\displaystyle\sum_{n=-\infin}^{\infin}f(nT_s)sa(\omega_c(t-nT_s))
f(t)=Tsπωcn=−∞∑∞f(nTs)sa(ωc(t−nTs))
3. 结果的展示
- 过采样（实际上这里有一点不太明白，为什么过采样恢复后信号的误差会比临界采样的大？？）
4. matlab 代码
%% matlab 完成Sa信号的采样和恢复
%% 取样（临界取样）
% 取样
figure(1);
wm = 1; %信号的最大频率
ws = 2 * wm; %信号的采样频率（根据奈奎斯特频率）
wc = 0.5 * ws;%滤波器的截止频率
Ts = 2*pi/ws;%采样间隔
N = 10;%时域采样点数
n = -N:N;
nTs = n * Ts;%采样数据的采样时间
fs = sinc(nTs/pi);%完成采样
subplot(311);
stem(nTs/pi,fs,'LineWidth',3);
xlabel("nTs");
ylabel("f(nTs)");
title("sa(t)的临界取样信号");
% 还原
Dt = 0.005;
t = -15:Dt:15;
fa = Ts*wc/pi * fs * sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));
subplot(312);
plot(t,fa,'LineWidth',3);
xlabel("t");
ylabel("f(t)");
title("由临界取样信号重构sa(t)");
% 展示误差
error = abs(fa-sinc(t/pi));
subplot(313);
plot(t,error,'LineWidth',3);
xlabel("t");
ylabel("error(t)");
title("重构信号与原信号的误差error(t)");
%% 取样（过取样）
% 取样
figure(2);
wm = 1; %信号的最大频率
ws = 4 * wm; %信号的采样频率（根据奈奎斯特频率）
wc = 0.5 * ws;%滤波器的截止频率
Ts = 2*pi/ws;%采样间隔
N = 20;%时域采样点数
n = -N:N;
nTs = n * Ts;%采样数据的采样时间
fs = sinc(nTs/pi);%完成采样
subplot(311);
stem(nTs/pi,fs,'LineWidth',3);
xlabel("nTs");
ylabel("f(nTs)");
title("sa(t)的过取样信号");
% 还原
Dt = 0.005;
t = -15:Dt:15;
fa = fs*Ts*wc/pi * sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));
subplot(312);
plot(t,fa,'LineWidth',3);
xlabel("t");
ylabel("f(t)");
title("由过取样信号重构sa(t)");
% 展示误差
error = abs(fa-sinc(t/pi));
subplot(313);
plot(t,error,'LineWidth',3);
xlabel("t");
ylabel("error(t)");
title("重构信号与原信号的误差error(t)");
%% 取样（欠取样）
% 取样
figure(3);
wm = 1; %信号的最大频率
ws = 1.5 * wm; %信号的采样频率（根据奈奎斯特频率）
wc = 0.5 * ws;%滤波器的截止频率
Ts = 2*pi/ws;%采样间隔
N = 7;%时域采样点数
n = -N:N;
nTs = n * Ts;%采样数据的采样时间
fs = sinc(nTs/pi);%完成采样
subplot(311);
stem(nTs/pi,fs,'LineWidth',3);
xlabel("nTs");
ylabel("f(nTs)");
title("sa(t)的欠取样信号");
% 还原
Dt = 0.005;
t = -15:Dt:15;
fa = fs*Ts*wc/pi * sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));
subplot(312);
plot(t,fa,'LineWidth',3);
xlabel("t");
ylabel("f(t)");
title("由欠取样信号重构sa(t)");
% 展示误差
error = abs(fa-sinc(t/pi));
subplot(313);
plot(t,error,'LineWidth',3);
xlabel("t");
ylabel("error(t)");
title("重构信号与原信号的误差error(t)");
现在已经转移到知乎，之后的文章会在知乎更新。
