% show the pcm encode and decode
clear all;
close all;

t=0:0.02:10;
vm1=-70:1:0;  % 輸入的正弦信號幅度不同
vm=10.^(vm1/20);
v=max(vm);
figure(1)
for k=1:length(vm)
        x=vm(k)*sin(2*pi*t);
        xx=x/v;  % normalize
        sxx=floor(xx*2048);
        y=pcm_encode(sxx);
        yy=pcm_decode(y,v);
        
        nq=sum((x-yy).*(x-yy))/length(x);
        sq=mean(yy.^2);
        snr=sq/nq;
        
        % drawnow
        subplot(211)
        plot(t,x);
        title('sample sequence');
        subplot(212);
        plot(t,yy);
        title('pcm decode sequence');
        snrq(k)=10*log10(snr);
end
figure(2)
plot(vm1,snrq);
axis([-60 0 0 60]);
xlabel('輸入正弦信號幅度的db值');
ylabel('量化信噪比')
grid;
function [out]=pcm_encode(x)
% x encode to pcm code
n=length(x);
% -2048<x<2048
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)=0;
        out(i,3)=0;
        out(i,4)=0;
        step=1;
        st=0;
    elseif 16<=abs(x(i)) & abs(x(i))<32
        out(i,2)=0;
        out(i,3)=0;
        out(i,4)=1;
        step=1;
        st=16;
    elseif 32<=abs(x(i)) & abs(x(i))<64
        out(i,2)=0;
        out(i,3)=1;
        out(i,4)=0;
        step=2;
        st=32;
    elseif 64<=abs(x(i)) & abs(x(i))<128
        out(i,2)=0;
        out(i,3)=1;
        out(i,4)=1;
        step=4;
        st=64;
    elseif 128<=abs(x(i)) & abs(x(i))<256
        out(i,2)=1;
        out(i,3)=0;
        out(i,4)=0;
        step=8;
        st=128;
    elseif 256<=abs(x(i)) & abs(x(i))<512
        out(i,2)=1;
        out(i,3)=0;
        out(i,4)=1;
        step=16;
        st=256;
    elseif 512<=abs(x(i)) & abs(x(i))<1024
        out(i,2)=1;
        out(i,3)=1;
        out(i,4)=0;
        step=32;
        st=512;
    elseif 1024<=abs(x(i)) & abs(x(i))<2048
        out(i,2)=1;
        out(i,3)=1;
        out(i,4)=1;
        step=64;
        st=1024;
    else
        out(i,2)=1;
        out(i,3)=1;
        out(i,4)=1;
        step=64;
        st=1024;
    end
    
    if(abs(x(i))>2048)
        out(i,2:8)=[1 1 1 1 1 1 1];
    else
        tmp=floor((abs(x(i))-st)/step);
        t=dec2bin(tmp,4)-48;  %  函數dec2bin輸出的是ASCII字符串，48對應0
        out(i,5:8)=t(1:4);
    end
end
out=reshape(out',1,8*n);