setf05.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function ErrCode=SetF05(FrameFile,FIRspec,a)
% SetF05    Set F as a general filter bank, using initial values from FIR filters
%           The FrameFile will be of Type 'g'
%
% ErrCode=SetF05(FrameFile,FIRspec);
% ErrCode=SetF05(FrameFile,FIRspec,a);
% -----------------------------------------------------------------------------------
% Arguments:
%  ErrCode   - 0 is returned if the function execute without error
%              1 an 'unexpected' error occurred, lasterr may explain it
%              2 or larger, another error occurred, 
%  FrameFile - the name of the mat-file used to store the frame
%  FIRspec   - the specifications for the filters, a matrix of size Jx4
%              where J is number of different filters to use.
%              The filters are designed using the fir1 function
%              FIRspec(j,1) is Lj, length (=order+1) of filter j
%              FIRspec(j,2) is nj, the upsampling factor for filter j
%              FIRspec(j,3) is w1, the lower frequency for the passband of 
%                 filter j, if w1 is 0 we have a low-pass filter
%              FIRspec(j,4) is w2, the higher frequency for the passband of 
%                 filter j, if w1 is 1 we have a high-pass filter
%              If w1=0 and w2=1 the filter coefficients are just filled with 
%              random values.
%  a         - indicate the kind of structure imposed on the filter bank
%              this affect the G matrix
%              0  - all non-zero variables are free, but zeros are fixed
%              a  - a is a real number and 0<a<1, then [F,G]=BuildG(Fg,a)
%                   a should be small!
%              1  - all filter coefficients are free, but equal filters are 
%                   kept equal.
%              a  - a is a real number and 1<a<2, then BuildG(Fg,a-1) is used
%                   but filters are kept equal, (a-1) should be small!
% -----------------------------------------------------------------------------------

%----------------------------------------------------------------------
% Copyright (c) 2001.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  11.01.2001  KS: function made
% Ver. 1.1  03.12.2002  KS: moved from ..\Frames to ..\FrameTools
%----------------------------------------------------------------------

Mfile='SetF05';
ErrCode=0; 

if (nargin < 3); a=0; end;   % the default
if (nargin < 2)
   disp([Mfile,': wrong number of arguments, see help.']);
   ErrCode=2;
   return
end

try
   if exist([FrameFile,'.mat'])
      disp([Mfile,': ',FrameFile,'.mat already exists, it will be overwritten.']);
   end
catch
   ErrCode=3;
   return
end

% initialize the variable in FrameFile
Class='';
Type='g';
Mdim=1;
F=[];
G=[];Dtab=[];Ctab=[];  
Fbest=F;
Savg=0;
Mdat=0;
PreProc=struct('Prog1','','Prog2','','Method',[],'arg1',[],'arg2',[],'arg3',[]);
VecSel=struct('Prog1','','arg1',[],'arg2',[],'arg3',[],'arg4',[],'arg5',[]);
InitialF=struct('Prog1',Mfile,'arg1',FrameFile,'arg2',FIRspec,'arg3',a,...
   'arg4',[],'arg5',[]);
History='';
SNRtot=[];

% check FIRspec
[J,temp]=size(FIRspec);
if temp~=4
   disp([Mfile,': illegal size of Firspec.']);
   ErrCode=4;
   return
end
if (min(FIRspec(:,1))<2)
   disp([Mfile,': a short filter is given in Firspec.']);
   ErrCode=5;
   return
end

% find N, K and P
N=FIRspec(1,2);       % N is least common multiple of upsampling factors
for j=2:J
   N=lcm(N,FIRspec(j,2));
end
K=sum(N./FIRspec(:,2));
P=max(ceil((FIRspec(:,1)-FIRspec(:,2))./N))+1;  % overlap factor
disp(['Some values are  N=',int2str(N),' K=',int2str(K),' P=',int2str(P)]);

%   make Fg
k=0;
Fg=zeros(N*P,K);
for j=1:J
   Lj=FIRspec(j,1);    % length of this filter
   nj=FIRspec(j,2);    % upsampling rate for this filter
   w1=FIRspec(j,3);    
   w2=FIRspec(j,4);    
   if (w1==0) & (w2>0) & (w2<1)
      f=fir1(Lj-1,w2);                % lowpass filter
   elseif (w1>0) & (w1<w2) & (w2<1)
      f=fir1(Lj-1,[w1,w2]);           % bandpass filter
   elseif (w2==1) & (w1>0) & (w1<1)   
      f=fir1(Lj-1,w1,'high');         % highpass filter
   else
      f=randn(1,Lj);                  % random values
   end
   f=f(1:Lj)';
   temp=f'*f;
   if (temp>0); f=f/sqrt(temp);  end;
   %
   n=floor((N*P-N+nj-Lj)/2);
   n=n+(1:Lj);   % the indexes to place f into
   for i=1:(N/nj)
      % place the filter into Fg
      k=k+1;
      Fg(n,k)=f;
      n=n+nj;
   end
end

% Make G
if a<1
   [F,G]=BuildG(Fg,a);          % F is now Qx1, and G is NPxK
elseif a<2  
   % let each filter be equal
   G=zeros(size(Fg));
   F=[];Q=0;
   k1=0;k2=0;
   for j=1:J
      k1=k2+1;                  % last value of k2
      k2=k2+N/FIRspec(j,2);     % new value of k2
      [f,G(:,k1:k2)]=BuildG(Fg(:,k1:k2),a-1);
      if a==1
         q=size(f,1)/(k2-k1+1);
      else
         q=size(f,1);
      end
      if rem(q,1)
         disp([Mfile,': non-integer value of q.']);
         q=ceil(q);
      end
      F=[F;f(1:q)];
      G(:,k1:k2)=sign(G(:,k1:k2)).*(abs(G(:,k1:k2))+Q);
      Q=Q+q;
      I=find(abs(G)>Q);
      while length(I)
         G(I)=sign(G(I)).*(abs(G(I))-q);
         I=find(abs(G)>Q);
      end
   end
else
   disp([Mfile,': strange value of fourth argument, a.']);
   ErrCode=6;
   return
end

[Q,temp]=size(F);
if temp~=1
   disp([Mfile,': possible logical error, size(F) is not correct']);
   ErrCode=7;
   return
end
if max(G(:))~=Q
   disp([Mfile,': possible logical error, max(G(:))~=Q.']);
   ErrCode=8;
   return
end

% reshape G from NPxK to NxKxP
temp=G;G=zeros(N,K,P);
for p=1:P; G(:,:,p)=temp((1:N)+(p-1)*N,:); end;

% find Ctab and Dtab and SizeF
[Ctab,Dtab]=MapCD(G);
SizeF=size(G);
Fbest=F;
History=char([Mfile,' ',datestr(now),', initialized the frame ',FrameFile],...
   ['using ',int2str(J),' different FIR filters, giving N=',...
    int2str(N),' K=',int2str(K),' P=',int2str(P),' Q=',int2str(Q)]);
% save the frame in a file 
try
   save(FrameFile,'Class','Type','Mdim','F','SizeF','G','Ctab','Dtab',...
        'Fbest','Savg','Mdat','PreProc','VecSel','InitialF','History','SNRtot');
catch
   disp([Mfile,': error saving ',FrameFile,'.']);
   ErrCode=9;
   return
end

return


% testing
clear all
FrameFile='test';
FIRspec=[40,8,0,0.05; 24,8,0,0.25; 20,4,0.20,0.50; 12,2,0.5,0.95];a=0;
ErrCode=SetF05(FrameFile,FIRspec,a);
load(FrameFile);
Fg=BuildFg(F,G);
disp(['The value of Q is ',int2str(max(G(:))),'==',int2str(size(F,1))]);
SizeF
PlotF(Fg);

% alternative filter specifications
FIRspec=[40,8,0,0.08; 32,8,0,0.15; 32,8,0.1,0.25; 
         20,8,0,0.4; 16,4,0.25,0.75; 16,4,0.2,0.5; 12,2,0.4,0.8; 9,2,0.4,1];
% a simple filter bank with Q=30 (all filters are symmetric)      
FIRspec=[24,8,0,0.05; 16,8,0,0.25; 12,4,0.25,0.75; 7,2,0.25,1];a=1+1e-10;
% a filter bank with P=5, Q=152 (all filters are free)      
FIRspec=[40,8,0,0.05; 24,8,0,0.25; 20,4,0.20,0.50; 12,2,0.5,0.95];a=0;
      
% a filter bank with P=5, Q= (all filters are free)      
FIRspec=[40,8,0,0.05; 40,8,0,0.15; 40,8,0.15,0.20; 40,8,0.25,0.30; ...
      24,4,0,0.25; 24,4,0.45,0.55;      16,4,0,0.5; 16,4,0.4,0.8; 8,2,0.5,1];
a=1;