http://s10/middle/59a069d5naf05485ecfc9&690

% 【例题6-1】设低通信号x(t)=0.1sin2πt +0.5cos4πt 。
% 1、 画出该低通信号的波形；
% 2、 画出抽样速率为 =4Hz的抽样序列；
% 3、 抽样序列恢复出原始信号。
% 低通抽样定理 dtchy.m

clear all;close all;
dt=0.01;
t=0:dt:10;
xt=0.1*sin(2*pi*t)+0.5*cos(4*pi*t);
[f,xf]=FFT_SHIFT(t,xt);
fs=4;
sdt=1/fs;
t1=0:sdt:10;
st=0.1*sin(2*pi*t1)+0.5*cos(4*pi*t1);
[f1,sf]=FFT_SHIFT(t1,st);
t2=-50:dt:50;
gt=sinc(fs*t2);
stt=INSERT0(st,sdt/dt);
xt_t=conv(stt,gt);
figure(1);
subplot(211);
plot(t,xt);title('原始信号');
subplot(212);
stem(t1,st);title('抽样信号');

%--------------------------------------

function[f,sf]=FFT_SHIFT(t,st)

% This function is FFT to calculate a signal's Fourier transform

% Input--  t:sampling time; st:signal data;time length must be greater than 2

% Output-- f:sampling frequence;sf:frequen

% Output is the frequency and the signal spectrum

dt=t(2)-t(1);

T=t(end);

df=1/T;

N=length(t);

f=[-N/2;N/2-1]*df;

sf=fft(st);

sf=T/N*fftshift(sf);

%------------------------------------
function [out]=INSERT0(d,M)
N=length(d);
out=zeros(1,M*N);
for i=0:N-1
     out(i*M+1)=d(i+1);
end;

http://s7/middle/59a069d5naf0dab3505f6&690
% Ex6_3.m

t=[0:0.1:2*pi];
s=sin(t);
[sqnr16,xqtz16,code16]=UniPcm(s,16);
plot(t,s,t+0.2,xqtz16,'*');
title('sin(t)的均匀量化（N=16）')

% -------------------------------------------

function [sqnr,x_qtz,code]=UniPcm(x,n)
xmax=max(abs(x));
x_qtz=x/xmax;
b_qtz=x_qtz;
delta=2/n;
q=delta*[0:n-1]-(n-1)/2*delta;
for i=1:n
    index=find((q(i)-delta/2<=x_qtz)&(x_qtz<=q(i)+delta/2));
    x_qtz(index)=q(i)*ones(1,length(index));
    b_qtz(find(x_qtz==q(i)))=(i-1)*ones(1,length(find(x_qtz==q(i))));
end
x_qtz=x_qtz*xmax;
nu=ceil(log2(n));
code=zeros(length(x),nu);
for i=1:length(x)
    for j=nu:-1:0
        if(fix(b_qtz(i)/(2^j))==1)
            code(i,nu-j)=1;
            b_qtz(i)=b_qtz(i)-2^j;
        end
    end
end
sqnr=20*log10(norm(x)./norm(x-x_qtz));

http://s3/middle/59a069d5naf0eff03a022&690
% 【例题6-5】画出A律曲线。
% demo for A law for quantize, filename:A_law.m
% A=87.6  y=Ax/(1+lnA) (0
%         y=(1+lnAx)/(1+lnA)
clear all; close all;
dx=0.01;
x=0:dx:1;
A=87.6;
for i=1:length(x)
 if abs(x(i))<1/A
  ya(i)=A*x(i)/(1+log(A));
 else
  ya(i)=sign(x(i))*(1+log(A*abs(x(i))))/(1+log(A));
 end
end
figure(1);
plot(x,ya,'k.:');
title('A律13折线压缩曲线');
xlabel('x');
ylabel('y');
grid on;
hold on;
xx=[0,1/128,1/64,1/32,1/16,1/8,1/4,1/2,1];
yy=[0,1/8,2/8,3/8,4/8,5/8,6/8,7/8,1];
plot(xx,yy);
stem(xx,yy);
http://s16/middle/59a069d5naf0da838d54f&690

% 【例题6-10】 已知输入信号为x(t)=sin20πt+0.4sin60πt，增量调制器的采样间隔
% 为1ms，量化阶距△=0.3，单位延迟器的初始值为0.求出前50个采样点时刻上的编码输出
% 序列以及解码样值波形。
Ts=1e-3; % 采样间隔
N=50;    % 采样点数
t=[0:N]*Ts; % 时间序列
x=sin(100*pi*t)+0.4*sin(200*pi*t); % 信号
delta=0.3;  % 量化间隔
code1=dltpcm(x,delta); % 编码
xe=depcm(code1,delta); % 解码
subplot(3,1,1);plot(t,x,'-o');axis([0 N*Ts,-2 2]);hold on;
title('原信号及其离散值')
subplot(3,1,2);stairs(t,code1);axis([0 N*Ts,-2 2]);
title('编码输出二进制序列值')
subplot(3,1,3);stairs(t,xe);hold on; % 解码输出
subplot(3,1,3);plot(t,x); % 原信号
title('解码信号与原信号对比')

%-----------------------------------

function codeout=dltpcm(x,delta)
xlen=length(x);
Dk=0;
for k=1:xlen
    err=x(k)-Dk;
    if(err>=0)
        qout=delta;
        codeout(k)=1;
    else
        qout=-delta;
        codeout(k)=0;
    end
    Dk=Dk+qout;
end

%-------------------------------

function xe=depcm(code,delta)
cdlen=length(code);
Dk=0;
for k=1:cdlen
    if(code(k)>0)
        err=delta;
    else
        err=-delta;
    end
    xe(k)=Dk+err;
    Dk=xe(k);
end

%--------------------------参考书 赵鸿图等. 通信原理MATLAB仿真教程[M].北京：人民邮电出版社，2010.