取样定理的matlab实现,数字信号处理实验1：采样定理

clear all

clc

%初始信号

t = -0.1 : 0.0005 : 0.1;

N = 1000;

k = -N : N;

W = k * 2000 / N;

x = sin(2 * pi * 60 * t) + cos(2 * pi * 50 * t) + sinc(t/pi);%原信号

y = x * exp(-1i * t' * W) * 0.0005;% 傅里叶变换

y = abs(y);

subplot(4, 2, 1); plot(t, x); title('原信号时域');

subplot(4, 2, 2); plot(W, y); title('原信号频域');

%信号的采样

sampling = 1/60;%60hz

t = -0.1 : sampling : 0.1;

x_60Hz = sin(2 * pi * 60 * t) + cos(2 * pi * 50 * t) + sinc(t/pi); %采样后的信号

y_60Hz = x_60Hz * exp(-1i * t' * W) * sampling; % 采样后的傅里叶变换

y_60Hz = abs(y_60Hz);

subplot(4, 2, 3); stem(t, x_60Hz); title('60Hz采样信号时域');

subplot(4, 2, 4); plot(W, y_60Hz); title('60Hz采样信号频域');

sampling = 1/125;%125hz

t = -0.1 : sampling : 0.1;

x_125Hz = sin(2 * pi * 60 * t) + cos(2 * pi * 50 * t) + sinc(t/pi);

y_125Hz = x_125Hz * exp(-1i * t' * W) * sampling;

y_125Hz = abs(y_125Hz);

subplot(4, 2, 5); stem(t, x_125Hz); title('125Hz采样信号时域');

subplot(4, 2, 6); plot(W, y_125Hz); title('125Hz采样信号频域');

sampling = 1/200; %200hz

t = -0.1 : sampling : 0.1;

x_200Hz = sin(2 * pi * 60 * t) + cos(2 * pi * 50 * t) + sin(2 * pi * 30 * t);

y_200Hz = x_200Hz * exp(-1i * t' * W) * sampling;

y_200Hz = abs(y_200Hz);

subplot(4, 2, 7); stem(t, x_200Hz); title('200Hz采样信号时域');

subplot(4, 2, 8); plot(W, y_200Hz); title('200Hz采样信号频域');

% 信号恢复

figure;

n = -200 : 200;

sampling = 1/60;%欠采样

sam = n * sampling;

t = -0.1 : 0.0005 : 0.1;

f1 = sin(2 * pi * 60 * sam) + cos(2 * pi * 50 * sam) + sinc(sam/pi);

f_new1 = f1 * sinc((1/sampling) * (ones(length(sam), 1) * t - sam' * ones(1, length(t))));%复原函数

subplot(3, 2, 1); plot(t, f_new1); title('60Hz信号恢复');

sampling = 1/125;%临界采样

sam = n * sampling;

t = -0.1 : 0.0005 : 0.1;

f1 = sin(2 * pi * 60 * sam) + cos(2 * pi * 50 * sam) + sinc(sam/pi);

f_new2 = f1 * sinc((1/sampling) * (ones(length(sam), 1) * t - sam' * ones(1, length(t))));

subplot(3, 2, 3); plot(t, f_new2); title('125Hz信号恢复');

sampling = 1/200;%过采样

sam = n * sampling;

t = -0.1 : 0.0005 : 0.1;

f1 = sin(2 * pi * 60 * sam) + cos(2 * pi * 50 * sam) + sinc(sam/pi);

f_covery3 = f1 * sinc((1/sampling) * (ones(length(sam), 1) * t - sam' * ones(1, length(t))));

subplot(3, 2, 5); plot(t, f_covery3); title('200Hz信号恢复');

%计算误差

t = -0.1 : 0.0005 : 0.1;

miss = abs(x - f_new1);%取误差绝对值

subplot(3, 2, 2);plot(t,miss);axis([-0.1 0.1 0 3 ]);title('60hz误差');

t = -0.1 : 0.0005 : 0.1;

miss = abs(x - f_new2);

subplot(3, 2, 4);plot(t,miss);axis([-0.1 0.1 0 0.03 ]);title('125hz误差');

t = -0.1 : 0.0005 : 0.1;

miss = abs(x - f_covery3);

subplot(3, 2, 6);plot(t,miss);axis([-0.1 0.1 0 0.003 ]);title('200hz误差');