📄 dfilters.m
字号:
% reconstruction filter
h0 = sqrt(2)*ftrans2(h,h_trans);
h1 = modulate2(h0,'b');
h0 = sqrt(2)*ftrans2(f,h_trans);
end
case 'optim9' % ------------------------------------------------------
% filter design by optimization method, synthesis and analysis
% fiters are the same centro-symmetric filters
% high pass filter is modulation of lowpass filter
h0 = 10^(-3)*[...
0.0744 -0.3361 -2.5823 0.7787 -2.0561 0.7966 -0.6866 -0.6434 0.0665
0.3058 4.7015 -7.7192 -5.4991 -2.1719 -10.6472 -3.4500 1.6892 0.6543
-2.0812 6.8068 12.5743 -34.8545 -27.5217 -43.0428 26.0018 3.6422 -0.9481
-0.5519 -2.7023 28.8237 74.4245 -278.9884 64.1302 45.6976 -10.8442 -0.9527
-2.3530 3.4648 -23.8808 266.2680 816.0931 266.2680 -23.8808 3.4648 -2.3530
-0.9527 -10.8442 45.6976 64.1302 -278.9884 74.4245 28.8237 -2.7023 -0.5519
-0.9481 3.6422 26.0018 -43.0428 -27.5217 -34.8545 12.5743 6.8068 -2.0812
0.6543 1.6892 -3.4500 -10.6472 -2.1719 -5.4991 -7.7192 4.7015 0.3058
0.0665 -0.6434 -0.6866 0.7966 -2.0561 0.7787 -2.5823 -0.3361 0.0744];
h1 = h0;
h1(2:2:end, 1:2:end) = -h1(2:2:end, 1:2:end);
h1(1:2:end, 2:2:end) = -h1(1:2:end, 2:2:end);
case 'optim11'
h0 = 10^(-3)*sqrt(2)*[ ...
-0.0457 0.4225 -0.1463 0.1600 1.1171 -0.4570 1.1171 0.1600 -0.1463 0.4225 -0.0457;
0.3557 0.0362 0.1650 -2.7032 -0.5488 -4.4739 -0.5488 -2.7032 0.1650 0.0362 0.3557;
0.0845 -0.1945 -7.3984 8.3702 5.9456 8.9636 5.9456 8.3702 -7.3984 -0.1945 0.0845;
0.1172 -2.7413 8.7215 21.4156 -43.1093 -7.2819 -43.1093 21.4156 8.7215 -2.7413 0.1172;
0.8169 -0.2165 5.4719 -40.6566 -31.7533 184.6749 -31.7533 -40.6566 5.4719 -0.2165 0.8169;
-0.2627 -3.6227 6.7719 -7.6965 183.1327 585.7587 183.1327 -7.6965 6.7719 -3.6227 -0.2627;
0.8169 -0.2165 5.4719 -40.6566 -31.7533 184.6749 -31.7533 -40.6566 5.4719 -0.2165 0.8169;
0.1172 -2.7413 8.7215 21.4156 -43.1093 -7.2819 -43.1093 21.4156 8.7215 -2.7413 0.1172;
0.0845 -0.1945 -7.3984 8.3702 5.9456 8.9636 5.9456 8.3702 -7.3984 -0.1945 0.0845;
0.3557 0.0362 0.1650 -2.7032 -0.5488 -4.4739 -0.5488 -2.7032 0.1650 0.0362 0.3557;
-0.0457 0.4225 -0.1463 0.1600 1.1171 -0.4570 1.1171 0.1600 -0.1463 0.4225 -0.0457];
h1 = 10^(-3)*sqrt(2)*[ ...
-0.0406 -0.4278 -0.1415 -0.1653 1.1221 0.4515 1.1221 -0.1653 -0.1415 -0.4278 -0.0406;
-0.3609 0.0412 -0.1705 -2.6982 0.5436 -4.4691 0.5436 -2.6982 -0.1705 0.0412 -0.3609;
0.0895 0.1893 -7.3936 -8.3754 5.9507 -8.9690 5.9507 -8.3754 -7.3936 0.1893 0.0895;
-0.1225 -2.7362 -8.7270 21.4207 43.1040 -7.2771 43.1040 21.4207 -8.7270 -2.7362 -0.1225;
0.8220 0.2113 5.4767 40.6513 -31.7482 -184.6804 -31.7482 40.6513 5.4767 0.2113 0.8220;
0.2575 -3.6176 -6.7774 -7.6915 -183.1379 585.7635 -183.1379 -7.6915 -6.7774 -3.6176 0.2575;
0.8220 0.2113 5.4767 40.6513 -31.7482 -184.6804 -31.7482 40.6513 5.4767 0.2113 0.8220;
-0.1225 -2.7362 -8.7270 21.4207 43.1040 -7.2771 43.1040 21.4207 -8.7270 -2.7362 -0.1225;
0.0895 0.1893 -7.3936 -8.3754 5.9507 -8.9690 5.9507 -8.3754 -7.3936 0.1893 0.0895;
-0.3609 0.0412 -0.1705 -2.6982 0.5436 -4.4691 0.5436 -2.6982 -0.1705 0.0412 -0.3609;
-0.0406 -0.4278 -0.1415 -0.1653 1.1221 0.4515 1.1221 -0.1653 -0.1415 -0.4278 -0.0406];
case 'meyer' % -------------------------------------------------------
% estimate meyer type diamond filter bank
N = 128;
[x1, x2] = meshgrid(0:2*pi/(N):2*pi,0:2*pi/(N):2*pi);
X = x1 + x2;
F1 = zeros(size(X));
F1( find( (X < 2*pi/3) ) )= 1;
int2 = find( (X >= 2*pi/3) & (X < 4*pi/3) );
F1( int2 ) = cos(pi/2*meyeraux(3/2/pi*X(int2)-1));
Fs = F1.^2+rot90(F1.^2,1)+rot90(F1.^2,2)+rot90(F1.^2,3);
Fs = Fs(1:N,1:N);
% F2 = fftshift(F1);
% Gs = F2.^2+rot90(F2.^2,1)+rot90(F2.^2,2)+rot90(F2.^2,3);
Gs = fftshift(Fs);
F = sqrt(Fs./(Fs+Gs));
G = sqrt(Gs./(Fs+Gs));
h0 = ifftshift(ifft2(F));
h1 = ifftshift(ifft2(G));
L = 30;
cen = N/2 + 1;
h0 = sqrt(2)*h0(cen-L:cen+L,cen-L:cen+L);
h1 = sqrt(2)*h1(cen-L:cen+L,cen-L:cen+L);
% nomralized , improve 1 to 2 db
% h0 = modulate2(h0,'c'); h0 = h0./sum(h0(:)); h0 = modulate2(h0,'c');
% h1 = modulate2(h1,'c'); h1 = h1./sum(h1(:)); h1 =
% modulate2(h1,'c');
case 'meyerh2' % estimate meyer type fan fb for the second level of the dual dfb
N = 128;
[x1, x2] = meshgrid(0:2*pi/(N):2*pi,0:2*pi/(N):2*pi);
X = x1 + x2;
F1 = zeros(size(X));
F1( find( (X < 2*pi/3) ) )= 1;
int2 = find( (X >= 2*pi/3) & (X < 4*pi/3) );
F1( int2 ) = cos(pi/2*meyeraux(3/2/pi*X(int2)-1));
Fs = F1.^2+rot90(F1.^2,1)+rot90(F1.^2,2)+rot90(F1.^2,3);
Fs = Fs(1:N,1:N);
Gs = fftshift(Fs);
%calculate the FG mask for the dual case, second level
alpha = 0.2*pi;
X2 = -abs(x1-x2)*(2*pi/(3*alpha))+ 4*pi/3;
X1 = fftshift(X2); X2(find(abs(x1-x2)> pi) ) = X1(find(abs(x1-x2)> pi) );
Fdm = zeros(size(X2)); % Fdual mask
Fdm( find( X2 < 2*pi/3) )= 1;
int2 = find( X2 >= 2*pi/3 );
Fdm( int2 ) = cos(pi/2*meyeraux(3/2/pi*X2(int2)-1));
F = sqrt(2*Fs./(Fs+Gs));
G = sqrt(2*Gs./(Fs+Gs));
F = F([N/2+1:end,1:N/2],:);
G = G([N/2+1:end,1:N/2],:);
% F and G are the FFT2 of Fan FB
% prepare the phase
[w1, w2] = meshgrid(-pi:2*pi/(N):pi, -pi:2*pi/(N):pi);
phase = ((w1+w2 )>= 0).*(-w1/2-w2/2+pi/2) + ...
((w1+w2 )< 0).*(-w1/2-w2/2-pi/2) ;
% make the phase border periodic
temp = phase(1,:)+phase(end,:);phase(1,:) = temp/2;phase(end,:) = temp/2;
temp = phase(:,1)+phase(:,end);phase(:,1) = temp/2;phase(:,end) = temp/2;
% make the phase symmetric
phase = phase - rot90(phase,2);
% cut the phase to NN
phase = 0.5*phase(1:N,1:N);
% not necessary, since the phase is the same is shifted by pi pi
phase = fftshift(phase);
% dual mask
Fdm = rot90(Fdm); Fdm = Fdm(1:N,1:N);
F = F.*exp(-i*Fdm.*phase);
G = G.*exp(-i*Fdm.*phase);
if lower(type(1)) == 'r'
F = conj(F); G = conj(G);
end
h1 = ifftshift(ifft2(F));
h0 = ifftshift(ifft2(G));
L = 30;
cen = N/2 + 1;
% h0 = real(fitmat(h0,2*L));
% h1 = real(fitmat(h1,2*L));
h0 = real(h0(cen-L:cen+L,cen-L:cen+L));
h1 = real(h1(cen-L:cen+L,cen-L:cen+L));
case 'meyerh3' % estimate meyer type fan fb for the second level of the dual dfb
N = 128;
[x1, x2] = meshgrid(0:2*pi/(N):2*pi,0:2*pi/(N):2*pi);
X = x1 + x2;
F1 = zeros(size(X));
F1( find( (X < 2*pi/3) ) )= 1;
int2 = find( (X >= 2*pi/3) & (X < 4*pi/3) );
F1( int2 ) = cos(pi/2*meyeraux(3/2/pi*X(int2)-1));
Fs = F1.^2+rot90(F1.^2,1)+rot90(F1.^2,2)+rot90(F1.^2,3);
Fs = Fs(1:N,1:N);
Gs = fftshift(Fs);
%calculate the FG mask for the dual case, second level
alpha = 0.2*pi;
X2 = -abs(x1-x2)*(2*pi/(3*alpha))+ 4*pi/3;
X1 = fftshift(X2); X2(find(abs(x1-x2)> pi) ) = X1(find(abs(x1-x2)> pi) );
Fdm = zeros(size(X2)); % Fdual mask
Fdm( find( X2 < 2*pi/3) )= 1;
int2 = find( X2 >= 2*pi/3 );
Fdm( int2 ) = cos(pi/2*meyeraux(3/2/pi*X2(int2)-1));
F = sqrt(2*Fs./(Fs+Gs));
G = sqrt(2*Gs./(Fs+Gs));
F = F([N/2+1:end,1:N/2],:);
G = G([N/2+1:end,1:N/2],:);
% F and G are the FFT2 of Fan FB
% prepare the phase
[w1, w2] = meshgrid(-pi:2*pi/(N):pi, -pi:2*pi/(N):pi);
phase = ((-w1+w2 )>= 0).*(w1/2-w2/2+pi/2) + ...
((-w1+w2 )< 0).*(w1/2-w2/2-pi/2) ;
% make the phase border periodic
temp = phase(1,:)+phase(end,:);phase(1,:) = temp/2;phase(end,:) = temp/2;
temp = phase(:,1)+phase(:,end);phase(:,1) = temp/2;phase(:,end) = temp/2;
% make the phase symmetric
phase = phase - rot90(phase,2);
% cut the phase to NN
phase = 0.5*phase(1:N,1:N);
% not necessary, since the phase is the same is shifted by pi pi
phase = fftshift(phase);
Fdm = Fdm(1:N,1:N);
F = F.*exp(-i*Fdm.*phase);
G = G.*exp(-i*Fdm.*phase);
if lower(type(1)) == 'r'
F = conj(F); G = conj(G);
end
h1 = ifftshift(ifft2(F));
h0 = ifftshift(ifft2(G));
L = 30;
cen = N/2 + 1;
% h0 = real(fitmat(h0,2*L));
% h1 = real(fitmat(h1,2*L));
h0 = real(h0(cen-L:cen+L,cen-L:cen+L));
h1 = real(h1(cen-L:cen+L,cen-L:cen+L));
otherwise
error('Unrecognized directional filter name');
end
%-----------------------------------------------------------------------
% Internal function
% ----------------------------------------------------------------------
function h_trans=trans_matrix(N)
% this function return the transformation matrix as in Tay and Kingsburry paper
%
% properties : size 2N+1*2N+1
% h_trans(k1,k2)= 0 when k1+k2 even
% h_trans(k1,k2) linear phase
% h_trans(w1,w2)=1 in the passband (diamond)
% and =-1 in the stopband (outside the dimamond)
[n1,n2]=meshgrid(-N:N,-N:N);
h_ideal=sinc((n1+n2)/2).*sinc((n1-n2)/2);
h_ideal(N+1,N+1)=0;
% choose Kaiser window
beta=2.5;
w = kaiser(2*N+1,beta);
% index matrix
n3=n1+n2;
n4=n1-n2;
% set up window
w1=zeros(2*N+1,2*N+1);
w2=zeros(2*N+1,2*N+1);
for i=-N:1:N
% Kaiser window for n1+n2
w1(abs(n3)==i)=w(i+N+1);
% Kaiser window for n1-n2
w2(abs(n4)==i)=w(i+N+1);
end
% This is the two dimension diamond Kaiser window
win=w1.*w2;
h_trans=h_ideal.*win;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -