function [pXinDao,D,RD] = myRD(Pu,S,d,tol)
%率失真函数迭代算法
%输入信源概率分布Pu，失真函数d，容忍度tol，拉氏乘子S
%输出信道矩阵，失真度D，在D下的率失真RD
%测试环境：MATLAB2018b

A = exp(S .* d); 
col = size(d,2);
Pv = ones(1,col)/col; %初始输出概率分布Pv

Tu = 0;
TL = 1;
while abs(Tu - TL) > tol        
    c = (Pu * A) ./ (Pv * A');
    Pv = Pv .* c;
    Tu = sum(Pv .* c .* log2(c));
    TL = max(log2(c));
end

pXinDao = (Pv .* A) ./ sum(Pv .* A,2)';
D = sum(sum(Pu .* pXinDao .* d));
RD = S*D - Pu * log2(sum(A.*Pv,2)) - sum(Pv .* c .* log2(c));

%myRDTest.m
%调用率失真函数
%测试环境：MATLAB2018b

format long
%Pu = [0.5 0.5]; %信源分布
%d = [0 1;1 0]; %失真函数
%tol = 1e-6; %容忍度
%meanD = 0.4; %失真限度

Pu = [0.25 0.25 0.25 0.25];
d = [0 1 1 1;1 0 1 1; 1 1 0 1; 1 1 1 0];
tol = 1e-6;
meanD = 0.2;

D = zeros(1,100);
RD = zeros(1,100);
pXinDao = zeros(size(d,1),size(d,2),10);

times = 1;
for S = -30:0.1:0
    [pXinDao(:,:,times),D(times),RD(times)] = myRD(Pu,S,d,tol);
    times = times+1;
end

RDmax = RD(1);
RDmin = RD(times-1);
Dmax = D(times-1);
Dmin = D(1);

plot(D,RD);
xlabel('失真度D');
ylabel('率失真函数');
title('信息率失真函数曲线图');
grid on

hold on
[~,index] = sort(abs(D-meanD));
plot(D(index(1)),RD(index(1)),'ro','LineWidth',2); %标记在失真限度处的值

xinDao = pXinDao(:,:,index(1));
%dim = [.4 .5 .6 .3];
%str = [num2str(pXinDao(:,:,index(1)))];
%annotation('textbox',dim,'String',str,'FitBoxToText','on'); %指明对应信道矩阵