语音信号PCM编译码 MATLAB程序

% clc

% clear all

% 产生信号

load handel

x=y(1:20000)%取前20000个采样点

sound(x,Fs)

%PCM编码

x1=x/0.8.*2048

yy=pcm_encode(x1)

figure

subplot(2,1,1)

stem(yy(1:80),'.')

title('PCM编码后的波形')

%加噪声

snr=10

sp=mean(yy.^2)

attn=sp./ 10^(snr/10)

attn = sqrt(attn)

noise=randn(1,length(yy)).*attn

np=mean(noise.^2)

snr1=10*log10(sp/np)

data=yy+noise

% data=yy%不加噪声

subplot(2,1,2)

stem(data(1:80),'.r')

title('PCM加噪声后波形')

%译码

demodata=data>0.5

zz=pcm_decode(demodata,0.8)

figure

subplot(2,1,1)

plot(x)

title('原始语音信号')

subplot(2,1,2)

plot(zz)

title('译码后的语音信号')

sound(zz,Fs)

figure

plot(x,'b')

hold on

plot(zz,'r')

legend('编译码前的语音信号','编译码后的语音信号')

title('编，译码前后的语音信号')

function y = pcm_encode( x )

y=zeros(length(x),8) %存储矩阵（全零）

z=sign(x) %判断x的正负

x=abs(x)%取绝对值

%%段落码判断段区间的取值范围为前开后闭区间

for k=1:length(x)

%符号位的判断

if z(k)>0

y(k,1)=1

elseif z(k)<0

y(k,1)=0

end

if x(k)>128 &x(k)<=2048 %在第五段与第八段之间，段位码第一位都为“1”

y(k,2)=1

end

if (x(k)>32 &x(k)<=128) || (x(k)>512 &x(k)<=2048)

y(k,3)=1 %在第三四七八段内，段位码第二位为“1”

end

if (x(k)>16&x(k)<=32)||(x(k)>64&x(k)<=128)||(x(k)>256&x(k)<=512)||(x(k)>1024&x(k)<=2048)

y(k,4)=1 %在二四六八段内，段位码第三位为“1”

end

end

%段内码判断程序

N=zeros(1,length(x))

for k=1:length(x)

N(k)=y(k,2)*4+y(k,3)*2+y(k,4)+1 %找到x位于第几段

end

a=[0,16,32,64,128,256,512,1024] %量化间隔

b=[1,1,2,4,8,16,32,64]%除以16，得到每段的最小量化间隔

for m=1:length(x)

q=ceil((x(m)-a(N(m)))/b(N(m))) %求出在段内的位置

if q==0

y(m,(5:8))=[0,0,0,0] %如果输入为零则输出“0”

else k=num2str(dec2bin(q-1,4)) %编码段内码为二进制

y(m,5)=str2num(k(1))

y(m,6)=str2num(k(2))

y(m,7)=str2num(k(3))

y(m,8)=str2num(k(4))

end

end

%将N行8列矩阵转换为1行8*N列的矩阵

y=y'

y=reshape(y,1,length(x)*8)

end

function x=pcm_decode(y,max)

%将1行8*N列的矩阵转换为N行8列矩阵

y=reshape(y,8,length(y)/8)

y=y'

%PCM译码

n=size(y,1) %求出输入码组的个数

a=[0,16,32,64,128,256,512,1024] %段落起点值

b=[1,1,2,4,8,16,32,64] %每段的最小量化间隔

for k=1:n

t1=y(k,1)%取符号

t2=y(k,2)*4+y(k,3)*2+y(k,4)+1%判断段落位置

t3=y(k,5)*8+y(k,6)*4+y(k,7)*2+y(k,8) %判断段内位置

if t3==0 %段内码为零时

m(k)=(a(t2)+1+0.5*b(t2))/2048*max

else

m(k)=(a(t2)+b(t2)*t3+0.5*b(t2))/2048*max%还原出量化后的电平值

end

%判断符号位

if t1==0

x(k)=-m(k)

else

x(k)=m(k)

end

end

end

function [f] = ecode(yn)

%本函数实现国际通用的PCM量化A律13折线特性近似，x为输入的序列，变换后的值赋给序列y

x = yn

ypcm=[]

x=x/max(x)%求出序列的最大值，同时归一化

s=sign(x) %求得每一序列值的符号

x=abs(x) %求序列的绝对值

for i=1:length(x)%直接将序列的绝对值量化

if((x(i)>=0)&&(x(i)<1/64))%序列值位于第1和第2折线

ypcm(i)=16*x(i)

else if(x(i)>=1/64&&x(i)<1/32) %序列值位于第3折线

ypcm(i)=8*x(i)+1/8

else if(x(i)>=1/32&&x(i)<1/16) %序列值位于第4折线

ypcm(i)=4*x(i)+2/8

else if(x(i)>=1/16&&x(i)<1/8)%序列值位于第5折线

ypcm(i)=2*x(i)+3/8

else if(x(i)>=1/8&&x(i)<1/4) %序列值位于第6折线

ypcm(i)=x(i)+4/8

else if(x(i)>=1/4&&x(i)<1/2) %序列值位于第7折线

ypcm(i)=1/2*x(i)+5/8

else if(x(i)>=1/2&&x(i)<=1)%序列值位于第8折线

ypcm(i)=1/4*x(i)+6/8

end

end

end

end

end

end

end

end

ypcm= s.*ypcm%重新将符号带回序列中

%本函数实现将输入的值（已量化好）编码，输入y为量化后的序列，其值应在1到1之间

f=zeros(length(ypcm),8) %定义一个二维数组，每一行的8位代表了对应的输入值

s=sign(ypcm)%得到输入序列的符号，确定编码的首位

ypcm=ypcm*128 %将序列值扩展到0到128之间，便于编码

ypcm=fix(ypcm) %将计算值取整

ypcm=abs(ypcm) %只计算绝对值的编码

for i=1:length(ypcm)

if(ypcm(i)==128) %如果输入为1，得到128，避免出现编码位为2的错误

ypcm(i)=127.999 %将其值近似为127.999

end

end

for i=1:length(ypcm)%下面的循环是将十进制转换为二进制数

for j=6:-1:0%分别计算序列值除以从64到1的数的商

f(i,8-j)=fix(ypcm(i)/(2^j))

ypcm(i)=mod(ypcm(i),(2^j))

end

end

for i=1:length(ypcm)

if(s(i)==-1)%输入值是负数

f(i,1)=0 %首位取0

else

f(i,1)=1 %输入值是正数，首位取1

end

end

%显示编码结果

我修改了一下，能够运行了，不知道是不是你想要的结果？

clear all

close all

t=0:0.01:10 %定义时间抽样点

%vm1=-70:1:0 %输入的信号幅度的db值

vm1=linspace(-70,0,1001)

vm=10.^(vm1/20) %输入信号幅度

figure(1)

for k=1:length(vm)

for m=1:2

%x=vm*sin(2*pi*t+2*pi*rand(t)) %输入语音信号

x=vm.*sin(2*pi*t) %输入语音信号

v=1

xx=x/v

sxx=floor(xx*4096)

y=pcm_encode(sxx) %PCM编码

yy=pcm_decode(y,v)%PCM译码后信号幅值

nq(m)=sum((x-yy).*(x-yy))/length(x)%噪音功率

sq(m)=mean(yy.^2)%信号均值

snr(m)=(sq(m)/nq(m)) %信噪比

drawnow

subplot(211)

plot(t,x)

title('采样序列')%画出采样序列的图形

subplot(212)

plot(t,yy)

title('解码序列')%画出PCM解码后的序列图

end

snrq(k)=10*log10(mean(snr)) %量化信噪比

end

figure(2)

plot(vm1,snrq)

axis([-60 0 0 60]) %X轴范围是（-60，0）Y轴范围是（0,60）

grid

function out=pcm_decode(in,v)

n=length(in)

in=reshape(in',8,n/8)'%将in值变换成8行

slot(1)=0

slot(2)=16

slot(3)=32

slot(4)=64

slot(5)=128

slot(6)=256

slot(7)=512

slot(8)=1024

step(1)=1

step(2)=1

step(3)=2

step(4)=4

step(5)=8

step(6)=16

step(7)=32

step(8)=64

for i=1:n/8

ss=2*in(i,1)-1

tmp=in(i,2)*4+in(i,3)*2+in(i,4)+1

st=slot(tmp)

dt=(in(i,5)*8+in(i,6)*4+in(i,7)*2+in(i,8))*step(tmp)+0.5*step(tmp)

out(i)=ss*(st+dt)/4096*v%量化输出值

end

function [out]=pcm_encode(x) %定义A率13折线压缩特性

n=length(x)

for i=1:n

if x(i)>0

out(i,1)=1 %代表正值

else

out(i,1)=0 %代表负值

end

if abs(x(i))>=0&&abs(x(i))<16

out(i,2)=0out(i,3)=0out(i,4)=0step=1st=0

elseif 16<=abs(x(i))&&abs(x(i))<32

out(i,2)=0out(i,3)=0out(i,4)=1step=1st=16

elseif 32<=abs(x(i))&&abs(x(i))<64

out(i,2)=0out(i,3)=1out(i,4)=0step=2st=32

elseif 64<=abs(x(i))&&abs(x(i))<128

out(i,2)=0out(i,3)=1out(i,4)=1step=4st=64

elseif 128<=abs(x(i))&&abs(x(i))<256

out(i,2)=1out(i,3)=0out(i,4)=0step=8st=128

elseif 256<=abs(x(i))&&abs(x(i))<512

out(i,2)=1out(i,3)=0out(i,4)=1step=16st=256

elseif 512<=abs(x(i))&&abs(x(i))<1024

out(i,2)=1out(i,3)=1out(i,4)=0step=32st=512

elseif 1024<=abs(x(i))&&abs(x(i))<2048

out(i,2)=1out(i,3)=1out(i,4)=1step=64st=1024 %由抽样值定义段落编码

end

if(abs(x(i))>=2048)

out(i,2:8)=[1 1 1 1 1 1 1]step=128st=2048

else

tmp=floor((abs(x(i))-st)/step)

if tmp<0

a=1

end

t=dec2bin(tmp,4)-48 %dec2bin函数表示输出ASCII字符串值，48表示0

out(i,5:8)=t(1:4)%输出段内码

end

end

out=reshape(out',1,8*n) %将out值变换成1行8n列

---