BM3D/CVBM3D.m

   1 function [Xdenoised] = CVBM3D(Xnoisy, sigma, Xorig)
   2 %  CVBM3D denoising of RGB videos corrupted with AWGN.
   3 %
   4 %
   5 %  [Xdenoised] = CVBM3D(Xnoisy, sigma, Xorig)
   6 %
   7 %  INPUTS:
   8 %
   9 %   1)  Xnoisy --> Either a filename of a noisy AVI RGB uncompressed video (e.g. 'SMg20.avi')
  10 %                  or a 4-D matrix of dimensions (M x N x 3 x NumberOfFrames)
  11 %                  The intensity range is [0,255]!
  12 %   2)  Sigma -->  Noise standard deviation (assumed intensity range is [0,255])
  13 %
  14 %   3)  Xorig     (optional parameter) --> Filename of the original video
  15 %
  16 %  OUTPUT: .avi files are written to the current matlab folder
  17 %
  18 %   1) Xdenoised --> A 4-D matrix with the denoised RGB-video
  19 %
  20 %  USAGE EXAMPLES:
  21 %   1) To denoise a video:
  22 %      CVBM3D('SMg20.avi', 20)
  23 %
  24 %   2) To denoise a video and print PSNR:
  25 %      CVBM3D('SMg20.avi', 20, 'SM.avi')
  26 %
  27 %   1) To denoise a 4-D matrix representing a noisy RGB video:
  28 %      CVBM3D(X_4D_matrix, 20)
  29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  30 %
  31 % Copyright © 2009 Tampere University of Technology. All rights reserved.
  32 % This work should only be used for nonprofit purposes.
  33 %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35
  36 % If no input argument is provided, then use the internal ones from below:
  37 if exist('sigma', 'var') ~= 1,
  38     Xnoisy = 'SMg20.avi';  sigma = 20;  ;
  39 end
  40
  41 % Whether or not to print information to the screen
  42 dump_information = 1;
  43
  44 % If the input is a 4-D matrix, then save it as AVI file that is used as
  45 % input to the denoising
  46 if ischar(Xnoisy) == 0;
  47     NumberOfFrames = size(Xnoisy,4);
  48
  49     if NumberOfFrames <= 1
  50         error('The input RGB video should be a 4-D matrix (M x N x 3 x NumberOfFrames)');
  51     end
  52     avi_filename = sprintf('ExternalMatrix_%.6d.avi', round(rand*50000));
  53     if exist(avi_filename, 'file') == 2,
  54         delete(avi_filename);
  55     end
  56     mov = avifile(avi_filename, 'Colormap', gray(256), 'compression', 'None', 'fps', 30);
  57     if mean2(Xnoisy) <= 1
  58         fprintf('Possible error: the input RGB-videos should be in range [0,255] and not in [0,1]!\n');
  59     else
  60         for ii = [1:NumberOfFrames],
  61             mov = addframe(mov, uint8(Xnoisy(:,:,:,ii)));
  62         end
  63     end
  64     mov = close(mov);
  65
  66     if dump_information == 1
  67         fprintf('The input 4-D matrix was written to: %s.\n', avi_filename);
  68     end
  69
  70     clear Xnoisy
  71     Xnoisy = avi_filename;
  72 end
  73
  74 % Read some properties of the noisy RGB video
  75 noi_avi_file_info = aviinfo(Xnoisy);
  76 NumberOfFrames = noi_avi_file_info.NumFrames;
  77
  78 %%% Read Xorig video --- needed if one wants to compute PSNR and ISNR
  79 if exist('Xorig', 'var') == 1,
  80     if ischar(Xorig) == 1;
  81         org_avi_file_info = aviinfo(Xorig);
  82         mo = aviread(Xorig);
  83         Xorig = zeros([size(mo(1).cdata), NumberOfFrames], 'single');
  84         for cf = 1:NumberOfFrames
  85             Xorig(:,:,:,cf) = single(mo(cf).cdata(:,:,:));
  86         end
  87         clear mo;
  88
  89         if (org_avi_file_info.NumFrames == noi_avi_file_info.NumFrames && org_avi_file_info.FramesPerSecond == noi_avi_file_info.FramesPerSecond && ...
  90                 org_avi_file_info.Width == noi_avi_file_info.Width && org_avi_file_info.Height == noi_avi_file_info.Height)
  91             dump_information = 1;
  92         end
  93     else
  94         Xorig = single(Xorig);
  95         if mean2(Xorig) <= 1
  96             fprintf('Possible error: the input RGB-videos should be in range [0,255] and not in [0,1]!\n');
  97         end
  98
  99     end
 100 end
 101
 102 denoiseFrames  = min(9, NumberOfFrames);
 103 denoiseFramesW = min(9, NumberOfFrames);
 104
 105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 106 %%%%  Quality/complexity trade-off
 107 %%%%
 108 %%%%  'np' --> Normal Profile (balanced quality)
 109 %%%%  'lc' --> Low Complexity Profile (fast, lower quality)
 110 %%%%
 111 if (exist('bm3dProfile') ~= 1)
 112     bm3dProfile         = 'np';
 113 end
 114
 115 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 116 %%%% Following are the parameters for the Normal Profile.
 117 %%%%
 118
 119 %%%% Select transforms ('dct', 'dst', 'hadamard', or anything that is listed by 'help wfilters'):
 120 transform_2D_HT_name     = 'bior1.5'; %% transform used for the HT filt. of size N1 x N1
 121 transform_2D_Wiener_name = 'dct';     %% transform used for the Wiener filt. of size N1_wiener x N1_wiener
 122 transform_3rd_dim_name   = 'haar'; %% tranform used in the 3-rd dim, the same for HT and Wiener filt.
 123
 124 %%%% Step 1: Hard-thresholding (HT) parameters:
 125 N1                  = 8;  %% N1 x N1 is the block size used for the hard-thresholding (HT) filtering
 126 Nstep               = 5;  %% sliding step to process every next refernece block
 127 N2                  = 8;  %% maximum number of similar blocks (maximum size of the 3rd dimension of the 3D groups)
 128 Ns                  = 7;  %% length of the side of the search neighborhood for full-search block-matching (BM)
 129 Npr                 = 3;  %% length of the side of the motion-adaptive search neighborhood, use din the predictive-search BM
 130 tau_match           = 3000; %% threshold for the block distance (d-distance)
 131 lambda_thr3D        = 2.7; %% threshold parameter for the hard-thresholding in 3D DFT domain
 132 dsub                = 13;  %% a small value subtracted from the distnce of blocks with the same spatial coordinate as the reference one
 133 Nb                  = 2;  %% number of blocks to follow in each next frame, used in the predictive-search BM
 134 beta                = 2.0; %% the beta parameter of the 2D Kaiser window used in the reconstruction
 135
 136
 137 %%%% Step 2: Wiener filtering parameters:
 138 N1_wiener           = 7;
 139 Nstep_wiener        = 4;
 140 N2_wiener           = 8;
 141 Ns_wiener           = 7;
 142 Npr_wiener          = 3;
 143 tau_match_wiener    = 1000;
 144 beta_wiener         = 2.0;
 145 dsub_wiener         = 1.5;
 146 Nb_wiener           = 2;
 147
 148 %%%% Block-matching parameters:
 149 stepFS              = 1; %% step that firces to switch to full-search BM, "1" implies always full-search
 150 stepFSW             = 1;
 151 thrToIncStep        = 8;  %% used in the HT filtering to increase the sliding step in uniform regions
 152
 153
 154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 155 %%%% Following are the parameters for the Low Complexity Profile.
 156 %%%%
 157 if strcmp(bm3dProfile, 'lc') == 1,
 158     lambda_thr3D = 2.8;
 159     denoiseFrames  = min(5, NumberOfFrames);
 160     denoiseFramesW = min(5, NumberOfFrames);
 161     N2_wiener = 4;
 162     N2 = 4;
 163     Ns = 3;
 164     Ns_wiener = 3;
 165     Nb = 1;
 166     Nb_wiener = 1;
 167 end
 168
 169 if strcmp(bm3dProfile, 'hi') == 1,
 170     Nstep        = 3;
 171     Nstep_wiener = 3;
 172 end
 173
 174 if sigma > 30,
 175     N1_wiener = 8;
 176     tau_match    = 4500;
 177     tau_match_wiener    = 3000;
 178 end
 179
 180 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 182 %%%% Note: touch below this point only if you know what you are doing!
 183 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 184 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 185
 186
 187 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 188 %%%% Create transform matrices, etc.
 189 %%%%
 190 decLevel                     = 0;    %% dec. levels of the dyadic wavelet 2D transform for blocks (0 means full decomposition, higher values decrease the dec. number)
 191 decLevel3                    = 0;    %% dec. level for the wavelet transform in the 3rd dimension
 192
 193 [Tfor, Tinv]   = getTransfMatrix(N1, transform_2D_HT_name, decLevel); %% get (normalized) forward and inverse transform matrices
 194 [TforW, TinvW] = getTransfMatrix(N1_wiener, transform_2D_Wiener_name); %% get (normalized) forward and inverse transform matrices
 195
 196 if (strcmp(transform_3rd_dim_name, 'haar') == 1 || strcmp(transform_3rd_dim_name(end-2:end), '1.1') == 1),
 197     %%% Fast internal transform is used, no need to generate transform
 198     %%% matrices.
 199     hadper_trans_single_den         = {};
 200     inverse_hadper_trans_single_den = {};
 201 else
 202     %%% Create transform matrices. The transforms are later computed by
 203     %%% matrix multiplication with them
 204     for hh = [1 2 4 8 16 32];
 205         [Tfor3rd, Tinv3rd]   = getTransfMatrix(hh, transform_3rd_dim_name, decLevel3);
 206         hadper_trans_single_den{hh}         = single(Tfor3rd);
 207         inverse_hadper_trans_single_den{hh} = single(Tinv3rd');
 208     end
 209 end
 210
 211
 212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 213 %%%% 2D Kaiser windows that scale the reconstructed blocks
 214 %%%%
 215 if beta_wiener==2 & beta==2 & N1_wiener==7 & N1==8 % hardcode the window function so that the signal processing toolbox is not needed by default
 216     Wwin2D = [ 0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924;
 217         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 218         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 219         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 220         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 221         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 222         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 223         0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924 ];
 224     Wwin2D_wiener = [ 0.1924    0.3151    0.4055    0.4387    0.4055    0.3151    0.1924;
 225         0.3151    0.5161    0.6640    0.7184    0.6640    0.5161    0.3151;
 226         0.4055    0.6640    0.8544    0.9243    0.8544    0.6640    0.4055;
 227         0.4387    0.7184    0.9243    1.0000    0.9243    0.7184    0.4387;
 228         0.4055    0.6640    0.8544    0.9243    0.8544    0.6640    0.4055;
 229         0.3151    0.5161    0.6640    0.7184    0.6640    0.5161    0.3151;
 230         0.1924    0.3151    0.4055    0.4387    0.4055    0.3151    0.1924 ];
 231 else
 232     Wwin2D           = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the aggregation of the HT part
 233     Wwin2D_wiener    = kaiser(N1_wiener, beta_wiener) * kaiser(N1_wiener, beta_wiener)'; % Kaiser window used in the aggregation of the Wiener filt. part
 234 end
 235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 236 %%%% Read an image, generate noise and add it to the image
 237 %%%%
 238
 239 if dump_information == 1
 240     fprintf('Input video: %s, sigma: %.1f\n', Xnoisy, sigma);
 241 end
 242
 243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 244 %%%% Determine unique filenames of intermediate avi files
 245 %%%%
 246
 247 HT_avi_file = sprintf('%s_cvbm3d_step1_0.avi', Xnoisy(1:end-4));
 248 Denoised_avi_file = sprintf('%s_cvbm3d_0.avi', Xnoisy(1:end-4));
 249 i = 1;
 250 while (exist(['./' HT_avi_file], 'file') ~= 0) | (exist(['./' Denoised_avi_file],'file') ~= 0)
 251     HT_avi_file = sprintf('%s_cvbm3d_step1_%d.avi', Xnoisy(1:end-4),i);
 252     Denoised_avi_file = sprintf('%s_cvbm3d_%d.avi', Xnoisy(1:end-4),i);
 253     i = i + 1;
 254 end
 255
 256 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 257 %%%% Initial estimate by hard-thresholding filtering
 258 HT_IO = {which(Xnoisy), HT_avi_file};
 259
 260 tic;
 261 bm3d_thr_video_c(HT_IO, hadper_trans_single_den, Nstep, N1, N2, 0,...
 262     lambda_thr3D, tau_match*N1*N1/(255*255), (Ns-1)/2, sigma/255, thrToIncStep,...
 263     single(Tfor), single(Tinv)', inverse_hadper_trans_single_den, single(ones(N1)),...
 264     'unused arg', dsub*dsub/255 * (sigma^2 / 255), ones(NumberOfFrames,1), Wwin2D,...
 265     (Npr-1)/2, stepFS, denoiseFrames, Nb, 0 );
 266 estimate_elapsed_time = toc;
 267
 268 if dump_information == 1
 269 %     mo = aviread(HT_avi_file);
 270 %     y_hat = zeros([size(mo(1).cdata(:,:,1)), 3, NumberOfFrames], 'single');
 271 %     for cf = 1:NumberOfFrames
 272 %         y_hat(:,:,:,cf) = single(mo(cf).cdata(:,:,:))/255;
 273 %     end
 274 %     clear  mo
 275 %
 276 %     PSNR_HT_ESTIMATE = 10*log10(1/mean2((Xorig-y_hat).^2));
 277 %     fprintf('HT ESTIMATE, PSNR: %.3f dB\n', PSNR_HT_ESTIMATE);
 278 %     clear y_hat;
 279      fprintf('STEP1 completed!\n');
 280 end
 281
 282 %%
 283 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 284 % %%%% Final estimate by Wiener filtering (using the hard-thresholding
 285 % initial estimate)
 286
 287 lut_ic = ClipComp16b(sigma/255);
 288
 289 WIE_IO = {which(Xnoisy), HT_avi_file, Denoised_avi_file};
 290
 291 tic;
 292 bm3d_wiener_video_c(WIE_IO, 'unused', hadper_trans_single_den, Nstep_wiener, N1_wiener, N2_wiener, ...
 293     'unused_arg', tau_match_wiener*N1_wiener*N1_wiener/(255*255), (Ns_wiener-1)/2, sigma/255, 'unused arg',...
 294     single(TforW), single(TinvW)', inverse_hadper_trans_single_den, 'unused arg', dsub_wiener*dsub_wiener/255*(sigma^2 / 255),...
 295     ones(NumberOfFrames,1), Wwin2D_wiener, (Npr_wiener-1)/2, stepFSW, denoiseFramesW, Nb_wiener, 0, lut_ic);
 296
 297 wiener_elapsed_time = toc;
 298
 299 if nargout == 1
 300     mo = aviread(Denoised_avi_file);
 301     Xdenoised = zeros([size(mo(1).cdata(:,:,1)), 3, NumberOfFrames], 'single');
 302     for cf = 1:NumberOfFrames
 303         Xdenoised(:,:,:,cf) = single(mo(cf).cdata(:,:,:));
 304     end
 305     clear  mo
 306 end
 307
 308 if dump_information == 1
 309     if nargout ~= 1
 310         mo = aviread(Denoised_avi_file);
 311         Xdenoised = zeros([size(mo(1).cdata(:,:,1)), 3, NumberOfFrames], 'single');
 312         for cf = 1:NumberOfFrames
 313             Xdenoised(:,:,:,cf) = single(mo(cf).cdata(:,:,:));
 314         end
 315         clear  mo
 316     end
 317
 318     PSNR_TEXT='';
 319     if exist('Xorig', 'var') == 1
 320         PSNR = 10*log10(255*255/mean((Xorig(:)-Xdenoised(:)).^2));
 321         PSNR_TEXT=sprintf(' PSNR: %.3f dB,', PSNR);
 322         New_Denoised_avi_file = sprintf('%s_PSNR%.2f.avi',Denoised_avi_file(1:end-4),PSNR);
 323         movefile(Denoised_avi_file, New_Denoised_avi_file);
 324         Denoised_avi_file = New_Denoised_avi_file;
 325     end
 326
 327 %     PSNRs = zeros(NumberOfFrames,1);
 328 %     for ii = 1:NumberOfFrames,
 329 %         PSNRs(ii) = 10*log10(1/mean2( (Xorig(:,:,:,ii)-Xdenoised(:,:,:,ii)).^2));
 330 %         fprintf('Frame: %d, PSNR: %.2f\n', ii, PSNRs(ii));
 331 %     end
 332     if nargout == 0
 333         clear Xdenoised
 334     end
 335
 336     fprintf('FILTERING COMPLETED (frames/sec: %.2f,%s denoised video saved as %s)\n', ...
 337         NumberOfFrames/(wiener_elapsed_time + estimate_elapsed_time), PSNR_TEXT, Denoised_avi_file);
 338
 339 end
 340
 341
 342 return;
 343
 344
 345 function [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 346 %
 347 % Create forward and inverse transform matrices, which allow for perfect
 348 % reconstruction. The forward transform matrix is normalized so that the
 349 % l2-norm of each basis element is 1.
 350 %
 351 % [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 352 %
 353 %  INPUTS:
 354 %
 355 %   N               --> Size of the transform (for wavelets, must be 2^K)
 356 %
 357 %   transform_type  --> 'dct', 'dst', 'hadamard', or anything that is
 358 %                       listed by 'help wfilters' (bi-orthogonal wavelets)
 359 %                       'DCrand' -- an orthonormal transform with a DC and all
 360 %                       the other basis elements of random nature
 361 %
 362 %   dec_levels      --> If a wavelet transform is generated, this is the
 363 %                       desired decomposition level. Must be in the
 364 %                       range [0, log2(N)-1], where "0" implies
 365 %                       full decomposition.
 366 %
 367 %  OUTPUTS:
 368 %
 369 %   Tforward        --> (N x N) Forward transform matrix
 370 %
 371 %   Tinverse        --> (N x N) Inverse transform matrix
 372 %
 373
 374 if exist('dec_levels') ~= 1,
 375     dec_levels = 0;
 376 end
 377
 378 if N == 1,
 379     Tforward = 1;
 380 elseif strcmp(transform_type, 'hadamard') == 1,
 381     Tforward    = hadamard(N);
 382 elseif (N == 8) & strcmp(transform_type, 'bior1.5')==1 % hardcoded transform so that the wavelet toolbox is not needed to generate it
 383     Tforward =  [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 384        0.219417649252501   0.449283757993216   0.449283757993216   0.219417649252501  -0.219417649252501  -0.449283757993216  -0.449283757993216  -0.219417649252501;
 385        0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846  -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284;
 386       -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284   0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846;
 387        0.707106781186547  -0.707106781186547                   0                   0                   0                   0                   0                   0;
 388                        0                   0   0.707106781186547  -0.707106781186547                   0                   0                   0                   0;
 389                        0                   0                   0                   0   0.707106781186547  -0.707106781186547                   0                   0;
 390                        0                   0                   0                   0                   0                   0   0.707106781186547  -0.707106781186547];
 391 elseif (N == 8) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
 392     Tforward = [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 393        0.490392640201615   0.415734806151273   0.277785116509801   0.097545161008064  -0.097545161008064  -0.277785116509801  -0.415734806151273  -0.490392640201615;
 394        0.461939766255643   0.191341716182545  -0.191341716182545  -0.461939766255643  -0.461939766255643  -0.191341716182545   0.191341716182545   0.461939766255643;
 395        0.415734806151273  -0.097545161008064  -0.490392640201615  -0.277785116509801   0.277785116509801   0.490392640201615   0.097545161008064  -0.415734806151273;
 396        0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274   0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274;
 397        0.277785116509801  -0.490392640201615   0.097545161008064   0.415734806151273  -0.415734806151273  -0.097545161008064   0.490392640201615  -0.277785116509801;
 398        0.191341716182545  -0.461939766255643   0.461939766255643  -0.191341716182545  -0.191341716182545   0.461939766255643  -0.461939766255643   0.191341716182545;
 399        0.097545161008064  -0.277785116509801   0.415734806151273  -0.490392640201615   0.490392640201615  -0.415734806151273   0.277785116509801  -0.097545161008064];
 400 elseif (N == 8) & strcmp(transform_type, 'dst')==1 % hardcoded transform so that the PDE toolbox is not needed to generate it
 401     Tforward = [ 0.161229841765317   0.303012985114696   0.408248290463863   0.464242826880013   0.464242826880013   0.408248290463863   0.303012985114696   0.161229841765317;
 402        0.303012985114696   0.464242826880013   0.408248290463863   0.161229841765317  -0.161229841765317  -0.408248290463863  -0.464242826880013  -0.303012985114696;
 403        0.408248290463863   0.408248290463863                   0  -0.408248290463863  -0.408248290463863                   0   0.408248290463863   0.408248290463863;
 404        0.464242826880013   0.161229841765317  -0.408248290463863  -0.303012985114696   0.303012985114696   0.408248290463863  -0.161229841765317  -0.464242826880013;
 405        0.464242826880013  -0.161229841765317  -0.408248290463863   0.303012985114696   0.303012985114696  -0.408248290463863  -0.161229841765317   0.464242826880013;
 406        0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863;
 407        0.303012985114696  -0.464242826880013   0.408248290463863  -0.161229841765317  -0.161229841765317   0.408248290463863  -0.464242826880013   0.303012985114696;
 408        0.161229841765317  -0.303012985114696   0.408248290463863  -0.464242826880013   0.464242826880013  -0.408248290463863   0.303012985114696  -0.161229841765317];
 409 elseif (N == 7) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
 410     Tforward =[ 0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227;
 411        0.521120889169602   0.417906505941275   0.231920613924330                   0  -0.231920613924330  -0.417906505941275  -0.521120889169602;
 412        0.481588117120063   0.118942442321354  -0.333269317528993  -0.534522483824849  -0.333269317528993   0.118942442321354   0.481588117120063;
 413        0.417906505941275  -0.231920613924330  -0.521120889169602                   0   0.521120889169602   0.231920613924330  -0.417906505941275;
 414        0.333269317528993  -0.481588117120063  -0.118942442321354   0.534522483824849  -0.118942442321354  -0.481588117120063   0.333269317528993;
 415        0.231920613924330  -0.521120889169602   0.417906505941275                   0  -0.417906505941275   0.521120889169602  -0.231920613924330;
 416        0.118942442321354  -0.333269317528993   0.481588117120063  -0.534522483824849   0.481588117120063  -0.333269317528993   0.118942442321354];
 417 elseif strcmp(transform_type, 'dct') == 1,
 418     Tforward    = dct(eye(N));
 419 elseif strcmp(transform_type, 'dst') == 1,
 420     Tforward    = dst(eye(N));
 421 elseif strcmp(transform_type, 'DCrand') == 1,
 422     x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x);
 423     if (Q(1) < 0),
 424         Q = -Q;
 425     end;
 426     Tforward = Q';
 427 else %% a wavelet decomposition supported by 'wavedec'
 428     %%% Set periodic boundary conditions, to preserve bi-orthogonality
 429     dwtmode('per','nodisp');
 430
 431     Tforward = zeros(N,N);
 432     for i = 1:N
 433         Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
 434     end
 435 end
 436
 437 %%% Normalize the basis elements
 438 Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))';
 439
 440 %%% Compute the inverse transform matrix
 441 Tinverse = inv(Tforward);
 442
 443 return;
 444
 445