数字信号处理 --- 信号的采样和奇妙的混叠（Aliasing）壹

信号的采样与混叠(时域)

在一定条件下，一个连续的时间信号完全可以用该信号在等时间间隔上的样本来表示，并且可以用这些样本的值把该信号完全恢复出来。这一非常重要的发现被称之为采样定理。

采样定理在连续时间信号（一个是连续的波形）和离散时间信号（一个是数组）之间架起了一座桥梁！！！（就好像微积分基本定理一样，把斜率和面积这两个毫不相干的概念紧密的联系在了一起。）

一维信号的采样：

通过一个固定周期的冲击串来实现信号的采样（即，连续信号的离散化）。如下图，一个本来是连续的一维信号被分割成了很多小块，这非常像黎曼积分（Riemann integral），定积分作黎曼和的极限！

上图中原来是连续的数据通过等间隔（采样周期为T）的采样，被分割成了等间隔，具有不同幅度的冲击串，见下图。

附：下图为黎曼和的极限，图中矩形的宽度就好比是采样间隔T！

二维信号的采样：

二维信号的采样常被用于图像的数字化，相对于一维采样，他只是增加了采样的维度。不仅在X方向采样还在Y方向采样，采样后的结果是一个二维矩阵，矩阵中的每一个点代表了该位置的响应灰度。

例如，在出版业里，一幅图是由很多很小的采样点构成的，如果这些采样点的间距很小，小到你无法用肉眼分辨的话。（或者说，即便报纸上的采样点间距不是很小，但是读者眼睛距离报纸不是很近的话。）那么，报纸上的照片在我们看来就是一幅连续的图片。但是，如果你拿着放大镜或者眼睛靠着很近观察的话，你就会发现报纸上的一幅图像其实是由一个个等间隔的离散化的点拼出来的。

下面我们用Matlab软件来仿真信号的采样，以及如果采样频率（采样周期的倒数）不能满足一定要求是所带来的混叠。

Tips : 本文的后面会有一个录像说明，现在的MATLAB仿真和录像说明有一个很大的不同就是。这里的仿真是保持原始信号的频率Fo不变，不断地改变采样频率Fs.

第一种情况 ------ 采样频率高于原始信号频率最高频的两倍（Fs > 2*Fo）

下图所绘的是一个频率为60Hz的连续的余弦波，图中的红点为采样点。当采样冲击串的频率（2000Hz）远高于原始信号时（原始信号的33.3倍），采样点在曲线上的分布非常密集，重建后的效果也很好。如果采样频率（400Hz）只是略高于原始信号频率

（原始信号的6.6倍）,则原曲线上的采样点明显不足，但不用担心，即便是这种看起来不太充足的采样还是可以重建原始的连续信号的。

点击图片放大

Matlab代码：

%% signal sampling and aliasing demo for CSDN
% Created: Late spring, 2018. (2018/05/03)
% Author: Z.Zhu, zzz0029@tigermail.auburn.edu
% Copy Rights Reserved.
% May not be copied, scanned, or duplicated, in whole or in part.
% Reference of algorithm by: Steve Eddins's Aliasing and a sampled cosine
% signal on his bolg "Steve on Image Processing"

close all;
clear all;
%% Signals with same frequency were sampled with different sampling frequencies.  
Omega = 60;                     % freq. of continous Signal 
Phi = 0;                        % Phase
A = 1;                          % Amplitude
tmin = -pi/55;                  % lower limit of X axis
tmax = pi/55;                   % uper limit of X axis
Bins = 1000;                    % Number of Bins
t = linspace(tmin, tmax, Bins); % Space vector

% Form a signal containing a 60(Hz) sinusoid of amplitude 1 with phase of 0. 
CosineSignal = A.*cos(2.*pi.*Omega.*t + Phi);

% Plot the signal in the Space domain.
subplot(2,2,1)
plot(t,CosineSignal,'-k','LineWidth',2);
title(['A continuous-time cosine signal at ',num2str(Omega),'Hz',',sampled with 2000Hz sampling frequency.'])
ylabel('Amplitude')
xlabel('t (seconds)')

% Let's sample it with a sampling frequency of 2000Hz.
Fs = 2000;   % Sampling Freq.
T = 1/Fs;   % Space vector
nmin = ceil(tmin / T);
nmax = floor(tmax / T);
n = nmin:nmax;
SampleSignal2000 = A.*cos(2.*pi.*Omega.*n.*T + Phi);
hold on
plot(n*T,SampleSignal2000,'o',...
    'MarkerSize',6.5,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor','r');
legend('Continuous-time signal','Samples','Location','northeast')
hold off
grid on

% Rebuild discrete sampled signals
subplot(2,2,3)
stem(n*T,SampleSignal2000,'filled','LineWidth',2);
title('Rebuild discrete sampled signals');

% The sampling frequency of 2000 Hz is well above 120 Hz, which is twice the frequency 
% of the cosine. And you can see that the samples are clearly capturing the oscillation 
% of the continuous-time cosine.

% Plot the signal in the Space domain.
subplot(2,2,2)
plot(t,CosineSignal,'-k','LineWidth',2);
title(['A continuous-time cosine signal at ',num2str(Omega),'Hz',',sampled with 400Hz sampling frequency.'])
ylabel('Amplitude')
xlabel('t (seconds)')

% Let's sample it with a sampling frequency of 400Hz.
Fs = 400;   % Sampling Freq.
T = 1/Fs;   % Space vector
nmin = ceil(tmin / T);
nmax = floor(tmax / T);
n = nmin:nmax;
SampleSignal400 = A.*cos(2.*pi.*Omega.*n.*T + Phi);
hold on
plot(n*T,SampleSignal400,'o',...
    'MarkerSize',6.5,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor','r');
legend('Continuous-time signal','Samples','Location','northeast')
hold off
grid on

% Rebuild discrete sampled signals
subplot(2,2,4)
stem(n*T,SampleSignal400,'filled','LineWidth',2);
title('Rebuild discrete sampled signals');

第二种情况 ------ 采样频率小于原始信号频率最高频的两倍，但大于原始信号的频率（2*Fo > Fs > Fo）

下面我们逐渐降低采样率，看看会发生什么。

先是用正好两倍于原始信号频率的采样频率（120Hz）来采样（原始信号的2倍），则采到的信号，正好是余弦波的波峰或者波谷即，1和-1。（注意：我这里所使用的是相位为0的余弦波，如果相位发生改变则采集到的点，不一定正好是波峰和波谷。）

接下来，我用小于原始信号频率两倍的采样频率（70Hz）来采时（原始信号的1.16倍），采集到的信号就会发生混叠，即采集到的是一个频率低于原始信号的错误的离散信号，也无法利用该错误信号重建原始的连续信号。

点击图片放大

Matlab代码：

% Let's sample it with a sampling frequency of 	120Hz.
Fs = 120;   % Sampling Freq.
T = 1/Fs;   % Space vector
nmin = ceil(tmin / T);
nmax = floor(tmax / T);
n = nmin:nmax;
SampleSignal120 = A.*cos(2.*pi.*Omega.*n.*T + Phi);
hold on
plot(n*T,SampleSignal120,'o',...
    'MarkerSize',6.5,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor','r');
legend('Continuous-time signal','Samples','Location','northeast')
hold off
grid on

% Rebuild discrete sampled signals
subplot(2,2,3)
stem(n*T,SampleSignal120,'filled','LineWidth',2);
title('Rebuild discrete sampled signals');

% The samples above are still adequately capturing the shape of the cosine.Now let's drop 
% the sampling frequency down to exactly 120 Hz, twice the frequency of the 60 Hz cosine. 

% Plot the signal in the Space domain.
subplot(2,2,2)
plot(t,CosineSignal,'-k','LineWidth',2);
title(['A continuous-time cosine signal at ',num2str(Omega),'Hz',',sampled with 70Hz sampling frequency.'])
ylabel('Amplitude')
xlabel('t (seconds)')

% Let's sample it with a sampling frequency of 70Hz.
Fs = 70;   % Sampling Freq.
T = 1/Fs;   % Space vector
nmin = ceil(tmin / T);
nmax = floor(tmax / T);
n = nmin:nmax;
SampleSignal70 = A.*cos(2.*pi.*Omega.*n.*T + Phi);
hold on
plot(n*T,SampleSignal70,'o',...
    'MarkerSize',6.5,...
    'MarkerEdgeColor','r',...
    'MarkerFaceColor','r');
legend('Continuous-time signal','Samples','Location','northeast')
hold off
grid on

% Rebuild discrete sampled signals
subplot(2,2,4)
stem(n*T,SampleSignal70,'filled','LineWidth',2);
title('Rebuild discrete sampled signals');

下图是利用上图中的采样点重建后的离散信号，图中我用红色的箭头大致的链接了一下图像中的离散点（真实重建的时候，不会使用直线来连接他们的，而是使用的Sinc函数。），下图的上半部分图像中的起伏和波动的频率和原始信号是相同的，而下半部分的图像的波动就明显慢了很多。