用MATLAB程序开发设计PCM编码调制系统

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

2、。二、PCM编解码原理简述均匀量化的基本原理在脉冲编码调制中，模拟信号首先以高于奈奎斯特的速率采样，然后将所的样本量化。假设模拟信号是以-max,max表示的区间内分布的，而量化电平数很大。量化电平可以是相等的或是不相等的；前者就属于均匀PCM，而后者就是非均匀PCM。关于量化的几个基本概念，量化间隔；量化误差；量化信噪比。 (1)相邻量化电平间距离称量化间隔， 用“”表示。(2)设抽样值为，量化后的值为， xq(kTs)与x(kTs)的误差称为量化误差，又称为量化噪声；量化误差不超过±/2，而量化级数目越多，值越小，量化误差也越小。(3)衡量量化的性能好坏最常用指标是量化信噪比(S

3、q/Nq),其中Sq表示量化信号值xq(kTs)产生的功率，Nq表示量化误差功率，量化信噪比越大，则量化性能越好。在均匀PCM中，长度为2Xmax的区间-max,max被划分为N个相等的子区间，每一子区间长度为=2Xmax/N。如果N足够大，那么在每一子区间内输入的密度函数就能认为是均匀的，产生的失真为D=2/12。如果N是2的幂次方即,那么就要求用比特来表示每个量化电平。这就意味着，如果模拟信号的带宽是，采样又是在奈奎斯特率下完成的，那么传输PCM信号所要求的带宽至少是（实际上1.5比较接近于实际）。这时失真由下式给出，如果模拟信号的功率用表示，则信号/量化噪声的比（SQNR）由下式给出式中

4、表示归一化输入，定义为以分贝（dB）计的SQNR为量化以后,这些已量化的电平用比特对每个已量化电平进行编码.编码通常使用自然二进制码(NBC),即最低电平映射为全0序列,最高电平映射为全1序列,全部其余的电平按已量化值的递增次序映射。非均匀量化的基本原理正变换：，其中是归一化输入（），是一个参数，在标准律的非线性中它等于255。反变换：在非均匀量化PCM中，输入信号首先通过一非线性环节以减小输入的动态范围，再将输出加到某一均匀PCM系统上。在接收端，输出再通过另一非线形环节，该环节是在发送端所用的非线性环节的逆特征。这样，总的效果就等效于一个在量化电平之间具有非均匀间隔的PCM系统。非均匀量化

5、时，量化器随着输入信号的大小采用不同的量化间隔，大信号时采用大的量化间隔，小信号时采用小的量化间隔，可以以较少的量化电平数达到输入动态范围的要求一般对语音信号传输来说，所使用的非线性可以是律的非线性，或是A律的非线性，中国和欧洲采用A率压缩特性(A=87.56)，北美和日本采用律压缩特性(=255)，压缩特性分别如下：律：A律：式中，x为归一化输入，y为归一化输出，A、为压缩系数，上式也可以表示为，由于A律在工程上不好实现，所以我们经常用近似的13折线压缩法去代替A率压缩下表列出了计算值与13折线时的值的比较。计算值与A律13折线时值的比较0101按折线分段时的01段落12345678斜率16

6、168421表中第二行的值是根据时计算得到的，第三行的值是13折线分段时的值。可见，13折线各段落的分界点与曲线十分逼近，同时按2的幂次分割有利于数字化。A律压扩特性是连续曲线，A律不同压扩特性也不同，在电路上实现这样的函数规律是相当复杂的。实际中，往往采用近似于A律函数规律的13折线（）的压扩特性。这样，它基本上保持了连续压扩特性曲线的优点，。本设计中所用到的PCM编码正是采用这种压扩特性来编码的。在13折线法中，无论输入信号是正是负，均按8段折线（8个段落）进行编码。若用8位折叠二进制码来表示输入信号的抽样量化值，其中用第一位表示量化值的极性，其余七位（第二位至第八位）则表示抽样量化值的绝

7、对大小。具体的做法是：用第二至第四位表示段落码，它的8种可能状态来分别代表8个段落的起点电平。其它四位表示段内码，它的16种可能状态来分别代表每一段落的16个均匀划分的量化级。这样处理的结果，8个段落被划分成128个量化级。段落码和8个段落之间的关系如表1所示；段内码与16个量化级之间的关系见表2所示。表1 段落码 表2 段内码段落序号段落码段落范围量化间隔段内码量化间隔段内码81111024-2048151111701117110512-1024141110601106101256-512131101501015100128-25612110040100401164-128111011300

8、11301032-6410101020010200116-32910011000110000-168100000000在13折线中，用8位的折叠二进制码表示信号量化值的具体步骤为：用第2到4位表示段落码，8个段落的起点电平由它的8种可能状态来分别表示。其他四位表示段内码，每一个段落它的16个均匀的划分量化级由它的16种可能状态来分别表示。这样就使得8个段落被划分为128个量化级。再加上负的，相当于一共有256种量化电平数。解码时，码字电平=段落起始电平+8M5+4M6+2M7+M8*k 解码电平=码字电平+k2三、MATLAB程序及仿真结果注：1、 为便于阅读，本节仅仅给出核心代码，详细代码请

9、见本文附录。2、 文中的图片若不够清晰，子文件夹中还有高清输出的图片（600ppi）1、 生成题中要求的动态信号题中要求输入模拟信源的最高频率限制在4KHz以内（可以是正弦信号，也可随机产生信号）；信号平均功率的动态范围30dB（自行设定幅度区间）。为了能够适应不同条件，本文生成的是频率范围在04k之间变化，动态范围在3035dB之间的随机信号。fana=104;%analog freqT=1/fana;%analog intervalSignalLength=64;tana=(0:SignalLength-1)*T;RealDynamicRange=0;while(RealDynamicRa

10、nge<=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 ra

11、nge-NFFT=2nextpow2(SignalLength);ini_fft=fft(signal,NFFT)/NFFT*2;%show reality amplitudef=fana/NFFT*(0:1:NFFT-1);abs_ini_fft=abs(ini_fft);%limit band below 4kk_stop_freq=floor(4000*NFFT*T);ini_fft(:,k_stop_freq:NFFT-k_stop_freq)=0;%lpfinisignal=real(ifft(ini_fft)*NFFT/2);%get the input signalRealDyn

12、amicRange=20*log10(max(abs(inisignal)/ min(abs(inisignal);end%check freq range againinisignal_fft=fft(inisignal,NFFT)/NFFT*2;%show reality amplitudeabs_inisignal_fft=abs(inisignal_fft);原始产生的随机信号，按照题目要求处理过之后的信号以及频谱如下图所示：分别为原始随机信号、原始随机信号的带宽；处理之后带限在04kHz之内的信号，处理之后的信号的带宽。生成信号的动态范围RealDynamicRange本次随机产生结

13、果为34.8dB。2、 信号采样本节对上述产生的信号进行采样，采样的频率为10000Hz。%sampling%fsam=105;%sampling freqTsam=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);ste

14、m(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=

15、2N_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

16、(sig_uniquancode,N_uniquan);fid=fopen('sig_uniquancode_s.txt','w');for i=1:length(inisignal) fprintf(fid,'%srn',sig_uniquancode_s(i,:);endfclose('all');量化之前的信号与量化之后恢复的信号分别为，本次信号均匀量化所需的量化位数为9。量化信噪比小于25dB。均匀编码结果如附2所示。4、PCM编解码%A_PCM%UniTo1=inisignal/max(abs(inisignal);%j

17、xm dlm dnmjxm=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)<

18、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(quan

19、numdlm(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; endendjxm_s=dec2bin(jxm);dlm_s=dec2bin(dlm,3);dn

20、m_s=dec2bin(dnm,4);fid=fopen('sig_pcmquancode_s.txt','w');for i=1:length(jxm); fprintf(fid,'%s%s%srn',jxm_s(i,:),dlm_s(i,:),dnm_s(i,:);endfclose('all');decodeV=zeros(1,length(inisignal);recodeV=zeros(1,length(inisignal);diffpcmquan=zeros(1,length(inisignal);for i=1:le

21、ngth(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); endend采用了A律13折线PCM编解码，所得的原始信号、解码信号、量化误差如下图所示：可见，很好地复原了原始信号，量化误差的最大值不超过30，编码结果在附3。附1：完整程序clc

22、;clear;tic;%Éú³É¶¯Ì¬ÐźÅ-fana=104;%analog freqT=1/fana;%analog intervalSignalLength=64;tana=(0:SignalLength-1)*T;RealDynamicRange=0;while(RealDynamicRange<=30|RealDynamicRange>=35) while(RealDynamicRange<=30|RealDynamicRange>=35

23、) 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); endsubplot(2,2,1);plot(tana,signal);%for better visual efftitle('random signal');xlabel('t/s&#

24、39;);ylabel('f(t)');%saveas(gcf,'random signal.bmp','bmp');%check freq range-NFFT=2nextpow2(SignalLength);ini_fft=fft(signal,NFFT)/NFFT*2;%show reality amplitudef=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

25、 visual efftitle('random signal band');xlabel('f/hz');ylabel('f(w)');%saveas(gcf,'random signal band.bmp','bmp');%limit band below 4kk_stop_freq=floor(4000*NFFT*T);ini_fft(:,k_stop_freq:NFFT-k_stop_freq)=0;%lpfinisignal=real(ifft(ini_fft)*NFFT/2);%get the inpu

26、t signalRealDynamicRange=20*log10(max(abs(inisignal)/ min(abs(inisignal);endsubplot(2,2,3);plot(tana,inisignal);title('initial signal');xlabel('t/s');ylabel('f(t)');%saveas(gcf,'initial signal.bmp','bmp');%check freq range againinisignal_fft=fft(inisignal,NFFT

27、)/NFFT*2;%show reality amplitudeabs_inisignal_fft=abs(inisignal_fft);subplot(2,2,4);plot(f(1:NFFT/2),abs_inisignal_fft(1:NFFT/2);%for better visual efftitle('initial signal band');xlabel('f/hz');ylabel('f(w)');print -r600 -djpeg inputsignal;%saveas(gcf,'initial signal ban

28、d.bmp','bmp');%check dynamic range again%sampling%fsam=105;%sampling freqTsam=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

29、 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

30、(min(abs(inisignal)/max(abs(inisignal)-4.77)/6.02);M_uniquan=2N_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(

31、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,'%srn',sig_uniquancode_s(i,:);endfclose('all');%A_PCM%UniTo1=inisignal/max(abs(inisignal);%jxm dlm dnm

32、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)

33、=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)&

34、gt;=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; endendjxm_s=dec2bin(jxm);dlm_s=dec2bin(dlm,3);dnm_s=dec2bi

35、n(dnm,4);fid=fopen('sig_pcmquancode_s.txt','w');for i=1:length(jxm); fprintf(fid,'%s%s%srn',jxm_s(i,:),dlm_s(i,:),dnm_s(i,:);endfclose('all');decodeV=zeros(1,length(inisignal);recodeV=zeros(1,length(inisignal);diffpcmquan=zeros(1,length(inisignal);for i=1:length(inisi

36、gnal) 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);endsubplot(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&#