03/14/0018:47:241/u/fxiao/WWW/ee392J/GetInterTable.m%%% Quantization matrix for inter quantization matrix%%% which is build by raising the intra matrix to a frame-rate %%% dependent powerfunction table= GetInterTable(frame_rate)power_ratio= 1+frame_rate/30;s_ratio=4;Intra_table= [ 8 16 19 22 26 27 29 34; 16 16 22 24 27 29 34 37; 19 22 26 27 29 34 34 38; 22 22 26 27 29 34 37 40; 22 26 27 29 32 35 40 48; 26 27 29 32 35 40 48 58; 26 27 29 34 38 46 56 69; 27 29 35 38 46 56 69 83];table= s_ratio*power( Intra_table/s_ratio, power_ratio);03/14/0018:47:241/u/fxiao/WWW/ee392J/GetIntraTable.m%%% MPEG default intra quantization matrixfunction table= GetIntraTabletable= [ 8 16 19 22 26 27 29 34; 16 16 22 24 27 29 34 37; 19 22 26 27 29 34 34 38; 22 22 26 27 29 34 37 40; 22 26 27 29 32 35 40 48; 26 27 29 32 35 40 48 58; 26 27 29 34 38 46 56 69; 27 29 35 38 46 56 69 83];%table= table/8;%table= power(table,1.7)*8;%table= round(table);03/14/0018:47:081/u/fxiao/WWW/ee392J/QP_effect.m%% show how SNR, MSE, DVQ varies with quantization factorclear all;close all;rows = 144;cols = 176;image_size=rows*cols;frame_size=image_size*3/2;xblocks= floor(cols/8);yblocks= floor(rows/8);orig_file= ’car.qcif’;fid = fopen(orig_file,’r’);[temp,count] = fread(fid,[cols,rows],’uchar’);fclose(fid);lum=temp’;block= zeros(8,8);block_dct=block;QP= 1:1:8;QP_levels= length(QP);dvq_mean=zeros(1,QP_levels);dvq_max=zeros(1,QP_levels);dvq_best=zeros(1,QP_levels);mse=zeros(1,QP_levels);snr=zeros(1,QP_levels);for k=1:QP_levels recon_lum= zeros(size(lum)); for y=1:yblocks; for x=1:xblocks; block=lum(y*8-7:y*8,x*8-7:x*8); block_dct=dct2(block); q_dct= block_dct/QP(k)/2; q_dct(q_dct>0)=q_dct(q_dct>0)+0.5; q_dct(q_dct<0)=q_dct(q_dct<0)-0.5; q_dct=round(q_dct); q_dct= q_dct*2*QP(k); if y==4 & x== 4 q_dct(1,1) = q_dct(1,1)*0.8; end q_dct= idct2(q_dct); recon_lum(y*8-7:y*8,x*8-7:x*8)=q_dct; end end figure; imshow(recon_lum,[0 255]); [dvq_mean(k),dvq_max(k),dvq_best(k)]=compute_dvq(lum,recon_lum); %mse(k)=compute_mse(lum,recon_lum); %snr(k)=compute_snr(lum,recon_lum);enddvq_meandvq_maxmsefigure; clf; plot(QP,dvq_mean); title(’DVQ Mean’);03/14/0018:47:082/u/fxiao/WWW/ee392J/QP_effect.mfigure; clf; plot(QP,dvq_max); title(’DVQ Max’);figure; clf; plot(QP,dvq_best); title(’DVQ Best’);figure; clf; plot(QP,mse); title(’MSE’);%figure; clf; plot(QP,snr); title(’SNR’);03/14/0018:47:091/u/fxiao/WWW/ee392J/SCSF.m%% Human spatial-temporal contrast sensitivity functionclose all;clear all;[fs,ft]= meshgrid(linspace(1,0,20),linspace(0,1,20));ratio_t=2; ratio_s= 3; ratio_st= 1;K= (fs-0.3).^2*ratio_s + ((ft-0.2).*(fs+0.5)).^2*ratio_t;K= exp(-K);K= exp(K); rat= max(max(log10(K)));exp_ft= exp(ft*log(10));exp_fs= exp(fs*log(10));surf(log10(exp_ft),log10(exp_fs),100*log10(K)/rat);xlabel(’Temporal Frequency (log unit)’);ylabel(’Spatial Frequency ( log unit)’);zlabel(’Contrast Sensitivity’);title(’Figure 1: Human Spatial-Temporal Contrast Sensitivity Function’);return;03/14/0018:48:031/u/fxiao/WWW/ee392J/block_filter.mfunction filtered_blocks= block_filter(filename,block_y,block_x,total_frames)rows = 144;cols = 176;total_frames=256;blocks= read_block(filename,block_y,block_x,total_frames);lc =zeros(total_frames, 8, 8);%%%% compute the local contrast of each block%%%% 1. DCT of each block. %%%% 2. DCT coeffeciences divided by its DC coeffeciences. for frameNo=1:total_frames block=squeeze(blocks(frameNo,:,:)); dct_block=dct2(block); ratio= power(dct_block(1,1)/1024,0.65); dc=(dct_block(1,1)+0.1)/ratio; lc(frameNo,:,:)=dct_block/dc;end%%%% second-order II temperal filterII_filter=[ linspace(0,1,10) linspace(1,0,10) ];a=[1];filtered_blocks= filter(II_filter,a,lc);03/14/0018:48:041/u/fxiao/WWW/ee392J/block_filter_hold.mfunction [filtered_blocks1,filtered_blocks2,mse]= block_filter_hold(filename,block_y,block_x,total_frames,drop_frames)rows = 144;cols = 176;total_frames=256;blocks= read_block(filename,block_y,block_x,total_frames);blocks_drop= frame_drop(blocks,drop_frames);mse= compute_mse(blocks,blocks_drop);lc =zeros(total_frames, 8, 8);%%%% compute the local contrast of each block%%%% 1. DCT of each block. %%%% 2. DCT coeffeciences divided by its DC coeffeciences. for frameNo=1:total_frames block=squeeze(blocks(frameNo,:,:)); dct_block=dct2(block); ratio= power(dct_block(1,1)/1024,0.65); dc=(dct_block(1,1)+0.1)/ratio; lc(frameNo,:,:)=dct_block/dc;end%%%% second-order II temperal filterII_filter=[linspace(0,1,10) linspace(1,0,10) ];a=[1];filtered_blocks1= filter(II_filter,a,lc);lc_drop=frame_drop(lc,drop_frames);filtered_blocks2= filter(II_filter,a,lc_drop);03/14/0018:47:521/u/fxiao/WWW/ee392J/compute_mse.mfunction mse= compute_mse(im1,im2)im_diff=im1-im2;im_diff=im_diff.^2;dims=ndims(im_diff);count=1;for dim=1:dims count=count*size(im_diff,dim);endfor dim=1:dims im_diff = sum(im_diff);endmse=im_diff/prod(size(im1));mse=sqrt(mse);03/14/0018:47:521/u/fxiao/WWW/ee392J/compute_snr.mfunction snr= compute_snr(im1,im2)mse= compute_mse(im1,im2);sq_mse= mse*mse;snr= 10*log10(255*255/sq_mse);03/14/0018:47:521/u/fxiao/WWW/ee392J/compute_vqm.mfunction [vqm, vqm_mean, vqm_max] = compute_vqm(im1, im2)[rows, cols]=size(im1);Intra_table= GetIntraTable;yblocks= floor(rows/8); xblocks= floor(cols/8);block= zeros(8,8);vqm_matrix=zeros(rows,cols);%im_diff=im1-im2;for y=1:yblocks for x=1:xblocks block=im1(y*8-7:y*8,x*8-7:x*8); block_dct=dct2(block); ratio= power(block_dct(1,1)/1024,0.65); dc=(block_dct(1,1)+0.1)/ratio; block_dct=block_dct/dc; %*ratio; %block_dct/dc; vqm1= 8*block_dct./Intra_table; block=im2(y*8-7:y*8,x*8-7:x*8); block_dct=dct2(block); ratio= power(block_dct(1,1)/1024,0.65); dc=(block_dct(1,1)+0.1)/ratio; block_dct=block_dct/dc; vqm2= 8*block_dct./Intra_table; vqm_matrix(y*8-7:y*8,x*8-7:x*8)= vqm2-vqm1; endendvqm_matrix= 1000*abs(vqm_matrix);vqm_mean= mean(mean(vqm_matrix));vqm_max= max(max(vqm_matrix));vqm= (vqm_mean+0.005*vqm_max)/1.01;03/14/0018:47:541/u/fxiao/WWW/ee392J/drop.m%% VQM distortion for dropping different amount of frames function vqm=drop(video_file,rows,cols)image_size=rows*cols;frame_size=image_size*3/2;xblocks= floor(cols/8);yblocks= floor(rows/8);fid =