📄 bindct.m
字号:
function [dctcoef, bindct, bindct4] = bindct(ver, scaling)% This file shows the coefficients and frequency responses of
% the 7 configuration of binDCT Type C, as shown in the
% SPIE paper, Config. 8 is the reduced dynamic range version of C1.
%
% Jie Liang, Feb 07, 2000
% Last modified date: 11/23/2000
%
%
%Set this flag to 0 to disable the scaling of binDCT coef, used to
%check the dynamic range of the binDCt coefs.
if nargin ~= 2
disp('Freq Response of various versions of binDCT-C:');
disp('Usage: [dctcoef, bindct, bindct4] = bindct(version, scaling):');
disp(' Version: 1-->8');
disp(' scaling: 1 to scaling the binDCT; 0: not scaling;'); disp(' dctcoef: DCT coefficients;');
disp(' bindct: binDCT coefficients.');
disp(' bindct4: 4-point binDCT.'); break;
end
M = 8;
%define the row and column indexes
ridx = [ 0 : M - 1]' ; % column vector
cidx = ridx' ; % row vector
%calculate the DCT coefficients
dctcoef = sqrt(2 / M) * cos(ridx * pi * (2 * cidx + 1) / (2 * M) );
dctcoef(1, : ) = dctcoef(1, : ) / sqrt(2);
close all;
liftcoef = zeros(5,3,20);
%binDCT-C1
% 9 shifts, 28 Adds.
%CG: 8.7686dB, Stopband: 8.2823dB, MSE: 2.3E-3.
%The fwd and inv coefs are in the format of x/128 and x/512.
%If we set u3 to 0, this will reduce to x/32, and x/128, see Config. 8
liftcoef(:,:,1) = [1, 1/2, 0; % pi/4 between X0 and X4.
-1/2, 1/2, 0; % p1 and u1 in the paper.
-1/2, 3/4, 1/2; % p4, u4, p5
-1/4, 1/4, 0; % p3 u3
1, -1/2 0; % p2 and u2
];
%binDCT-C2
% 14 shifts, 33 Adds.
%CG: 8.8033dB, Stopband: 8.9536dB, MSE: 5.7799E-4.
liftcoef(:,:,2) = [1, 1/2, 0;
% -1/2, 3/8, 0;
-3/8, 1/4, 0; -7/16, 3/4, 3/8;
-3/16, 1/4, 0;
7/8, -1/2 0;
];
%binDCT-C3
% 17 shifts, 36 Adds.
%CG: 8.8159dB, Stopband: 8.1849dB, MSE: 4.1928E-4.
liftcoef(:,:,3) = [1, 1/2, 0;
-3/8, 3/8, 0;
-7/16, 11/16, 3/8;
-3/16, 3/16, 0;
7/8, -1/2 0;
];
%binDCT-C4:
%---> 19 shifts, 37 Adds.
%CG: 8.822dB, stopband: 10.0214dB, MSE: 8.4897E-5
liftcoef(:,:,4) = [1, 1/2, 0;
-7/16, 3/8, 0; %7/16 replaces 3/8, no change in complexity.
-7/16, 11/16, 3/8;
-3/16, 3/16, 0;
5/8, -7/16 0; %5/8 replaces 7/8, 7/16 replaces 1/2
%but 15/32 replace 7/16 will be worse.!!!
];
%binDCT-C5
%21 shifts, 40 Adds.
%CG: 8.8233dB, stopband: 9.5246dB, MSE: 3.3725E-5
liftcoef(:,:,5) = [1, 1/2, 0;
-13/32, 11/32, 0; %13/32 rep 3/8: +1S, 1A; 11/32 rep 3/8: +1S, 1A.
-7/16, 11/16, 3/8;
-3/16, 3/16, 0;
11/16, -15/32, 0; %11/16 rep 7/8: +1S(1-1/4-1/16 vs 1-1/8),1A; 15/32 rep 1/2: +1S, 1A.
];
%binDCT-C6:
%21 shifts, 39 Adds.
%CG: 8.8240dB, stopband: 10.0923dB, MSE: 5.6968E-5
liftcoef(:,:,6) = [1, 1/2, 0;
-7/16, 3/8, 0;
-13/32, 11/16, 13/32; % + 2S, 2A than 19-shift version.
-3/16, 3/16, 0;
5/8, -7/16, 0;
];
%binDCT-C7
%23 shifts, 42 Adds.
%CG: 8.8251dB, stopband: 9.5813dB, MSE: 1.1069E-5
liftcoef(:,:,7) = [1, 1/2, 0;
-13/32, 11/32, 0;
-13/32, 11/16, 13/32;
-3/16, 3/16, 0;
11/16, -15/32, 0;
];
%binDCT-C0
% set u3(or p3?) to 0, this has CG of 8.7483dB.
% The dynamic range reduces to 32 and 128.
% 8 shifts, 27 Adds.
%CG: 8.7483dB, Stopband: 7.8452dB, MSE: 2.27E-3.
liftcoef(:,:,8) = [1, 1/2, 0; % pi/4 between X0 and X4.
-1/2, 1/2, 0; % p1 and u1 in the paper.
-1/2, 3/4, 1/2; % p4, u4, p5
% -1/2, 1, 1/2; % denominator: 8 and 32. 8.63dB.
-1/4, 0, 0;
1, -1/2 0; % p2 and u2
];
%lift (:,:,9): binDCT-A: but
%Don't switch X[2] and X[6],X[1] and X[7].
%In this case, the freq resps of the 6th, 7th rows are not good
%
%MSE: 0.00116.
liftcoef(:,:,9) = [1, 1/2, 0;
2.40625, -3/8, 0;
% -3/8, 3/8, 0;
% 2.375, -3/8, 0;
-7/16, 11/16, 3/8;
5+1/64+1/128, -3/16, 0;
% 5+1/32, -3/16, 0;
% -3/16, 3/16, 0;
7/8, -1/2 0;
];
%binDCT-C10 modification:
%CG: 8.8257dB, stopband: 9.9883dB, MSE: 1.038E-6
liftcoef(:,:,10) = [1, 1/2, 0;
-53/128, 45/128, 0;
-53/128, 181/256, 53/128;
-25/128, 49/256, 0;
85/128, -59/128, 0;
];
%8shifts
%8.6886dB
%gcd of fwd: 8, gcd of inv: 32. Use this in paper.
liftcoef(:,:,11) = [1, 1/2, 0; % pi/4 between X0 and X4.
-1/2, 1/2, 0; % p1 and u1 in the paper.
-1/2, 1, 1/2; % p4, u4, p5
-1/4, 0, 0; % p3 u3
1, -1/2 0; % p2 and u2
];
%5 shifts
%8.4083dB
%gcd of fwd: 8, gcd of inv: 32.
liftcoef(:,:,12) = [1, 1/2, 0; % pi/4 between X0 and X4.
-1, 1/2, 0; % p1 and u1 in the paper.
0, 1/2, 1/2; % p4, u4, p5
0, 0, 0; % p3 u3
-1, 1/2, 0; % p2 and u2
];
%How to get WHT?%liftcoef(:,:,13) = [1, 1/2, 0; % pi/4 between X0 and X4.
% -1, 1/2, 0; % p1 and u1 in the paper.
% 0, 0, 0; % p4, u4, p5
% -1, 1/2, 0; % p3 u3
% 1, -1/2, 0; % p2 and u2
%];
%Show the frequency response
freq = linspace(0, pi, 256); % Equally spaced points between [0, 2pi].
allplot = figure;
title('DCT Frequency Response');
for i = 1 : M
h = abs( freqz(dctcoef(i, :), 1, freq) );
% h = h / max(h);
figure(allplot);
hold on;
plot(freq, 20*log10(h));
fid(i) = figure;
plot(freq, 20*log10(abs(h)),'-.');
v=axis;
v(1) = 0;
v(2) = pi;
v(3) = -50;
v(4) = 10;
axis(v);
s = sprintf('Frequency Response of DCT Basis Function No: %d', i);
set(gcf, 'Name', s);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');
end
figure(allplot);
v=axis;
v(1) = 0;
v(2) = pi;
v(3) = -50;
v(4) = 10;
axis(v);
s = sprintf('Frequency Response of All DCT Basis Functions');
set(gcf, 'Name', s);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');
%
% Calculate the matrix of the new BinDCT
%
%First stage: 8-row butterfly
i4 = eye(4);
j4 = fliplr(i4);
stage1 = [i4, j4;
j4, -i4];
%4-row butterfly
i2 = eye(2);
j2 = fliplr(i2);
butfly4 = [i2, j2;
j2, -i2];
stage2 = eye(8);
stage2(1:4,1:4) = butfly4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upper branch rotation of pi/4: up 1; down 1/2, min x4
% for X[0] and X[4]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
predictor = eye(8);
predictor(1, 2) = liftcoef(1,1,ver); %1
updator = eye(8);
updator(2, 1) = liftcoef(1,2,ver); %1/2;
updator(2, 2) = -1;
rot_pi4 = updator * predictor; %without considering the scaling factors
%verify rotation pi/4
factor = eye(8);
factor(1,1) = sqrt(2)/2;
factor(2,2) = sqrt(2);
true_rot_pi4 = factor * updator * predictor; % with scaling factors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%X[6] and X[2]:
%Upper branch rotation of 3pi/8
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
predictor = eye(8);
%predictor(3, 4) = 1 - sqrt(2);
predictor(3, 4) = liftcoef(2,1,ver); %-3/8; %-7/16;
updator = eye(8);
%updator(4, 3) = sqrt(2) / 4;
updator(4, 3) = liftcoef(2,2,ver); %3/8; %11/32;
rot_pi38 = updator * predictor;
%verify rotation 3pi/8
factor = eye(8);
factor(3,3) = -sin(3*pi/8);
factor(4,4) = 1 / sin(3*pi/8);
true_rot_pi38 = factor * rot_pi38; % with scaling factors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%First rotation in matrix V: lower branch
%between the butterfly:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vlift_p1 = eye(8);
Vlift_p1(6, 7) = liftcoef(3,1,ver); %-7/16; %-13/32;
Vlift_u = eye(8);
Vlift_u(7, 6) = liftcoef(3,2,ver); %11/16; %23/32;
Vlift_p2 = eye(8);
Vlift_p2(6, 6) = -1;
Vlift_p2(6, 7) = liftcoef(3,3,ver); %3/8; %13/32;
low_rot_pi4 = Vlift_p2 * Vlift_u * Vlift_p1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Two 2-row butterfly in V
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tmp1 = [1 1;
1 -1];
tmp2 = [-1 1;
1 1];
butfly2 = eye(8);
butfly2(5:6, 5:6) = tmp1;
butfly2(7:8, 7:8) = tmp2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rotation 7pi/16: X{7] and X[1]
%switch X[1] and X[7], as which leads bad freq response of the 8th rows
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lift716_p = eye(8);
%lift716_p(5, 8) = - cos(7*pi/16) / sin(7*pi/16);
lift716_p(5, 8) = liftcoef(4,1,ver); %-3/16;
lift716_u = eye(8);
%lift716_u(8, 5) = sin(7*pi/16) * cos(7*pi/16);
lift716_u(8, 5) = liftcoef(4,2,ver); %3/16;
rot_pi716 = lift716_u * lift716_p;
%verify rotation 7pi/16
factor = eye(8);
factor(5,5) = -sin(7*pi/16);
factor(8,8) = 1 / sin(7*pi/16);
true_rot_pi716 = factor * rot_pi716; % with scaling factors
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rotation 3pi/16
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lift316_p = eye(8);
lift316_p(6, 7) = liftcoef(5,1,ver); %7/8; %better choice: 5/8;
lift316_u = eye(8);
lift316_u(7, 6) = liftcoef(5,2,ver); %-1/2; %-7/16;
rot_pi316 = lift316_u * lift316_p;
%verify rotation 3pi/16
factor = eye(8);
factor(6,6) = cos(3*pi/16);
factor(7,7) = 1 / cos(3*pi/16);
true_rot_pi316 = factor * rot_pi316; % with scaling factors
%
%multiply and get the overall matrix,
%without consifering the final scaling parameters.
%
bindct = rot_pi316 * rot_pi716 * butfly2 * low_rot_pi4 * rot_pi38 * rot_pi4 * stage2 * stage1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4-point binDCT%%%%%%%%%%%%%%%%%%%%%%%%%%%%bindct4 = rot_pi38 * rot_pi4 * stage2;bindct4 = bindct4(1:4,1:4);%rearrange:tmp = bindct4;bindct4(1,:) = tmp(1,:);bindct4(2,:) = tmp(4,:);bindct4(3,:) = tmp(2,:);bindct4(4,:) = -tmp(3,:); %The minus sign comes from the implementation of this program. If use the new binDCT structure, it's absorbed in the lifting structure.if scaling == 1 %normalization bindct4 = diag(1./sqrt(sum((bindct4.^2)'))) * bindct4;endbplot(dct(eye(4)), bindct4);afrplot(bindct4, inv(bindct4)');
%Get the inverse transform matrix
if ver == 1
%pi/4
predictor = eye(8);
predictor(2, 1) = 1/2;
predictor(2,2) = -1;
updator = eye(8);
updator(1, 2) = -1;
inv_pi4 = updator * predictor;
%3pi/8
predictor = eye(8);
predictor(4, 3) = -3/8;
updator = eye(8);
updator(3, 4) = 3/8;
inv_pi38 = updator * predictor;
stage2 = eye(8,8);
stage2(1,1)=1; stage2(1,4)= 1;
stage2(4,1)= 1; stage2(4,4)= -1;
stage2(2,2) =1; stage2(2,3) = 1;
stage2(3,2) =1; stage2(3,3) = -1;
% 3pi/16
predictor = eye(8);
predictor(7, 6) = 1/2;
updator = eye(8);
updator(6, 7) = -7/8;
inv_pi316 = updator * predictor;
% 7pi/16
predictor = eye(8);
predictor(8, 5) = -3/16;
updator = eye(8);
updator(5, 8) = 3/16;
inv_pi716 = updator * predictor;
stage2_low = eye(8);
stage2_low(5,5)=1; stage2_low(5,6) =1;
stage2_low(6,5) =1; stage2_low(6,6) =-1;
stage2_low(7,7)=-1; stage2_low(7,8) =1;
stage2_low(8,7) =1; stage2_low(8,8) =1;
%pi/4
predictor = eye(8);
predictor(6, 6) = -1;
predictor(6,7) =3/8;
updator = eye(8);
updator(7, 6) = -11/16;
updator(7, 7) = 1;
predictor2 = eye(8);
predictor2(6, 6) = 1;
predictor2(6,7) =7/16;
inv_pi4_2 = predictor2 * updator * predictor;
inv_bin = stage1 * stage2 * inv_pi4_2 * stage2_low * inv_pi716 * inv_pi316 * inv_pi4 * inv_pi38;
inv_bin * bindct
end
%final scaling
if scaling == 1
bindct(1,:) = bindct(1,:) * sqrt(2) / 2;
bindct(2,:) = bindct(2,:) * (sqrt(2));
if ver == 9
%don't switch X[2] and X[6], X[1] and X[7].
% to show that the sensitivity of the result to angles close to pi/2.
bindct(3,:) = bindct(3,:) * cos(3*pi/8);
bindct(4,:) = bindct(4,:) / cos(3*pi/8);
bindct(5,:) = bindct(5,:) * cos(7*pi/16);
bindct(8,:) = bindct(8,:) / cos(7*pi/16);
else
bindct(3,:) = bindct(3,:) * (-1) *sin(3*pi/8);
bindct(4,:) = bindct(4,:) / sin(3*pi/8);
bindct(5,:) = bindct(5,:) * (-sin(7*pi/16));
bindct(8,:) = bindct(8,:) / sin(7*pi/16);
end
bindct(6,:) = bindct(6,:) * cos(3*pi/16);
bindct(7,:) = bindct(7,:) / cos(3*pi/16);
bindct = bindct / 2; %inheritated from factorization.
end
%rearrange the rows
tmp = bindct;
bindct = zeros(8,8);
bindct(1,:) = tmp(1,:);
bindct(5,:) = tmp(2,:);
bindct(6,:) = tmp(6,:);
bindct(4,:) = tmp(7,:);
if ver == 9
bindct(3,:) = tmp(3,:);
bindct(7,:) = tmp(4,:);
bindct(2,:) = tmp(5,:);
bindct(8,:) = tmp(8,:);
else
bindct(7,:) = tmp(3,:);
bindct(3,:) = tmp(4,:);
bindct(8,:) = tmp(5,:);
bindct(2,:) = tmp(8,:);
end
bplot(dctcoef, bindct);
%Show the frequency response
freq = linspace(0, pi, 256); % Equally spaced points between [0, 2pi].
allplot = figure;
title('binDCT Frequency Response');
for i = 1 : M
h = abs( freqz(bindct(i, :), 1, freq) );
% h = h / max(h);
figure(allplot);
hold on;
plot(freq, 20*log10(h));
figure(fid(i));
hold on;
plot(freq, 20*log10(abs(h)));
v=axis;
v(1) = 0;
v(2) = pi;
v(3) = -50;
v(4) = 10;
axis(v);
s = sprintf('Frequency Response of DCT Basis Function No: %d', i);
set(gcf, 'Name', s);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');
legend('DCT','binDCT');
end
figure(allplot);
v=axis;
v(1) = 0;
v(2) = pi;
v(3) = -50;
v(4) = 10;
axis(v);
s = sprintf('Frequency Response of All DCT Basis Functions');
set(gcf, 'Name', s);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');
afrplot(bindct,inv(bindct)');
disp('error between DCT and binDCT:');
error = dctcoef - bindct;
disp('MSE between DCT and binDCT:');
mse = sum(sum(error .* error)) / M / M;
Rxx=toeplitz(0.95.^[0:7]');
mse = trace( error * Rxx * error') / M
mse = trace( error * error') / M;
break;
figure;
for r = 1 : 2
for c = 1 : 4
subplot(2,4,(r-1)*4+c);
h = abs( freqz(dctcoef((r-1)*4+c, :), 1, freq) );
plot(freq, 20*log10(h));
hold on;
h = abs( freqz(bindct((r-1)*4+c, :), 1, freq) );
plot(freq, 20*log10(abs(h)),'-.');
v=axis;
v(1) = 0;
v(2) = pi;
v(3) = -50;
v(4) = 10;
axis(v);
end
end
subplot(2,4,5);
ylabel('Magnitude Response(dB)');
xlabel('Frequency(rad/sec)');
subplot(2,4,8);
legend('DCT','binDCT C4');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% End of file
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -