homework2_2.m

该代码实现了A率和mu率PCM算法

%mu律的非均匀量化过程图、量化误差图及其余相关参数
close all
clear all
k=input('采样点数N=2^k，请输入k的值，默认为6:');%设置采样点数
if isempty(k)
    k=6;
end
mu=255;%设定mu律的参数mu为255
N=2^k;
%非均匀量化级数都为256
n1=8;
n2=16;
n3=64;
s=linspace(0,2*pi,N);%产生采样点
a=2*sin(s);%产生输入信号
[f1,amax]=mulaw(a,mu);%对输入信号a进行mu律压缩
[f2,amax]=mulaw(a,mu);%对输入信号a进行mu律压缩
[f3,amax]=mulaw(a,mu);%对输入信号a进行mu律压缩
%对压缩结果进行归一化
f1_max=max(abs(f1));
f2_max=max(abs(f2));
f3_max=max(abs(f3));
f1_quan=f1/f1_max;
f2_quan=f2/f2_max;
f3_quan=f1/f3_max;
b1_quan=f1_quan;
b2_quan=f2_quan;
b3_quan=f3_quan;
d1=2/n1;%量化间隔
d2=2/n2;%量化间隔
d3=2/n3;%量化间隔
q1=d1.*[0:n1-1];
q1=q1-((n1-1)/2)*d1;%设置量化电平
q2=d2.*[0:n2-1];
q2=q2-((n2-1)/2)*d2;%设置量化电平
q3=d3.*[0:n3-1];
q3=q3-((n3-1)/2)*d3;%设置量化电平
for i=1:n1
    f1_quan(find((q1(i)-d1/2<=f1_quan) & (f1_quan<=q1(i)+d1/2)))=...%定位第i个量化间隔码子
    q1(i).*ones(1,length(find((q1(i)-d1/2<=f1_quan)&(f1_quan<=q1(i)+d1/2))));%赋值为相应的量化电平
    b1_quan(find(f1_quan==q1(i))) = (i-1) .* ones(1,length(find(f1_quan==q1(i))));%对量化电平进行编号
end
for i=1:n2
    f2_quan(find((q2(i)-d2/2<=f2_quan) & (f2_quan<=q2(i)+d2/2)))=...%定位第i个量化间隔码子
    q2(i).*ones(1,length(find((q2(i)-d2/2<=f2_quan)&(f2_quan<=q2(i)+d2/2))));%赋值为相应的量化电平
    b2_quan(find(f2_quan==q2(i))) = (i-1) .* ones(1,length(find(f2_quan==q2(i))));%对量化电平进行编号
end
for i=1:n3
    f3_quan(find((q3(i)-d3/2<=f3_quan) & (f3_quan<=q3(i)+d3/2)))=...%定位第i个量化间隔码子
    q3(i).*ones(1,length(find((q3(i)-d3/2<=f3_quan)&(f3_quan<=q3(i)+d3/2))));%赋值为相应的量化电平
    b3_quan(find(f3_quan==q3(i))) = (i-1) .* ones(1,length(find(f3_quan==q3(i))));%对量化电平进行编号
end
%恢复原值（量化后）
f1_quan=f1_quan*f1_max;
f2_quan=f2_quan*f2_max;
f3_quan=f3_quan*f3_max;
a1_quan=invmulaw(f1_quan,mu);%对量化完毕的信号进行mu律的扩张
a2_quan=invmulaw(f2_quan,mu);%对量化完毕的信号进行mu律的扩张
a3_quan=invmulaw(f3_quan,mu);%对量化完毕的信号进行mu律的扩张
%由于扩张后得到的结果是对原信号的归一化，因此还需要再次进行恢复
a1_quan=a1_quan*amax;
a2_quan=a2_quan*amax;
a3_quan=a3_quan*amax;
nu1=ceil(log2(n1));%向上取整，得到码字的长度
nu2=ceil(log2(n2));%向上取整，得到码字的长度
nu3=ceil(log2(n3));%向上取整，得到码字的长度
%初始化码序列为全0矩阵
code1=zeros(length(a),nu1);
code2=zeros(length(a),nu2);
code3=zeros(length(a),nu3);
for i=1:length(a)
    %由于非均匀量化编码为折叠码，因此在编码之前要对量化电平的编号进行一定的处理
    if b1_quan(i)<n1/2
        b1_quan(i)=n1/2-b1_quan(i)-1;
    end
    if b2_quan(i)<n2/2
        b2_quan(i)=n2/2-b2_quan(i)-1;
    end
    if b3_quan(i)<n3/2
        b3_quan(i)=n3/2-b3_quan(i)-1;
    end
    for j=nu1:-1:0%从高向低编码
        if fix(b1_quan(i)/(2^j))==1
            code1(i,(nu1-j))=1;
            b1_quan(i)=b1_quan(i)-2^j;
        end
    end
    for j=nu2:-1:0%从高向低编码
        if fix(b2_quan(i)/(2^j))==1
            code2(i,(nu2-j))=1;
            b2_quan(i)=b2_quan(i)-2^j;
        end
    end
    for j=nu3:-1:0%从高向低编码
        if fix(b3_quan(i)/(2^j))==1
            code3(i,(nu3-j))=1;
            b3_quan(i)=b3_quan(i)-2^j;
        end
    end
end
%求量化信噪比
sqnr1=20*log10(norm(a)/norm(a-a1_quan));%norm()为求欧式距离的函数
sqnr2=20*log10(norm(a)/norm(a-a2_quan));
sqnr3=20*log10(norm(a)/norm(a-a3_quan));
input('下面将显示\mu律下量化级数为8时前十个输入信号的量化值，请按回车键继续:')
a1_quan(1:10)
input('下面将显示\mu律下量化级数为8时前十个输入信号的编码，请按回车键继续:')
code1(1:10,:)
input('下面将显示\mu律下量化级数为16时前十个输入信号的量化值，请按回车键继续:')
a2_quan(1:10)
input('下面将显示\mu律下量化级数为16时前十个输入信号的编码，请按回车键继续:')
code2(1:10,:)
input('下面将显示\mu律下量化级数为64时前十个输入信号的量化值，请按回车键继续:')
a3_quan(1:10)
input('下面将显示\mu律下量化级数为64时前十个输入信号的编码，请按回车键继续:')
code3(1:10,:)
input('下面将显示图线，请按回车键继续:')
figure(1)%开启一个编号为1的绘图窗口
set(1,'position',[0,400,340,300])%设置窗口1的位置
hold on%使得MATLAB不清除图形窗口，从而可以在同一图形窗口内画多条图线
grid on%开启网格
plot(s,a,'k')%画原信号的图线
plot(s,a,'.r')%画出采样点的位置
stairs(s,a1_quan)%画出量化后信号的阶梯图
axis([0,2*pi,-2,2])%设置坐标轴的观察范围
text(0.5,-0.4,'信噪比为:','Fontsize',10.5)
text(0.5,-0.7,num2str(sqnr1),'Fontsize',10.5)
title('\mu律非均匀量化过程曲线图(量化级数为8)')
xlabel('t(ms)')
ylabel('幅度值(V)')
figure(2)%开启一个编号为2的绘图窗口
set(2,'position',[340,400,340,300])%设置窗口2的位置
hold on%使得MATLAB不清除图形窗口，从而可以在同一图形窗口内画多条图线
grid on%开启网格
plot(s,a,'k')%画原信号的图线
plot(s,a,'.r')%画出采样点的位置
stairs(s,a2_quan)%画出量化后信号的阶梯图
axis([0,2*pi,-2,2])%设置坐标轴的观察范围
text(0.5,-0.4,'信噪比为:','Fontsize',10.5)
text(0.5,-0.7,num2str(sqnr2),'Fontsize',10.5)
title('\mu律非均匀量化过程曲线图(量化级数为16)')
xlabel('t(ms)')
ylabel('幅度值(V)')
figure(3)%开启一个编号为3的绘图窗口
set(3,'position',[680,400,340,300])%设置窗口1的位置
hold on%使得MATLAB不清除图形窗口，从而可以在同一图形窗口内画多条图线
grid on%开启网格
plot(s,a,'k')%画原信号的图线
plot(s,a,'.r')%画出采样点的位置
stairs(s,a3_quan)%画出量化后信号的阶梯图
axis([0,2*pi,-2,2])%设置坐标轴的观察范围
text(0.5,-0.4,'信噪比为:','Fontsize',10.5)
text(0.5,-0.7,num2str(sqnr3),'Fontsize',10.5)
title('\mu律非均匀量化过程曲线图(量化级数为64)')
xlabel('t(ms)')
ylabel('幅度值(V)')
figure(4)%开启一个编号为4的绘图窗口
set(4,'position',[0,50,340,300])%设置窗口3的位置
grid on%开启网格
stem(s,a-a1_quan)%画出量化信噪比的离散序列数据柄状图
axis([0,2*pi,min(a-a1_quan),max(a-a1_quan)])%设置坐标轴的观察范围
title('\mu律非均匀量化误差曲线图(量化级数为8)')
xlabel('t(ms)')
ylabel('幅度值(V)')
figure(5)%开启一个编号为5的绘图窗口
set(5,'position',[340,50,340,300])%设置窗口3的位置
grid on%开启网格
stem(s,a-a2_quan)%画出量化信噪比的离散序列数据柄状图
axis([0,2*pi,min(a-a2_quan),max(a-a2_quan)])%设置坐标轴的观察范围
title('\mu律非均匀量化误差曲线图(量化级数为16)')
xlabel('t(ms)')
ylabel('幅度值(V)')
figure(6)%开启一个编号为6的绘图窗口
set(6,'position',[680,50,340,300])%设置窗口4的位置
grid on%开启网格
stem(s,a-a3_quan)%画出量化信噪比的离散序列数据柄状图
axis([0,2*pi,min(a-a3_quan),max(a-a3_quan)])%设置坐标轴的观察范围
title('\mu律非均匀量化误差曲线图(量化级数为64)')
xlabel('t(ms)')
ylabel('幅度值(V)')