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 ));