matlab 曲线拟合--视频编码中PSNR计算及码率计算

matlab 曲线拟合分为多项式拟合和一般曲线拟合

 一、多项式拟合

用到的函数为：

a=polyfit(xdata,ydata,n);

n表示多项式的最高阶数；

（我遇到的问题是要拟合一般曲线，因此多项式拟合带过）；

 二、一般曲线拟合

 [para,resnorm] = lsqcurvefit(fun, x0, xdata, ydata);

 其中para便是我想要的一般曲线中的系数；

比方说我要拟合如下函数的系数，a,b,c,d

y = (a + b*x + c*x²)/(x + d) 

那么para=[a,b,c,d];

 x0: 是给定的一个para的初始值；这里有自己来定，我看有人给的全是1，我也照此做，是可行的。

xdata：就是自变量x的已知系数；

ydata：就是自变量y的已知系数；

 

 举例：

%myfun.m

function SNR=myfun(a,bit)
SNR=(a(1)+a(2)*bit+a(3)*(bit.^2))./(bit+a(4)); 

 %curve.m

 bs_bit=[47304.38,19553.66,8973.04,4603.99];

bs_SNR=[47.661,45.280,43.107,40.882];
lb=[];
ub=[];
[para,res]=lsqcurvefit(@myfun,ones(1,4),bs_bit,bs_SNR,lb,ub,optimset('MaxFunEvals',16000,'MaxIter',800));

 这里我们看到需要为lsqcurvefit定义两个option；

MaxFunEvals和MaxIter

这里我认为算法用来拟合数据所需要的迭代次数；

在上一篇文章中，我们根据测出的四组psnr与bit rate之间的关系，根据这个关系以及VCEG-M34的提议拟合出相对的曲线。

现在需要计算你所采用的方法对psnr和bit rate 的影响。那么就需要求对应相同的psnr时候，码流的变化，和对应相同的码流是对应的psnr的变化。对于后者，我们根据上篇文章拟合出来的公式即可获得。

对于前者，我们需要知道其反函数才能求出对应的值。 

然而，对于符号表达的参数如何求其反函数呢？matalb提供了相应的解决方法。首先先需要下面这三个函数。

syms: 指定符号变量，如，syms x,t,y

sym(finverse(myfun(par0,para1)));

subs( ) 

 下面通过举例来求得反函数计算：

%fmyfun.m
function y=fmyfun(para,snr)
syms bit
f=sym(finverse(myfun(para,bit))) ;
y=subs(f,'bit',snr);

 

这样，就可以实现在已知snr的情况下，求得bit rate

下面给出一个应用例子：

SNR_org = [
  %blue_sky.yuv
  47.658,45.278,43.129,40.932;
  %pedestrian_area.yuv
  47.362,44.574,42.903,41.470;
  %riverbed.yuv
  46.845,43.768,41.320,39.330;
  %rush_hour.yuv
  47.000,44.486,43.223,42.086;
  %station2.yuv
  46.700,43.697,41.978,40.400;
  %sunflower.yuv
  47.156,45.147,43.680,41.922;
  %tractor.yuv
  46.903,43.875,41.375,39.215;
];
bit_org = [
  %blue_sky.yuv
  46552.18,19108.14,8675.51,4435.90;
  %pedestrain_area.yuv
  62321.42,24765.34,9946.34,5424.77;
  %river_bed
  125062.30,75294.25,44999.84,27501.96;
  %rush_hour
  60778.31,21838.84,8316.55,4418.00;
  %station2
  67541.70,20353.78,4219.07,1900.40;
  %sunflower
  40896.95,12166.38,5558.72,2951.06;
  %tractor
  97024.69,49946.76,23004.60,11202.05;
];
bit_org = bit_org/1000;

SNR_25candi_us = [
  %blue_sky
  47.655,45.270,43.122,40.918;
  %pedestrain_area
  47.356,44.569,42.902,41.467;
  %riverbed
  46.847,43.764,41.316,39.323;
  %rush hour
  46.997,44.479,43.212,42.076;
  %station2
  46.690,43.688,41.965,40.388;
  %sunflower
  47.142,45.130,43.670,41.917;
  %tractor
  46.894,43.855,41.355,39.198;
];

bit_25candi_us = [
  %blue_sky
  46785.35,19217.93,8763.00,4475.66
  %pedestrain_area
  62328.03,24823.76,9958.28,5435.45;
  %riverbed
  124661.84,75122.54,44870.01,27449.57;
  %rush_hour
  60824.93,21945.10,8330.94,4427.84;
  %station2
  67662.64,20556.07,4223.99,1892.38;
  %sunflower
  41094.02,12226.13,5552.23,2954.66;
  %tractor
  97280.39,50153.28,23102.80,11236.82;
];
bit_25candi_us = bit_25candi_us/1000;

av_snr = zeros(1,7);
av_bit = zeros(1,7);
for seq = 1:1:7
  bit_T8x8off = bit_25candi_us(seq,:);
  snr_T8x8off = SNR_25candi_us(seq,:);
  [para_T8x8off,res]=lsqcurvefit(@myfun,ones(1,4),bit_T8x8off,snr_T8x8off);

  bit_T8x8on = bit_org(seq,:);
  snr_T8x8on = SNR_org(seq,:);
  [para_T8x8on,res]=lsqcurvefit(@myfun,ones(1,4),bit_T8x8on,snr_T8x8on);

  bit_min = max(bit_T8x8off(1),bit_T8x8on(1));
  bit_max = min(bit_T8x8off(4),bit_T8x8on(4));

  bit     = bit_max:0.1:bit_min;
  snr_off = myfun(para_T8x8off,bit);
  snr_on  = myfun(para_T8x8on, bit);

  snr_min = max(snr_T8x8off(1),snr_T8x8on(1));
  snr_max = min(snr_T8x8off(4),snr_T8x8on(4));
  snr     = snr_max:0.1:snr_min;
  bit_off = fmyfun(para_T8x8off,snr);
  bit_on  = fmyfun(para_T8x8on, snr);

  er_snr  = snr_on - snr_off;
  er_bit  = (bit_on - bit_off)./bit_on;

  av_snr(1,seq) = sum(er_snr)/(size(er_snr,2));
  av_bit(1,seq) = sum(er_bit)/(size(er_bit,2));