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信號恢復');