用MATLAB程序开发设计PCM编码调制系统教学总结.docx

用MATLAB程序开发设计PCM编码调制系统
目录
一、作业要求 1
二、PCM编解码原理简述 1
三、MATLAB程序及仿真结果 6

一、作业要求
用MATLAB程序开发设计PCM编码调制系统。系统参数设置如下：
1． 输入模拟信源的最高频率限制在4KHz以内（可以是正弦信号，也可随机产生信号）；
2． 信号平均功率的动态范围≥30dB（自行设定幅度区间）。
要求完成以下工作：
1． 采用均匀量化，若要求信噪比不应低于25dB，不考虑线路衰减和损耗，则设计该均匀量化器，并显示原始输入信号，量化输出信号，以及编码结果。
2． 若采用A律13折线PCM编码，显示原始输入信号图及解码后的信号图；打印输出编码结果，并比较译码后的量化误差。

二、PCM编解码原理简述
均匀量化的基本原理
在脉冲编码调制中，模拟信号首先以高于奈奎斯特的速率采样，然后将所的样本量化。假设模拟信号是以［-Ⅹmax,Ⅹmax］表示的区间内分布的，而量化电平数很大。量化电平可以是相等的或是不相等的；前者就属于均匀PCM，而后者就是非均匀PCM。
关于量化的几个基本概念，量化间隔；量化误差；量化信噪比。
(1)相邻量化电平间距离称量化间隔， 用“Δ”表示。
(2)设抽样值为，量化后的值为， xq(kTs)与x(kTs)的误差称为量化误差，又称为量化噪声；量化误差不超过±Δ/2，而量化级数目越多，Δ值越小，量化误差也越小。
(3)衡量量化的性能好坏最常用指标是量化信噪比(Sq/Nq),其中Sq表示量化信号值xq(kTs)产生的功率，Nq表示量化误差功率，量化信噪比越大，则量化性能越好。
在均匀PCM中，长度为2Xmax的区间［-Ⅹmax,Ⅹmax］被划分为N个相等的子区间，每一子区间长度为△=2Xmax/N。如果N足够大，那么在每一子区间内输入的密度函数就能认为是均匀的，产生的失真为D=△2/12。
如果N是2的幂次方即,那么就要求用比特来表示每个量化电平。这就意味着，如果模拟信号的带宽是，采样又是在奈奎斯特率下完成的，那么传输PCM信号所要求的带宽至少是（实际上1.5比较接近于实际）。
这时失真由下式给出，

如果模拟信号的功率用表示，则信号/量化噪声的比（SQNR）由下式给出

式中表示归一化输入，定义为

以分贝（dB）计的SQNR为

量化以后,这些已量化的电平用比特对每个已量化电平进行编码.编码通常使用自然二进制码(NBC),即最低电平映射为全0序列,最高电平映射为全1序列,全部其余的电平按已量化值的递增次序映射。
非均匀量化的基本原理
正变换：，其中是归一化输入（），是一个参数，在标准律的非线性中它等于255。
反变换：
在非均匀量化PCM中，输入信号首先通过一非线性环节以减小输入的动态范围，再将输出加到某一均匀PCM系统上。在接收端，输出再通过另一非线形环节，该环节是在发送端所用的非线性环节的逆特征。这样，总的效果就等效于一个在量化电平之间具有非均匀间隔的PCM系统。
非均匀量化时，量化器随着输入信号的大小采用不同的量化间隔，大信号时采用大的量化间隔，小信号时采用小的量化间隔，可以以较少的量化电平数达到输入动态范围的要求
一般对语音信号传输来说，所使用的非线性可以是律的非线性，或是A律的非线性，中国和欧洲采用A率压缩特性(A=87.56)，北美和日本采用律压缩特性(=255)，压缩特性分别如下：
μ律：
A律：
式中，x为归一化输入，y为归一化输出，A、为压缩系数，上式也可以表示为，

由于A律在工程上不好实现，所以我们经常用近似的13折线压缩法去代替A率压缩
下表列出了计算值与13折线时的值的比较。
计算值与A律13折线时值的比较

0

1

0

1
按折线
分段时的
0

1
段落
1
2
3
4
5
6
7
8
斜率
16
16
8
4
2
1

表中第二行的值是根据时计算得到的，第三行的值是13折线分段时的值。可见，13折线各段落的分界点与曲线十分逼近，同时按2的幂次分割有利于数字化。
A律压扩特性是连续曲线，A律不同压扩特性也不同，在电路上实现这样的函数规律是相当复杂的。实际中，往往采用近似于A律函数规律的13折线（）的压扩特性。这样，它基本上保持了连续压扩特性曲线的优点，。本设计中所用到的PCM编码正是采用这种压扩特性来编码的。
在13折线法中，无论输入信号是正是负，均按8段折线（8个段落）进行编码。若用8位折叠二进制码来表示输入信号的抽样量化值，其中用第一位表示量化值的极性，其余七位（第二位至第八位）则表示抽样量化值的绝对大小。具体的做法是：用第二至第四位表示段落码，它的8种可能状态来分别代表8个段落的起点电平。其它四位表示段内码，它的16种可能状态来分别代表每一段落的16个均匀划分的量化级。这样处理的结果，8个段落被划分成128个量化级。段落码和8个段落之间的关系如表1所示；段内码与16个量化级之间的关系见表2所示。

表1 段落码 表2 段内码
段落序号
段落码
段落范围

量化间隔
段内码

量化间隔
段内码
8
111
1024-2048
15
1111
7
0111
7
110
512-1024
14
1110
6
0110
6
101
256-512
13
1101
5
0101
5
100
128-256
12
1100
4
0100
4
011
64-128
11
1011
3
0011
3
010
32-64
10
1010
2
0010
2
001
16-32
9
1001
1
0001
1
000
0-16
8
1000
0
0000
在13折线中，用8位的折叠二进制码表示信号量化值的具体步骤为：用第2到4位表示段落码，8个段落的起点电平由它的8种可能状态来分别表示。其他四位表示段内码，每一个段落它的16个均匀的划分量化级由它的16种可能状态来分别表示。这样就使得8个段落被划分为128个量化级。再加上负的，相当于一共有256种量化电平数。
解码时，码字电平=段落起始电平+8M5+4M6+2M7+M8*Δk
解码电平=码字电平+Δk2

三、MATLAB程序及仿真结果
注：
1、 为便于阅读，本节仅仅给出核心代码，详细代码请见本文附录。
2、 文中的图片若不够清晰，子文件夹中还有高清输出的图片（600ppi）

1、 生成题中要求的动态信号
题中要求输入模拟信源的最高频率限制在4KHz以内（可以是正弦信号，也可随机产生信号）；信号平均功率的动态范围≥30dB（自行设定幅度区间）。
为了能够适应不同条件，本文生成的是频率范围在0~4k之间变化，动态范围在30~35dB之间的随机信号。
fana=10^4;%analog freq
T=1/fana;%analog interval
SignalLength=64;
tana=(0:SignalLength-1)*T;
RealDynamicRange=0;
while(RealDynamicRange<=30||RealDynamicRange>=35)
while(RealDynamicRange<=30||RealDynamicRange>=35)
signal=randi([-300,300],1,SignalLength)/10;%raise accuracy
idx=find(signal==0); % find all 0
signal(idx)=1; % set 0.1 to these indexes
RealDynamicRange=20*log10(max(abs(signal))/ min(abs(signal)));
end
%check freq range-------------------
NFFT=2^nextpow2(SignalLength);
ini_fft=fft(signal,NFFT)/NFFT*2;%show reality amplitude
f=fana/NFFT*(0:1:NFFT-1);
abs_ini_fft=abs(ini_fft);
%limit band below 4k
k_stop_freq=floor(4000*NFFT*T);
ini_fft(:,k_stop_freq:NFFT-k_stop_freq)=0;%lpf
inisignal=real(ifft(ini_fft)*NFFT/2);%get the input signal
RealDynamicRange=20*log10(max(abs(inisignal))/ min(abs(inisignal)));
end
%check freq range again
inisignal_fft=fft(inisignal,NFFT)/NFFT*2;%show reality amplitude
abs_inisignal_fft=abs(inisignal_fft);
原始产生的随机信号，按照题目要求处理过之后的信号以及频谱如下图所示：

分别为原始随机信号、原始随机信号的带宽；处理之后带限在0~4kHz之内的信号，处理之后的信号的带宽。
生成信号的动态范围RealDynamicRange本次随机产生结果为34.8dB。
2、 信号采样
本节对上述产生的信号进行采样，采样的频率为10000Hz。
%%%%%%%%%%%%%%%%%%%%sampling%%%%%%%%%%%%%%%%%%%%%%%%%
fsam=10^5;%sampling freq
Tsam=1/fsam;
tsam=0:Tsam:max(tana);
sigsampling=interp1(tana,inisignal,tsam,'linear');
subplot(2,1,1);
plot(tana,inisignal);
title('analog signal');
xlabel('t/s');
ylabel('f(t)');
subplot(2,1,2);
stem(tsam,sigsampling,'Marker','none');
title('sampling signal');
xlabel('t/s');
ylabel('fs(t)');
print -r600 -djpeg AnalogAndSamplingSignal;
采样前后信号为

3、信号的均匀量化
%%%%%%%%%%%%%%%%%%%uniform quantization%%%%%%%%%%%%%%%%%%%%%%
N_uniquan=round((25-20*log10(min(abs(inisignal))/max(abs(inisignal)))-4.77)/6.02);
M_uniquan=2^N_uniquan;
sig_uniquan=quant(inisignal,M_uniquan,max(abs(inisignal)))';
plot(tana,inisignal,'b');
hold on;
plot(tana,sig_uniquan,'r');
print -r600 -djpeg UniQuan;
delta=2*(max(abs(sig_uniquan)))/(M_uniquan);
sig_uniquancode=round((sig_uniquan-min(sig_uniquan))/delta);
sig_uniquancode_s=dec2bin(sig_uniquancode,N_uniquan);
fid=fopen('sig_uniquancode_s.txt','w');
for i=1:length(inisignal)
fprintf(fid,'%s\r\n',sig_uniquancode_s(i,:));
end
fclose('all');
量化之前的信号与量化之后恢复的信号分别为，本次信号均匀量化所需的量化位数为9。量化信噪比小于25dB。均匀编码结果如附2所示。

4、PCM编解码
%%%%%%%%%%%%%%%%%%%%%%A_PCM%%%%%%%%%%%%%%%%%%%%%%%%
UniTo1=inisignal/max(abs(inisignal));
%jxm dlm dnm
jxm=zeros(1,length(UniTo1));
dlm=jxm;dnm=jxm;
quannumdlm=jxm;
for i=1:length(UniTo1)
%jxm
if(UniTo1(i))<0
jxm(i)=0;
else jxm(i)=1;
end
%dlm
tmpsig=abs(UniTo1(i));
quannumdlm(i)=round(tmpsig/(1/2048));
startquanV=[0,16,32,64,128,256,512,1024];
quaninterval=[1,1,2,4,8,16,32,64];
if(quannumdlm(i)<16)
dlm(i)=0;
elseif(quannumdlm(i)>=16&&quannumdlm(i)<32)
dlm(i)=1;
elseif(quannumdlm(i)>=32&&quannumdlm(i)<64)
dlm(i)=2;
elseif(quannumdlm(i)>=64&&quannumdlm(i)<128)
dlm(i)=3;
elseif(quannumdlm(i)>=128&&quannumdlm(i)<256)
dlm(i)=4;
elseif(quannumdlm(i)>=256&&quannumdlm(i)<512)
dlm(i)=5;
elseif(quannumdlm(i)>=512&&quannumdlm(i)<1024)
dlm(i)=6;
else
dlm(i)=7;
end
%dnm
dnm(i)=floor((quannumdlm(i)-(startquanV(dlm(i)+1)))/quaninterval(dlm(i)+1));
if(dnm(i)==16)
dnm(i)=15;
end
end
jxm_s=dec2bin(jxm);
dlm_s=dec2bin(dlm,3);
dnm_s=dec2bin(dnm,4);

fid=fopen('sig_pcmquancode_s.txt','w');
for i=1:length(jxm);
fprintf(fid,'%s%s%s\r\n',jxm_s(i,:),dlm_s(i,:),dnm_s(i,:));
end
fclose('all');

decodeV=zeros(1,length(inisignal));
recodeV=zeros(1,length(inisignal));
diffpcmquan=zeros(1,length(inisignal));
for i=1:length(inisignal)
decodeV(i)=startquanV(dlm(i)+1)+(dnm(i)+0.5)*quaninterval(dlm(i)+1);
diffpcmquan(i)=quannumdlm(i)-decodeV(i);
recodeV(i)=decodeV(i)/2048*max(abs(inisignal));
if(jxm(i)==0)
recodeV(i)=-recodeV(i);
end
end
采用了A律13折线PCM编解码，所得的原始信号、解码信号、量化误差如下图所示：
可见，很好地复原了原始信号，量化误差的最大值不超过30，编码结果在附3。

附1：完整程序
clc;clear;tic;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Éú³É¶¯Ì¬ÐźÅ-------------
fana=10^4;%analog freq
T=1/fana;%analog interval
SignalLength=64;
tana=(0:SignalLength-1)*T;
RealDynamicRange=0;
while(RealDynamicRange<=30||RealDynamicRange>=35)
while(RealDynamicRange<=30||RealDynamicRange>=35)
signal=randi([-300,300],1,SignalLength)/10;%raise accuracy
idx=find(signal==0); % find all 0
signal(idx)=1; % set 0.1 to these indexes
RealDynamicRange=20*log10(max(abs(signal))/ min(abs(signal)));
end
subplot(2,2,1);
plot(tana,signal);%for better visual eff
title('random signal');
xlabel('t/s');
ylabel('f(t)');
%saveas(gcf,'random signal.bmp','bmp');
%check freq range-------------------
NFFT=2^nextpow2(SignalLength);
ini_fft=fft(signal,NFFT)/NFFT*2;%show reality amplitude
f=fana/NFFT*(0:1:NFFT-1);
abs_ini_fft=abs(ini_fft);
subplot(2,2,2);
plot(f(1:NFFT/2),abs_ini_fft(1:NFFT/2));%for better visual eff
title('random signal band');
xlabel('f/hz');
ylabel('f(w)');
%saveas(gcf,'random signal band.bmp','bmp');
%limit band below 4k
k_stop_freq=floor(4000*NFFT*T);
ini_fft(:,k_stop_freq:NFFT-k_stop_freq)=0;%lpf
inisignal=real(ifft(ini_fft)*NFFT/2);%get the input signal
RealDynamicRange=20*log10(max(abs(inisignal))/ min(abs(inisignal)));
end
subplot(2,2,3);
plot(tana,inisignal);
title('initial signal');
xlabel('t/s');
ylabel('f(t)');
%saveas(gcf,'initial signal.bmp','bmp');
%check freq range again
inisignal_fft=fft(inisignal,NFFT)/NFFT*2;%show reality amplitude
abs_inisignal_fft=abs(inisignal_fft);
subplot(2,2,4);
plot(f(1:NFFT/2),abs_inisignal_fft(1:NFFT/2));%for better visual eff
title('initial signal band');
xlabel('f/hz');
ylabel('f(w)');
print -r600 -djpeg inputsignal;
%saveas(gcf,'initial signal band.bmp','bmp');
%check dynamic range again
%%%%%%%%%%%%%%%%%%%%sampling%%%%%%%%%%%%%%%%%%%%%%%%%
fsam=10^5;%sampling freq
Tsam=1/fsam;
tsam=0:Tsam:max(tana);
sigsampling=interp1(tana,inisignal,tsam,'linear');
%RealDynamicRange=20*log10(max(abs(sigsampling))/ min(abs(sigsampling)));
subplot(2,1,1);
plot(tana,inisignal);
title('analog signal');
xlabel('t/s');
ylabel('f(t)');
subplot(2,1,2);
stem(tsam,sigsampling,'Marker','none');
title('sampling signal');
xlabel('t/s');
ylabel('fs(t)');
print -r600 -djpeg AnalogAndSamplingSignal;
%%%%%%%%%%%%%%%%%%%uniform quantization%%%%%%%%%%%%%%%%%%%%%%
N_uniquan=round((25-20*log10(min(abs(inisignal))/max(abs(inisignal)))-4.77)/6.02);
M_uniquan=2^N_uniquan;
sig_uniquan=quant(inisignal,M_uniquan,max(abs(inisignal)))';
plot(tana,inisignal,'b');
hold on;
plot(tana,sig_uniquan,'r');
print -r600 -djpeg UniQuan;
delta=2*(max(abs(sig_uniquan)))/(M_uniquan);
sig_uniquancode=round((sig_uniquan-min(sig_uniquan))/delta);
sig_uniquancode_s=dec2bin(sig_uniquancode,N_uniquan);
fid=fopen('sig_uniquancode_s.txt','w');
for i=1:length(inisignal)
fprintf(fid,'%s\r\n',sig_uniquancode_s(i,:));
end
fclose('all');
%%%%%%%%%%%%%%%%%%%%%%A_PCM%%%%%%%%%%%%%%%%%%%%%%%%
UniTo1=inisignal/max(abs(inisignal));
%jxm dlm dnm
jxm=zeros(1,length(UniTo1));
dlm=jxm;dnm=jxm;
quannumdlm=jxm;
for i=1:length(UniTo1)
%jxm
if(UniTo1(i))<0
jxm(i)=0;
else jxm(i)=1;
end
%dlm
tmpsig=abs(UniTo1(i));
quannumdlm(i)=round(tmpsig/(1/2048));
startquanV=[0,16,32,64,128,256,512,1024];
quaninterval=[1,1,2,4,8,16,32,64];
if(quannumdlm(i)<16)
dlm(i)=0;
elseif(quannumdlm(i)>=16&&quannumdlm(i)<32)
dlm(i)=1;
elseif(quannumdlm(i)>=32&&quannumdlm(i)<64)
dlm(i)=2;
elseif(quannumdlm(i)>=64&&quannumdlm(i)<128)
dlm(i)=3;
elseif(quannumdlm(i)>=128&&quannumdlm(i)<256)
dlm(i)=4;
elseif(quannumdlm(i)>=256&&quannumdlm(i)<512)
dlm(i)=5;
elseif(quannumdlm(i)>=512&&quannumdlm(i)<1024)
dlm(i)=6;
else
dlm(i)=7;
end
%dnm
dnm(i)=floor((quannumdlm(i)-(startquanV(dlm(i)+1)))/quaninterval(dlm(i)+1));
if(dnm(i)==16)
dnm(i)=15;
end
end
jxm_s=dec2bin(jxm);
dlm_s=dec2bin(dlm,3);
dnm_s=dec2bin(dnm,4);

fid=fopen('sig_pcmquancode_s.txt','w');
for i=1:length(jxm);
fprintf(fid,'%s%s%s\r\n',jxm_s(i,:),dlm_s(i,:),dnm_s(i,:));
end
fclose('all');

decodeV=zeros(1,length(inisignal));
recodeV=zeros(1,length(inisignal));
diffpcmquan=zeros(1,length(inisignal));
for i=1:length(inisignal)
decodeV(i)=startquanV(dlm(i)+1)+(dnm(i)+0.5)*quaninterval(dlm(i)+1);
diffpcmquan(i)=quannumdlm(i)-decodeV(i);
recodeV(i)=decodeV(i)/2048*max(abs(inisignal));
if(jxm(i)==0)
recodeV(i)=-recodeV(i);
end
%diffpcmquan(i)=round(diffpcmquan(i)/5);
end
subplot(3,1,1);
plot(tana,inisignal,'b');
title('sampling signal');
xlabel('t/s');
ylabel('f(t)');
subplot(3,1,2);
plot(tana,recodeV,'r');
title('decode signal');
xlabel('t/s');
ylabel('f(t)');
subplot(3,1,3);
plot(tana,diffpcmquan,'g');
title('diffpcmquan');
xlabel('t/s');
ylabel('delta');
print -r600 -djpeg pcm_decode_diffquan;
toc;
附2：均匀编码结果（保存在sig_uniquancode_s.txt）
010110011
100011111
100100110
011100100
010110101
100001101
101111001
011000010
001101011
101001010
100100010
001011001
001111011
011000001
100000101
101000101
100100101
100011100
011000001
001110010
101100000
110111111
100101001
110000101
110000010
100011011
110001111
101010001
100011001
101011111
010111000
001111100
001110100
001011001
101100110
101110010
011001010
100101110
010001110
000011000
101001100
110001101
011110010
010101111
001011011
011100100
110111000
011110111
001110101
110010101
110111101
010111100
011101110
101111011
010111101
000000000
010101101
101100000
100001011
100111101
110110011
100001101
001100000
001100001
附3：PCM编码结果（保存在sig_pcmquancode_s.txt中）
01100011
11001110
11010010
01001011
01100010
10111010
11101110
01011110
01110010
11100010
11010000
01110100
01110000
01011111
10100010
11100001
11010010
11001011
01011111
01110001
11100111
11110111
11010100
11110000
11110000
11001010
11110001
11100100
11001000
11100111
01100010
01110000
01110001
01110100
11101001
11101100
01011010
11010110
01101100
01111100
11100010
11110001
00111100
01100100
01110100
01001011
11110110
00110001
01110001
11110010
11110111
01100000
01000001
11101110
01100000
01111111
01100100
11100111
10110110
11011110
11110110
10111001
01110100
01110011