一个频谱在区间 ( − ω 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) 确定。
取样频率不能太低,必须 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。
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)");