PCM抽样的MATLAB程序设计按如下步骤进行：
（1）确定输入的模拟信号为sa(200t)；
（2）根据输入的模拟信号，确定抽样频率，对输入信号进行抽样，并将正常抽样和会产生失真的抽样进行对比，对抽样定理加以验证；
（3）编写程序，画出满足采样定理和不满足的时、频域图形。
MATLAB源程序：
function sample()
t0=10; %定义时间长度
ts=0.001; fs=1/ts;
t=[-t0/2:ts:t0/2]; %定义时间序列
df=0.5; %定义频率分辨率
x=sin(200*t); m=x./(200*t+eps);
w=t0/(2*ts)+1; %确定t=0的点
m(w)=1; %修正t=0点的信号值
m=m.*m;
[M,mn,dfy]=fft_seq(m,ts,df); %傅立叶变换
M=M/fs;
f=[0:dfy:dfy*length(mn)-dfy]-fs/2; %定义频率序列
figure(1)
subplot(2,1,1); plot(t,m);
xlabel('时间');ylabel('幅值');title('原始信号(fh=200/2piHz)的波形');
axis([-0.15,0.15,0,1.5]);
subplot(2,1,2);
plot(f,abs(fftshift(M)));
xlabel('频率');ylabel('幅值');
axis([-500,500,0,0.03]);title('原始信号的频谱');
t0=10; %信号持续的时间
ts1=0.005; %满足抽样条件的抽样间隔
fs1=1/ts1;
t1=[-t0/2:ts1:t0/2]; %定义满足抽样条件的时间序列
x1=sin(200*t1); m1=x1./(200*t1+eps); w1=t0/(2*ts1)+1;
m1(w1)=1; %修正t=0时的信号值
m1=m1.*m1; %定义信号
[M1,mn1,df1]=fft_seq(m1,ts1,df); %对满抽样条件的信号进行傅立叶变换
M1=M1/fs1;N1=[M1,M1,M1,M1,M1,M1,M1,M1,M1,M1,M1,M1,M1];
f1=[-7*df1*length(mn1):df1:6*df1*length(mn1)-df1]-fs1/2;
figure(2)
subplot(2,1,1); stem(t1,m1);
xlabel('时间');ylabel('幅值');
title('抽样正常(fs=200Hz)时的信号波形');axis([-0.15,0.15,0,1]);
subplot(2,1,2)
plot(f1,abs(fftshift(N1)));
xlabel('频率');ylabel('幅值');axis([-500,500,0,0.05]);
title('抽样正常时的信号频谱');axis([-500,500,-0.01,0.03]);
t0=10; %信号持续的时间
ts2=0.01; %不满足抽样条件的抽样间隔
fs2=1/ts2;
t2=[-t0/2:ts2:t0/2]; %定义不满足抽样条件的时间序列
x2=sin(200*t2); m2=x2./(200*t2+eps); w2=t0/(2*ts2)+1;
m2(w2)=1; %修正t=0时的信号值
m2=m2.*m2; %定义信号
[M2,mn2,df2]=fft_seq(m2,ts2,df);%对不满足抽样条件的信号进行傅立叶变换
M2=M2/fs2;N2=[M2,M2,M2,M2,M2,M2,M2,M2,M2,M2,M2,M2,M2];
f2=[-7*df2*length(mn2):df2:6*df2*length(mn2)-df2]-fs2/2;
figure(3)
subplot(2,1,1); stem(t2,m2);
xlabel('时间');ylabel('幅值');title('抽样失真(fs=100Hz)时的信号波形');
axis([-0.15,0.15,0,1]);subplot(2,1,2)
plot(f2,abs(fftshift(N2)));
xlabel('频率');ylabel('幅值');axis([-500,500,0,0.02]);
title('抽样失真时的信号频谱');axis([-500,500,0.005,0.02]);
function [M,m,df]=fft_seq(m,ts,df)
fs=1/ts;
if nargin==2 n1=0
else n1=fs/df
end
n2=length(m);n=2^(max(nextpow2(n1),nextpow2(n2)));
M=fft(m,n);m=[m,zeros(1,n-n2)];df=fs/n
PCM抽样仿真结果：
原始信号的波形和频谱
PCM正常抽样时信号的波形及频谱
PCM抽样失真时信号的波形及频谱
3.3.2 PCM均匀量化的仿真
PCM均匀量化的MATLAB程序设计按如下步骤进行：
（1）确定输入模拟信号为sin(t)；
（2）根据均匀量化的原理均匀量化的算法程序；
（3）绘制并比较模拟输入信号与量化输出的波形。
PCM均匀量化源程序如下：
function average()
t=[0:0.01:4*pi];
y=sin(t);
w=jylh(y,1,64);
subplot(2,1,1);
plot(t,y);
xlabel('时间');
ylabel('幅度');
axis([0,4*pi,-1.1,1.1]);
title('原始信号');
subplot(2,1,2);
plot(t,w);
xlabel('时间');
ylabel('幅度');
axis([0,4*pi,-1.1,1.1]);
title('均匀量化后的信号');
function h=jylh(f,V,L)
n=length(f);t=2*V/L;
p=zeros(1,L+1);
for i=1:L+1,p(i)=-V+(i-1)*t;end
for i=1:n
if f(i)>V,h(i)=V;end
if f(i)<=-V,h(i)=-V;end
flag=0;
for j=2:L/2+1
if(flag==0)
if(f(i)<p(j))
h(i)=p(j-1);
flag=1;
end;
end;
end;
for j=L/2+2:L+1
if(flag==0)
if(f(i)<p(j))
h(i)=p(j);
flag=1;
end
end
end
end
nq=V^2/(3*L^2);
仿真结果：
3.3.3 PCM A律非均匀量化的MATLAB实现
PCM A律非均匀量化的MATLAB程序设计按如下步骤进行：
（1）确定输入模拟信号；
（2）根据非均匀量化的原理确定A律非均匀量化的算法程序；
（3）绘制并比较模拟输入信号与量化输出的波形。
源程序如下：
function a_quantize()
t=0:0.00000125:0.0005;
y=sin(8000*pi*t);
figure
subplot(2,1,1)
plot(t,y)
axis([0 0.0005 -1.2 1.2])
xlabel('时间') ylabel('幅度')
title('原始信号')
z=a_pcm(y,87.6);
subplot(2,1,2) plot(t,z)
axis([0 0.0005 -1.2 1.2])
xlabel('时间')
ylabel('幅度')
title('A律量化后的信号')
function y=a_pcm(x,a)
t=1/a;
for i=1:length(x)
if x(i)>=0
if(x(i)<=t)
y(i)=(a*x(i))/(1+log(a));
else
y(i)=(1+log(a*x(i)))/(1+log(a));
end
else
if(x(i)>=-t)
y(i)=-(a*-x(i))/(1+log(a));
else
y(i)=-(1+log(a*-x(i)))/(1+log(a));
end
end
end
仿真结果：
PCM均匀量化的MATLAB程序设计按如下步骤进行：
（1）确定输入模拟信号；
（2）根据给均匀量化的原理确定非均匀量化的算法程序；
（3）将上述编码的十进制数转化成8位二进制数。
源程序如下：
function a_13code()
t=0:0.000025:0.00025;
y=sin(8000*pi*t)
z=line13(y)
c=pcmcode(z)
function y=line13(x)
x=x/max(x);
z=sign(x);
x=abs(x);
for i=1:length(x)
if((x(i)>=0)&(x(i)<1/64))
y(i)=16*x(i);
else
if((x(i)>=1/64)&(x(i)<1/32))
y(i)=8*x(i)+1/8;
else
if((x(i)>=1/32)&(x(i)<1/16))
y(i)=4*x(i)+2/8;
else
if((x(i)>=1/16)&(x(i)<1/8))
y(i)=2*x(i)+3/8;
else
if((x(i)>=1/8)&(x(i)<1/4))
y(i)=x(i)+4/8;
else
if((x(i)>=1/4)&(x(i)<1/2))
y(i)=1/2*x(i)+5/8;
else
if((x(i)>=1/2)&(x(i)<=1))
y(i)=1/4*x(i)+6/8;
end
end
end
end
end
end
end
end
y=z.*y;
function f=pcmcode(y)
f=zeros(length(y),8);
z=sign(y);
y=y.*128;
y=fix(y);
y=abs(y);
for i=1:length(y)
if(y(i)==128)
y(i)=127.999;
end
end
for i=1:length(y)
for j=6:-1:0
f(i,8-j)=fix(y(i)/2^j);
y(i)=mod(y(i),(2^j));
end
end
for i=1:length(y);
if(z(i)==1)
f(i,1)=0;
else
f(i,1)=1;
end
end
程序运行结果：
y =
0 0.5878 0.9511 0.9511 0.5878 0.0000 -0.5878 -0.9511 -0.9511 -0.5878 -0.0000
z =
0 0.9045 1.0000 1.0000 0.9045 0.0000 -0.9045 -1.0000 -1.0000 -0.9045 -0.0000
c =
1 0 0 0 0 0 0 0
0 1 1 1 0 0 1 1
0 1 1 1 1 1 1 1
0 1 1 1 1 1 1 1
0 1 1 1 0 0 1 1
0 0 0 0 0 0 0 0
1 1 1 1 0 0 1 1
1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1
1 1 1 1 0 0 1 1
1 0 0 0 0 0 0 0
4.信道编码和译码
在数字电视和通信系统中，为提高信息传输可靠性，广泛使用了具有一定纠错能力的信道编码技术，如奇偶校验码、行列监督码、恒比码、汉明码、循环码（CRC）等编码技术。信道编码的本质是增加通信的可靠性，或者说增加整个系统的抗干扰性。对信道编码有以下要求：1.透明性：要求对所传消息的内容不加任何限制；2.有纠错能力；3.效率高：为了与信道频谱匹配和具有纠错能力，通常要向原信号添加一些码，要求加入最少的比特数而得到最大的利益；4.包含适当的定时信息。在这些要求中，除编码的必须信息外，所作的处理主要有两条：一是要求码列的频谱特性适应通道的频谱特性从而使传输过程中能量损失最小，提高信噪比。减少发生差错的可能性；二是增加纠错能力，使得即便出现差错，也能得到纠正。
4.1差错控制的基本概念
4.1.1差错的特点
由于通信线路上总有噪声存在，噪声和有用信息中的结果，就会出现差错。噪声可分为两类，一类是热噪声，另一类是冲激噪声，热噪声引起的差错是一种随机差错，亦即某个码元的出错具有独立性，与前后码元无关。冲激噪声是由短暂原因造成的，例如点击的启动、停止，电器设备的放弧等。冲击噪声引起的差错是成群的，其差错持续时间称为突发错的长度。
衡量信道传输性能的指标之一是误码率PO，PO=错误接收的码元数/接收的总码元数。目前普通电话线路中，当传输速率在600~2400bit/s时，PO在之间，对于大多数通信系统，PO在之间，而计算机之间的数据传输则要求误码率低于
。
4.1.2 差错控制的基本方式
差错控制方式基本上分为两类，一类称为“反馈纠错”，另一类称为“前向纠错”。在这两类基础上又派生出一种称为“混合纠错”。
(1)反馈纠错
这种方式在是发信端采用某种能发现一定程度传输差错的简单编码方法对所传信息进行编码，加入少量监督码元，在接收端则根据编码规则收到的编码信号进行检查，一量检测出(发现)有错码时，即向发信端发出询问的信号，要求重发。发信端收到询问信号时，立即重发已发生传输差错的那部分发信息，直到正确收到为止。所谓发现差错是指在若干接收码元中知道有一个或一些是错的，但不一定知道错误的准确位置。
(2)前向纠错
这种方式是发信端采用某种在解码时能纠正一定程度传输差错的较复杂的编码方法，使接收端在收到信码中不仅能发现错码，还能够纠正错码。采用前向纠错方式时，不需要反馈信道，也无需反复重发而延误传输时间，对实时传输有利，但是纠错设备比较复杂。
(3)混合纠错
混合纠错的方式是：少量纠错在接收端自动纠正，差错较严重，超出自行纠正能力时，就向发信端发出询问信号，要求重发。因此，“混合纠错”是“前向纠错”及“反馈纠错”两种方式的混合。
对于不同类型的信道，应采用不同的差错控制技术，否则就将事倍功半。
反馈纠错可用于双向数据通信，前向纠错则用于单向数字信号的传输，例如广播数字电视系统，因为这种系统没有反馈通道。
4.1.3 误码控制编码的分类
随着数字通信技术的发展，研究开发了各种误码控制编码方案，各自建立在不同的数学模型基础上，并具有不同的检错与纠错特性，可以从不同的角度对误码控制编码进行分类。
按照误码控制的不同功能，可分为检错码、纠错码和纠删码等。检错码仅具备识别错码功能而无纠正错码功能；纠错码不仅具备识别错码功能，同时具备纠正错码功能；纠删码则不仅具备识别错码和纠正错码的功能，而且当错码超过纠正范围时可把无法纠错的信息删除。
按照误码产生的原因不同，可分为纠正随机错误的码与纠正突发性错误的码。前者主要用于产生独立的局部误码的信道，而后者主要用于产生大面积的连续误码的情况，例如磁带数码记录中磁粉脱落而发生的信息丢失。按照信息码元与附加的监督码元之间的检验关系可分为线性码与非线性码。如果两者呈线性关系，即满足一组线性方程式，就称为线性码；否则，两者关系不能用线性方程式来描述，就称为非线性码。
按照信息码元与监督附加码元之间的约束方式之不同，可以分为分组码与卷积码。在分组码中，编码后的码元序列每n位分为一组，其中包括k位信息码元和r位附加监督码元，即n=k+r，每组的监督码元仅与本组的信息码元有关，而与其他组的信息码元无关。卷积码则不同，虽然编码后码元序列也划分为码组，但每组的监督码元不但与本组的信息码元有关，而且与前面码组的信息码元也有约束关系。
按照信息码元在编码之后是否保持原来的形式不变，又可分为系统码与非系统码。在系统码中，编码后的信息码元序列保持原样不变，而在非系统码中，信息码元会改变其原有的信号序列。由于原有码位发生了变化，使译码电路更为复杂，故较少选用。
根据编码过程中所选用的数字函数式或信息码元特性的不同，又包括多种编码方式。对于某种具体的数字设备，为了提高检错、纠错能力，通常同时选用几种误码控制编码方式。以下，以线性分组码为例，对几种简单的编码方式进行介绍。
4.2循环码（CRC）
循环码是线性码的一个重要的子类，它有以下两大特点：第一，码的结构可以用代数方法来构造和分析，并且可以找到各种实用的译码方法；第二，由于其循环特性，编码运算和伴随式计算，可用反馈移位寄存器来实现，硬件实现简单。
循环码是在严密的代数学理论基础上建立起来的，这种码的编码和解码设备都不太复杂，而且纠错的能力较强。顾名思义，循环码除具有线性码的一般性质之外，还具有循环性，即任一码组循环一位以后，仍为该码中的一个码组。
5.2 2FSK二进制数字调制的仿真
附录程序：标签：subplot,plot,码元,title,数字通信,系统,信号,ylabel From： https://blog.51cto.com/u_15815923/5744771
i=10; %基带信号码元数
j=5000;
a=round(rand(1,i)); %产生随机序列
t=linspace(0,5,j);
f1=10; %载波1频率
f2=5; %载波2频率
fm=i/5; %基带信号频率
B1=2*f1; %载波1带宽
B2=2*f2; %载波2带宽
%产生基带信号
st1=t;
for n=1:10
if a(n)<1;
for m=j/i*(n-1)+1:j/i*n
st1(m)=0;
end
else
for m=j/i*(n-1)+1:j/i*n
st1(m)=1;
end
end
end
st2=t;
%基带信号求反
for n=1:j;
if st1(n)>=1;
st2(n)=0;
else
st2(n)=1;
end
end;
figure(1);
subplot(411);
plot(t,st1);
title('基带信号');
axis([0,5,-1,2]);
subplot(412);
plot(t,st2);
title('基带信号反码');
axis([0,5,-1,2]);
%载波信号
s1=cos(2*pi*f1*t);
s2=cos(2*pi*f2*t);
subplot(413)
plot(s1);
title('载波信号1');
subplot(414),
plot(s2);
title('载波信号2');
%调制
F1=st1.*s1; %加入载波1
F2=st2.*s2; %加入载波2
figure(2);
subplot(311);
plot(t,F1); %s1*st1信号
title('s1*st1');
subplot(312);
plot(t,F2);
title('s2*st2'); %s2*st2信号
e_fsk=F1+F2;
subplot(313);
plot(t,e_fsk);
title('二进制频移键控调制信号') %二进制频移键控调制信号
4.3.2 调制解调仿真程序：
fs=2000; %采样频率
dt=1/fs;
f1=20;
f2=120; %两个信号的频率
a=round(rand(1,10)); %随机信号
g1=a
g2=~a; %信号反转，和g1反向
g11=(ones(1,2000))'*g1; %抽样
g1a=g11(:)';
g21=(ones(1,2000))'*g2;
g2a=g21(:)';
t=0:dt:10-dt;
t1=length(t);
fsk1=g1a.*cos(2*pi*f1.*t);
fsk2=g2a.*cos(2*pi*f2.*t);
fsk=fsk1+fsk2; %产生的信号
no=0.01*randn(1,t1); %噪声
sn=fsk+no;
subplot(311);
plot(t,no); %噪声波形
title('噪声波形')
ylabel('幅度')
subplot(312);
plot(t,fsk);
title('产生的波形')
ylabel('幅度')
subplot(313);
plot(t,sn);
title('将要通过滤波器的波形')
ylabel('幅度的大小')
xlabel('t')
figure(2) %FSK解调
b1=fir1(101,[10/800 20/800]);
b2=fir1(101,[90/800 110/800]); %设置带通参数
H1=filter(b1,1,sn);
H2=filter(b2,1,sn); %经过带通滤波器后的信号
subplot(211);
plot(t,H1);
title('经过带通滤波器f1后的波形')
ylabel('幅度')
subplot(212);
plot(t,H2);
title('经过带通滤波器f2后的波形')
ylabel('幅度')
xlabel('t')
sw1=H1.*H1;
sw2=H2.*H2; %经过相乘器
figure(3)
subplot(211);
plot(t,sw1);
title('经过相乘器h1后的波形')
ylabel('幅度')
subplot(212);
plot(t,sw2);
title('经过相乘器h2后的波形')
ylabel('•幅度')
xlabel('t')
bn=fir1(101,[2/800 10/800]); %经过低通滤波器
figure(4)
st1=filter(bn,1,sw1);
st2=filter(bn,1,sw2);
subplot(211);
plot(t,st1);
title('经过低通滤波器sw1后的波形')
ylabel('幅度')
subplot(212);
plot(t,st2);
title('经过低通滤波器sw2后的波形')
ylabel('幅度')
xlabel('t')
%判决
for i=1:length(t)
if(st1(i)>=st2(i))
st(i)=0;
else st(i)=st2(i);
end
end
figure(5)
st=st1+st2;
subplot(211);plot(t,st); title('经过抽样判决器后的波形')ylabel('幅度') subplot(212);plot(t,sn); title('原始的波形')
ylabel('幅度') xlabel('t')