• 建議配合國寶老師的視頻食用，
  信號與線性系統分析 吳大正 郭寶龍

(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π1?F(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?πωc??sa(ωc?t)，之所以選取 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?πωc??n=?∞∑∞?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)");