Matlab实现抽样定理

Matlab实现抽样定理

正弦信号的抽样：
首先时间跨度选择 -0.2 到 0.2，间隔0.0005取一个点，原信号取 sin⁡(2π*60t) ，则频率为60Hz。
由于需要输出原始信号的波形，我选择了手动编写代码进行傅里叶变换，有公式origin_F = origin * exp(-1i * t' * W) * 0.0005 傅里叶变换后的值，并取绝对值。
采样则调整取点的间隔就ok了。
恢复波形不太懂，在网上找到的方法：
f_covery = f_uncovery * sinc((1/Nsampling) * (ones(length(n_sam), 1) * t - n_sam' * ones(1, length(t))));
出处为：
http://www.mathsky.cn/wiki/index.php?search-fulltext-title-%D0%C5%BA%C5%BB%D6%B8%B4--all-0-within-time-desc-1
最后则可以输出波形和原始信号进行对比分析。

混合信号的抽样：
这里和正弦信号相比，只是待抽样信号不同了而已，但是混合信号我用的是正弦和余弦的叠加
sin(2 * pi * 60 * t) + cos(2 * pi * 25 * t) + sin(2 * pi * 30 * t)
由于抽样频率没变，依然是80Hz、121Hz、150Hz，所以得到的结果和上面的是不一样的
下面的结果图会有相应的分析

实现效果

正弦信号：

恢复的波形为

对比80Hz的信号和121Hz的信号可知，原信号为60Hz的信号，至少需要120Hz才能不失真地恢复信号，由图可得，80Hz的信号虽然还是正弦信号，但是相位信息已经失真了。121Hz和150Hz的抽样信息则准确地恢复了原信号

混合信号：

因为只有原信号和下面的代码不一样，所以节省国家树木资源便不全部截图代码了。
不一样的地方为原信号为混合信号：

%% 设置原始信号
t = -0.2 : 0.0005 : 0.2;
N = 1000;
k = -N : N;
W = k * 2000 / N;
origin = sin(2 * pi * 60 * t);% 原始信号为正弦信号
origin_F = origin * exp(-1i * t' * W) * 0.0005;% 傅里叶变换
origin_F = abs(origin_F);% 取正值
figure;
subplot(4, 2, 1); plot(t, origin); title('原信号时域');
subplot(4, 2, 2); plot(W, origin_F); title('原信号频域');

运行效果图：

这个信号明显地可以看出80Hz采样的失真情况。由于混合信号中频率最高的那个信号为60Hz，因此也是至少需要120Hz才能不失真地恢复原始信号。

代码实现

clear all
clc
%% 设置原始信号
t = -0.2 : 0.0005 : 0.2;
N = 1000;
k = -N : N;
W = k * 2000 / N;
origin = sin(2 * pi * 60 * t) + cos(2 * pi * 25 * t) + sin(2 * pi * 30 * t);% 原始信号为正弦信号叠加
origin_F = origin * exp(-1i * t' * W) * 0.0005;% 傅里叶变换
origin_F = abs(origin_F);% 取正值
figure;
subplot(4, 2, 1); plot(t, origin); title('原信号时域');
subplot(4, 2, 2); plot(W, origin_F); title('原信号频域');
%% 对原始信号进行80Hz采样率采样
Nsampling = 1/80; % 采样频率
t = -0.2 : Nsampling : 0.2;
f_80Hz = sin(2 * pi * 60 * t) + cos(2 * pi * 25 * t) + sin(2 * pi * 30 * t); %采样后的信号
F_80Hz = f_80Hz * exp(-1i * t' * W) * Nsampling; % 采样后的傅里叶变换
F_80Hz = abs(F_80Hz);
subplot(4, 2, 3); stem(t, f_80Hz); title('80Hz采样信号时域');
subplot(4, 2, 4); plot(W, F_80Hz); title('80Hz采样信号频域');
%% 对原始信号进行121Hz采样率采样
Nsampling = 1/121; % 采样频率
t = -0.2 : Nsampling : 0.2;
f_80Hz = sin(2 * pi * 60 * t) + cos(2 * pi * 25 * t) + sin(2 * pi * 30 * t); %采样后的信号
F_80Hz = f_80Hz * exp(-1i * t' * W) * Nsampling; % 采样后的傅里叶变换
F_80Hz = abs(F_80Hz);
subplot(4, 2, 5); stem(t, f_80Hz); title('121Hz采样信号时域');
subplot(4, 2, 6); plot(W, F_80Hz); title('121Hz采样信号频域');
%% 对原始信号进行150Hz采样率采样
Nsampling = 1/150; % 采样频率
t = -0.2 : Nsampling : 0.2;
f_80Hz = sin(2 * pi * 60 * t) + cos(2 * pi * 25 * t) + sin(2 * pi * 30 * t); %采样后的信号
F_80Hz = f_80Hz * exp(-1i * t' * W) * Nsampling; % 采样后的傅里叶变换
F_80Hz = abs(F_80Hz);
subplot(4, 2, 7); stem(t, f_80Hz); title('150Hz采样信号时域');
subplot(4, 2, 8); plot(W, F_80Hz); title('150Hz采样信号频域');
%% 恢复原始信号
% 从80Hz采样信号恢复
figure;
n = -100 : 100;
Nsampling = 1/80;
n_sam = n * Nsampling;
f_uncovery = sin(2 * pi * 60 * n_sam) + cos(2 * pi * 25 * n_sam) + sin(2 * pi * 30 * n_sam);
t = -0.2 : 0.0005 : 0.2;
f_covery = f_uncovery * sinc((1/Nsampling) * (ones(length(n_sam), 1) * t - n_sam' * ones(1, length(t))));
subplot(3, 1, 1); plot(t, f_covery); title('80Hz信号恢复');
% 从121Hz采样信号恢复
Nsampling = 1/121;
n_sam = n * Nsampling;
f_uncovery = sin(2 * pi * 60 * n_sam) + cos(2 * pi * 25 * n_sam) + sin(2 * pi * 30 * n_sam);
t = -0.2 : 0.0005 : 0.2;
f_covery = f_uncovery * sinc((1/Nsampling) * (ones(length(n_sam), 1) * t - n_sam' * ones(1, length(t))));
subplot(3, 1, 2); plot(t, f_covery); title('121Hz信号恢复');
% 从150Hz采样信号恢复
Nsampling = 1/150;
n_sam = n * Nsampling;
f_uncovery = sin(2 * pi * 60 * n_sam) + cos(2 * pi * 25 * n_sam) + sin(2 * pi * 30 * n_sam);
t = -0.2 : 0.0005 : 0.2;
f_covery = f_uncovery * sinc((1/Nsampling) * (ones(length(n_sam), 1) * t - n_sam' * ones(1, length(t))));
subplot(3, 1, 3); plot(t, f_covery); title('150Hz信号恢复');