用Matlab計算BD-PSNR

%%
%  這是Matlab版的BD-PSNR計算程式碼
%  與C版本的計算器avsnr4, 存在0.0x的誤差
%%

clc;clear;close all; % 關閉視窗與清除先前資料

%定義區間,給空陣列 [ ] 會自動取值
range = [442.36 66612];
range = log(range); % for BD-PSNR =>  bitrate

%我的方法
A =  ...
[66612 74.59
29822 46.08
6239   32.65
442    23.72];

% 比較方法A
B =  ...
[68774 56.65
31866 41.97
7083 29.43
416    21.58];

% 比較方法B
B = ...
[71331 44.86
21995 38.70
3593 29.31
148    13.28];

%比較方法C
B = ...
[69671 45.43
23729 40.03
4384 31.59
413    24.78];

%%

Ax = A(:,1);
Ay = A(:,2);
Axl = log(Ax);
Bx = B(:,1);
By = B(:, 2);
Bxl = log(Bx);

figure;plot(Axl,Ay,Bxl,By);title('PSNR-BR');
figure;plot(Ax,Ay,Bx,By);title('PSNR-log(BR)');


% 多項式fitting
rdcA = polyfit(Axl,Ay,3);
rdcB = polyfit(Bxl,By,3);

%自動取值或是在限定範圍內取值
if isempty(range)
    maxX = min(Axl(1), Bxl(1));
    minX = max(Axl(end), Bxl(end));
else
    maxX = range(2);
    minX = range(1);
end

x = minX:0.001:maxX;
figure;plot(x,polyval(rdcA,x), x,polyval(rdcB,x));
title('Estimated curve in PSNR-log(BR)');


% 積分,計算BD

funA = @(x) polyval(rdcA,x);
funB = @(x) polyval(rdcB,x);

sumA = quad(funA,minX,maxX);
sumB = quad(funB,minX,maxX);

BdPsnr = (sumA-sumB)/(maxX-minX)