kamyk

close all;
clear all;
%初始化，n为量化级数
mu=input('Please enter the value MU:');
if isempty(mu)
    mu=255;
end;
n=input('Please enter the qualify level n:');
if isempty(n)
    n=8;
end;
dt=0.01;                                %时域分辨率的初始化（单位：ms）
N=500;                                  %采样点个数
T=dt*N;                                 %时域宽度
t=linspace(-T/2,T/2,N);                 %时间轴值
fx=0.04;                                %模拟信号的频率(单位:kHz)
a=sin(2*pi*fx*t);                       %绘制曲线(约能画出0.5个周期)
%a=randn(1,500);
%a=[20*randn(1,20),rand(1,480)];
a_ori=a;                                %保留原始的信号
lt=5;                                   %每次抽样需要处理的采样点个数
te=linspace(-T/2,T/2,N/lt);
for i=1:lt:length(t)
    if (i+lt)>length(t)
        a([i:length(a)])=a(i+1);        %防止超标
    else
        a(i+[0:lt])=a(i+1);
    end
end
amax=max(abs(a));
a_quan=a/amax;                          %归一化
a_quan=mulaw(a_quan,mu);                %进行mu率算法
b_quan=a_quan;
d=2/n;                                  %量化间隔
q=d.*[0:n-1];
q=q-((n-1)/2)*d;                        %量化电平
%进行量化
for i=1:n
    a_quan(find((q(i)-d/2<=a_quan) & (a_quan<=q(i)+d/2)))=...%定位第i个量化间隔码子
    q(i).*ones(1,length(find((q(i)-d/2<=a_quan)&(a_quan<=q(i)+d/2))));
    %赋值为相应的量化电平
    b_quan(find(a_quan==q(i))) =(i-1).*ones(1,length(find(a_quan==q(i))));%转化为对应的整数
end
a_quan=invmulaw(a_quan,mu);
a_quan=a_quan*amax;                     %恢复原值，去规一化（量化后）
nu=ceil(log2(n));                       %编码位数
code=zeros(length(a),nu);
for i=1:length(a)
    for j=nu:-1:0                       %从高向低编码
        if (fix(b_quan(i)/(2^j))==1)
        code(i,(nu-j))=1;
        b_quan(i)=b_quan(i)-2^j;
        end
    end
end

code_temp=zeros(10,nu);                 %存储前10个量化编码的中间变量
quan_temp=zeros(10,1);                  %存储前10个量化值的中间变量
count=1;                                %计数器
for i=1:length(a)/lt                    %循环找到前10个输入信号量化值
    count=count+1;
    quan_temp(i,1)=a_quan(1,i*lt);
    code_temp(i,:)=code(i*lt,:);
    if count==11                        %控制循环次数
        break;
    end
end

sqnr=20*log10(norm(a)/norm(a-a_quan))  %求量化信噪比
error=zeros(1,N/lt);
for i=1:N/lt
    error(i)=a(i*lt)-a_quan(i*lt);
end
%显示前10个值
quan_temp
code_temp
%画图:
figure(1)
set(1,'position',[10,50,350,250])%设置窗口位置
plot(t,a_quan,'b',t,a_ori,'r')
grid on
axis([-T/2,T/2,min(a),max(a)])%设置坐标格式
xlabel('t(ms)')%给X轴加标签
ylabel('s(t)') %给Y轴加标签
title('量化级数为n时的u率量化曲线')

figure(2)
set(2,'position',[350,50,350,250])%设置窗口位置
plot((a_ori/max(a_ori)),(a_quan/max(a_ori)),'b')
grid on
axis([-1,1,-1,1])%设置坐标格式
xlabel('input')%给X轴加标签
ylabel('output') %给Y轴加标签
title('量化级数为n时的输入－输出曲线')

figure(3)
set(3,'position',[700,50,350,250])%设置窗口位置
stem(te,error,'b')
grid on
axis([-T/2,T/2,min(a),max(a)])%设置坐标格式
xlabel('t(ms)')%给X轴加标签
ylabel('Δs') %给Y轴加标签
title('量化级数为n时的非均匀量化误差图')