📄 bindct16v2.m
字号:
function [dctcoef, bindct] = bindct16v2(ver,scaling)%16-point binDCT based on Loeffler's factorization.
% Usage: [dctcoef, bindct] = bindct16v2(ver,scaling, plotall)% Input:% ver: binDCT configuration: 1-4;% 1: 49 shifts, coding gain 9.397 dB;% 2: 51 shifts, coding gain 9.4499 dB% 3: 53 shifts, coding gain 9.45 dB.% 4: 56 shifts, coding gain 9.4515dB.% scaling: 1: include the final scaling. The output bindct is floating-point;% 0: no final scaling. The output bindct is dyadic.% Return:% dctcoef: DCT coefficient;% bindct: binDCT coefficient;%% Notice:% This version uses Loeffler's 8-point DCT in the even channels.% See Ref 1, Fig. 10.
%% Reference:% 1. J. Liang, T. D. Tran, Fast Multiplierless Approximations of the DCT with the Lifting% Scheme, IEEE Trans. Signal Processing, Vol. 49, No. 12, pp. 3032-3044, Dec. 2001. % 2. C. Loeffler, A. Lightenberg, G. Moschytz, Practical fast 1-D DCT algorithms with% 11 multiplications, Prc. IEEE ICASSP-89, Vol. 2, pp. 988-991, Feb. 1999.% % Trac D. Tran and Jie Liang% ECE Department, The Johns Hopkins University% 3400 North Charles Street, 105 Barton Hall,% Baltimore, MD 21218% E-mail: trac@jhu.edu, jieliang@jhu.edu% Dec. 2000% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -% Copyright (c) 2000 Trac D Tran and Jie Liang% This program is Copyright (c) by Trac D Tran and Jie Liang.% It may not be redistributed without the consent of the copyright% holders. In no circumstances may the copyright notice be removed.% The program may not be sold for profit nor may they be incorporated% in commercial programs without the written permission of the copyright% holders. This program is provided as is, without any express or% implied warranty, without even the warranty of fitness for a% particular purpose.%-----------------------------------------------------------------------close all;
if nargin == 0
help bindct16v2; break;
end
if nargin == 1 scaling = 1;end if ver == 1 disp('16-point binDCT-1: 49 shifts, coding gain 9.397 dB;');elseif ver == 2 disp('16-point binDCT-2: 51 shifts, coding gain 9.4499 dB.');elseif ver == 3 disp('16-point binDCT-3: 53 shifts, coding gain 9.45 dB.');elseif ver == 4 disp('16-point binDCT-4: 56 shifts, coding gain 9.4515dB.');endliftcoef = zeros(10,3,10);
%
% Floating DCT: CG: 9.4555 dB, Stop: 10.3668 dB.
%
%
% 49 shifts. CG: 9.397 dB,
liftcoef(:,:,1) = [ 1, 1/2, 0;
-7/16, 3/8, 0; %switch X12 and X4: C12
7/16, -3/8, 0; %switch X13 and X3: C-12
7/16, -3/8, 0; %switch X5 and X11: C-12
-5/16, 9/16, -1/4; %C6: 3pi/16, use result in Fig. 1 instead of Fig. 9.
-3/32, 3/16, -3/32; %C2: pi/16
-11/32, 5/8, -11/32; % C7
-19/32, 7/8, -19/32; % C-11. Remember: reverse direction of P and U !!!
-1/4, 5/16, -1/4; % C3
-1+3/32, 1, -1+3/32; % C-15: Remember: reverse direction of P and U !!!
];
% 51 shifts, CG: 9.4499 dB
liftcoef(:,:,2) = [ 1, 1/2, 0;
-7/16, 3/8, 0; %switch X12 and X4: C12
7/16, -3/8, 0; %switch X13 and X3: C-12
7/16, -3/8, 0; %switch X5 and X11: C-12
-5/16, 9/16, -1/4; %C6: 3pi/16, use result in Fig. 1 instead of Fig. 9.
-3/32, 3/16, -3/32; %C2: pi/16
-11/32, 5/8, -11/32; % C7
-19/32, 7/8, -19/32; % C-11. Remember: reverse direction of P and U !!!
-5/32, 5/16, -5/32; % C3
-1+3/32, 1, -1+3/32; % C-15: Remember: reverse direction of P and U !!!
];
% 53 shifts, coding gain 9.45 dB.
liftcoef(:,:,3) = [ 1, 1/2, 0;
-13/32, 3/8, 0; %switch X12 and X4: C12
13/32, -3/8, 0; %switch X13 and X3: C-12
13/32, -3/8, 0; %switch X5 and X11: C-12
-5/16, 9/16, -1/4; %C6: 3pi/16, use result in Fig. 1 instead of Fig. 9.
-3/32, 3/16, -3/32; %C2: pi/16
-11/32, 5/8, -11/32; % C7
-19/32, 7/8, -19/32; % C-11. Remember: reverse direction of P and U !!!
-5/32, 5/16, -5/32; % C3
-1+3/32, 1, -1+3/32; % C-15: Remember: reverse direction of P and U !!!
];
% 56 shifts, coding gain 9.4515dB.
liftcoef(:,:,4) = [ 1, 1/2, 0;
-13/32, 11/32, 0; %switch X12 and X4: C12
13/32, -11/32, 0; %switch X13 and X3: C-12
13/32, -11/32, 0; %switch X5 and X11: C-12
-5/16, 9/16, -1/4; %C6: 3pi/16, use result in Fig. 1 instead of Fig. 9.
-3/32, 3/16, -3/32; %C2: pi/16
-11/32, 5/8, -11/32; % C7
-19/32, 7/8, -19/32; % C-11. Remember: reverse direction of P and U !!!
-5/32, 5/16, -5/32; % C3
-1+3/32, 1, -1+3/32; % C-15: Remember: reverse direction of P and U !!!
];
M = 16;
%define the row and column indexesridx = [ 0 : M - 1]' ; % column vectorcidx = ridx' ; % row vector%calculate the DCT coefficientsdctcoef = dct(eye(M));afrplot(dctcoef,inv(dctcoef)');
set(gcf, 'name', 'DCT Frequency Response');bplot(dctcoef, dctcoef);set(gcf, 'name', 'DCT Basis Functions');%
% Calculate the matrix of the new BinDCT
%
stage1 = eye(16);
stage2 = eye(16);
stage3 = eye(16);
stage4 = eye(16);
stage5 = eye(16);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 1:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i8 = eye(8);
j8 = fliplr(i8);
stage1 = [i8, j8;
j8, -i8];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 2:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i4 = eye(4);
j4 = fliplr(i4);
butfly8 = [i4, j4;
j4, -i4];
stage2(1:8,1:8) = butfly8;
%-------
% X8 and X15: angle is C7
%-------
r7_p1 = eye(16);
r7_p1(16, 9) = liftcoef(7,1,ver);
r7_u = eye(16);
r7_u(9, 16) = liftcoef(7,2,ver);
r7_p2 = eye(16);
r7_p2(16, 9) = liftcoef(7,3,ver);
r7 = r7_p2 * r7_u * r7_p1;
%-------
% 9 and 14: C-11
%-------
r11_p1 = eye(16);
r11_p1(10, 15) = liftcoef(8,1,ver);
r11_u = eye(16);
r11_u(15, 10) = liftcoef(8,2,ver);
r11_p2 = eye(16);
r11_p2(10, 15) = liftcoef(8,3,ver);
r11 = r11_p2 * r11_u * r11_p1;
%-------
% 10 and 13: C3
%-------
r3_p1 = eye(16);
r3_p1(14, 11) = liftcoef(9,1,ver);
r3_u = eye(16);
r3_u(11, 14) = liftcoef(9,2,ver);
r3_p2 = eye(16);
r3_p2(14, 11) = liftcoef(9,3,ver);
r3 = r3_p2 * r3_u * r3_p1;
%-------
% 11 and 12: C-15
%-------
r15_p1 = eye(16);
r15_p1(12, 13) = liftcoef(10,1,ver);
r15_u = eye(16);
r15_u(13, 12) = liftcoef(10,2,ver);
r15_p2 = eye(16);
r15_p2(12, 13) = liftcoef(10,3,ver);
r15 = r15_p2 * r15_u * r15_p1;
%----------------
%stage2_low: cmbination of the four angles above.
%----------------
foo = r7 * r11 * r3 * r15;
stage2(9:16, 9:16) = foo(9:16, 9:16);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 3:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
i2 = eye(2);
j2 = fliplr(i2);
butfly4 = [i2, j2;
j2, -i2];
stage3(1:4,1:4) = butfly4;
%-------
% C6: This is from the 8-point DCT in Fig. 1 of Loeffler's paper.
%in 8-point, the angle is denoted C3, in 16-point, it should be C6.
%-------
r6_p1 = eye(4);
r6_p1(4, 1) = liftcoef(5,1,ver);
r6_u = eye(4);
r6_u(1, 4) = liftcoef(5,2,ver);
r6_p2 = eye(4);
r6_p2(4, 1) = liftcoef(5,3,ver);
r6 = r6_p2 * r6_u * r6_p1;
%-------
% C2: This is from the 8-point DCT in Fig. 1 of Loeffler's paper.
%in 8-point, the angle is denoted C1, in 16-point, it should be C2.
%-------
r2_p1 = eye(4);
r2_p1(3, 2) = liftcoef(6,1,ver);
r2_u = eye(4);
r2_u(2, 3) = liftcoef(6,2,ver);
r2_p2 = eye(4);
r2_p2(3, 2) = liftcoef(6,3,ver);
r2 = r2_p2 * r2_u * r2_p1;
stage3(5:8, 5:8) = r2 * r6;
%------
% stage 3 lower part
%------
i2 = eye(2);
j2 = fliplr(i2);
butfly4 = [i2, j2;
j2, -i2];
stage3(9:12, 9:12) = butfly4;
stage3(13:16, 13:16) = butfly4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 4:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-------
% X0, X8: pi/4
%-------
predictor = eye(2);
predictor(1, 2) = liftcoef(1,1,ver); %1
updator = eye(2);
updator(2, 1) = liftcoef(1,2,ver); %1/2;
updator(2, 2) = -1;
stage4(1:2, 1:2) = updator * predictor; %without considering the scaling factors
%-------
% X12, X4: 3pi/8
%-------
predictor = eye(2);
predictor(1, 2) = liftcoef(2,1,ver); %-3/8; %-7/16;
updator = eye(2);
updator(2, 1) = liftcoef(2,2,ver); %3/8; %11/32;
stage4(3:4, 3:4) = updator * predictor;
stage4(5:8, 5:8) = [ 1 0 1 0;
0 -1 0 1;
1 0 -1 0;
0 1 0 1];
%-------------------------
% stage4 lower part
%-------------------------
stage4(12,12) = -1;
stage4(14,14) = -1;
stage4(15,15) = -1;
stage4(16,16) = -1;
stage4(9,15) = 1;
stage4(15,9) = 1;
stage4(10,16) = 1;
stage4(16,10) = 1;
stage4(11,12) = 1;
stage4(12,11) = 1;
stage4(13,14) = 1;
stage4(14,13) = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 5:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Y14 and Y2: use a lifting to ensure 3 butterflies for each channel.
stage5(5:8, 5:8) = [-1/2 0 0 1/2;
0 1 0 0;
0 0 1 0;
1 0 0 1];
%-------
% Y13, Y3: C-12
%-------
predictor = eye(2);
predictor(1, 2) = liftcoef(3,1,ver); %7/16;
updator = eye(2);
updator(2, 1) = liftcoef(3,2,ver); %-3/8;
stage5(9:10, 9:10) = updator * predictor;
% Y15 and Y1: use a lifting to ensure 3 bf for each channel.
stage5(12:13, 12:13) = [ 1, 1; 1/2, -1/2];
predictor = eye(2);
predictor(1, 2) = liftcoef(4,1,ver); %7/16;
updator = eye(2);
updator(2, 1) = liftcoef(4,2,ver); %-3/8;
stage5(15:16, 15:16) = updator * predictor;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Get the matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bindct = stage5 * stage4 * stage3 * stage2 * stage1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%final scaling
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if scaling == 1
bindct(1,:) = bindct(1,:) / 4;
bindct(2,:) = bindct(2,:) / 2;
bindct(3,:) = bindct(3,:) * (-sin(3*pi/8)) * sqrt(2) / 4; %sqrt(2) from where?
bindct(4,:) = bindct(4,:) / sin(3*pi/8) * sqrt(2) / 4; %sqrt(2) from where?
bindct(5,:) = bindct(5,:) / 2; % lifting!!!
bindct(6,:) = bindct(6,:) * sqrt(2) / 4; %sqrt(2) from Fig. 1
bindct(7,:) = bindct(7,:) * sqrt(2) / 4; %sqrt(2) from Fig. 1
bindct(8,:) = bindct(8,:) / 4;
bindct(9,:) = bindct(9,:) * sin(3*pi/8) * sqrt(2) / 4; %sqrt(2) from Fig. 9
bindct(10,:) = bindct(10,:) / (-sin(3*pi/8)) * sqrt(2) / 4; %sqrt(2) from Fig. 9
bindct(11,:) = bindct(11,:) * sqrt(2) / 4;
bindct(12,:) = bindct(12,:) / 4;
bindct(13,:) = bindct(13,:) * (-1) / 2; %-x[1] is from Fig. 9. 2 from lifting.
bindct(14,:) = bindct(14,:) * sqrt(2) / 4;
bindct(15,:) = bindct(15,:) * sin(3*pi/8) * sqrt(2) / 4; %sqrt(2) from Fig. 1
bindct(16,:) = bindct(16,:) / (-sin(3*pi/8)) * sqrt(2) * (-1) / 4; %sqrt(2) from Fig. 1. -1 from Fig. 9. -X11.
end
%rearrange the rows
tmp = bindct;
bindct = zeros(16,16);
bindct(1,:) = tmp(1,:);
bindct(9,:) = tmp(2,:);
bindct(13,:) = tmp(3,:);
bindct(5,:) = tmp(4,:);
bindct(15,:) = tmp(5,:);
bindct(7,:) = tmp(6,:);
bindct(11,:) = tmp(7,:);
bindct(3,:) = tmp(8,:);
bindct(14,:) = tmp(9,:);
bindct(4,:) = tmp(10,:);
bindct(10,:) = tmp(11,:);
bindct(16,:) = tmp(12,:);
bindct(2,:) = tmp(13,:);
bindct(8,:) = tmp(14,:);
bindct(6,:) = tmp(15,:);
bindct(12,:) = tmp(16,:);
afrplot(bindct,inv(bindct)');
set(gcf, 'name', 'binDCT Frequency Response');bplot(bindct, bindct);set(gcf, 'name', 'binDCT Basis Functions');if scaling == 1
error = dctcoef - bindct;
Rxx=toeplitz(0.95.^[0:M-1]');
mse = trace( error * Rxx * error') / M
end
break;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of file%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -