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時的非均勻量化誤差圖')