%用三角波函數進行頻域低通抽樣定理的證明
clc
clear
a=-5;%信號的起始時間
b=5;%信號的終止時間
tc=2;%sinc函數的最高頻分量值
fs=1/(32*tc);%為使信號能較平滑的畫出，設置的采樣時間間隔
ts=1/fs;%由上面的采樣時間間隔得到的采樣頻率
f=a:fs:b;
t=-ts/2+ts/length(f):ts/length(f):ts/2;
y=sinc(tc*f).^2;
y_spect=fftshift(abs(fft(y)));
%------
ts1=tc*2;%臨界采樣頻率
fs1=1/ts1;%臨界采樣的時間間隔
n1=fs1/fs;
f1=a:fs1:b;%臨界采樣時刻

x1=sinc(tc*f1).^2;%臨界采樣得到的采樣值
y1=zeros(1,n1*(length(x1)-1)+1);
y1(1:n1:n1*length(x1))=x1;
y1_spect=fftshift(abs(fft(y1)));

z1=convn(y1,sinc(f/fs1),'same');%重建信號
z1_spect=fftshift(abs(fft(z1)));%重建信號的頻譜
%-----
ts2=tc*1.5;%欠采樣頻率
fs2=1/ts2;%欠采樣的時間間隔
n2=fs2/fs;
f2=a:fs2:b;%欠采樣時刻

x2=sinc(tc*f2).^2;%欠采樣得到的采樣值
y2=zeros(1,n1*(length(x2)-1)+1);
y2(1:n2:n2*length(x2))=x2;
y2_spect=fftshift(abs(fft(y2)));

z2=convn(sinc(f/fs2),y2,'same');%重建信號
z2_spect=fftshift(abs(fft(z2)));%重建信號的頻譜
%----
ts3=tc*3;%過采樣頻率
fs3=1/ts3;%過采樣的時間間隔
n3=fs3/fs;
f3=a:fs3:b;%過采樣時刻

x3=sinc(tc*f3).^2;%過采樣得到的采樣值
y3=zeros(1,n3*(length(x3)-1)+1);
y3(1:n3:n3*length(x3))=x3;
y3_spect=fftshift(abs(fft(y3)));

z3=convn(y3,sinc(f/fs3),'same');%重建信號
z3_spect=fftshift(abs(fft(z3)));%重建信號的頻譜
%-------繪圖
figure(1)
subplot(3,2,1)
plot(f,y);title('原始信號的頻譜');grid
subplot(3,2,2)
plot(t,y_spect);title('原始信號');grid
axis([min(t),max(t),min(y_spect),max(y_spect)])
subplot(3,2,3)
stem(f1,x1);title('臨界采樣得到的抽樣點');grid
subplot(3,2,4)
plot(t,y1_spect);title('抽樣點的反傅里葉變換');grid
axis([min(t),max(t),min(y1_spect),max(y1_spect)])
subplot(3,2,5)
plot(f,z1);title('用臨界采樣得到的抽樣點重建的信號');grid
subplot(3,2,6)
plot(t,z1_spect);title('重建信號的反傅里葉變換');grid
axis([min(t),max(t),min(z1_spect),max(z1_spect)])

figure(2)
subplot(3,2,1)
plot(f,y);title('原始信號的頻譜');grid
subplot(3,2,2)
plot(t,y_spect);title('原始信號');grid
axis([min(t),max(t),min(y_spect),max(y_spect)])
subplot(3,2,3)
stem(f2,x2);title('欠采樣得到的抽樣點');grid
subplot(3,2,4)
plot(t,y2_spect);title('抽樣點的反傅里葉變換');grid
axis([min(t),max(t),min(y2_spect),max(y2_spect)])
subplot(3,2,5)
plot(f,z2);title('用欠采樣得到的抽樣點重建的信號');grid
subplot(3,2,6)
plot(t,z2_spect);title('重建信號的反傅里葉變換');grid
axis([min(t),max(t),min(z2_spect),max(z2_spect)])

figure(3)
subplot(3,2,1)
plot(f,y);title('原始信號的頻譜');grid
subplot(3,2,2)
plot(t,y_spect);title('原始信號');grid
axis([min(t),max(t),min(y_spect),max(y_spect)])
subplot(3,2,3)
stem(f3,x3);title('過采樣得到的抽樣點');grid
subplot(3,2,4)
plot(t,y3_spect);title('抽樣點的反傅里葉變換');grid
axis([min(t),max(t),min(y3_spect),max(y3_spect)])
subplot(3,2,5)
plot(f,z3);title('用過采樣得到的抽樣點重建的信號');grid
subplot(3,2,6)
plot(t,z3_spect);title('重建信號的反傅里葉變換');grid
axis([min(t),max(t),min(z3_spect),max(z3_spect)])