[580282]: wavelets / fwtinit.m Maximize Restore History

Download this file

fwtinit.m    307 lines (269 with data), 8.1 kB

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
function [w] = fwtinit(wavname,varargin)
%FWTINIT Fast Wavelet Transform Filterbank Structure Initialization
% Usage: w = fwtinit(wavname,'a',a);
%
% Input parameters:
% wavname : Wavelet filters generating function name (without the prefix).
%
% Output parameters:
% w : Structure defining the filters.
%
% `w=fwtinit(wavname,...)` produces a structure describing the analysis
% and synthesis parts of a one level perfect reconstruction wavelet-type
% filterbank. The *wavname* can be a string defining a function name (without
% the prefix) to be called or it can be a cell-array with the first element
% beeing the function name and the remaining elements are passed further
% to the `wfilt_` function.
%
% The function is a wrapper for calling all the functions starting with
% `wfilt_` defined in the LTFAT wavelets directory. The structure which
% the function produces can be directly passed to all
% functions instead of the cell-arrays with wavelet filters which all
% `wfilt_` functions produces.
%
% The structure has the following fields:
%
% `w.filts`
% analysis or synthesis filterbank
%
% `w.h`
% analysis filter bank
%
% `w.g`
% synthesis filter bank
%
% `w.a`
% implicit subsampling factors
%
%
% Choosing wavelet filters
% ------------------------
%
%
% orthogonality/biorthogonality/frame
% number of vanishing moments psi
% symmetry/linear phase
%
% support of psi
% smoothnes and regularity of psi
%
%
%
% See also: fwt, ifwt, wfilt_db
%
% References: ma08wt
% chached last filterbank
persistent cachw;
% cached function parameters passed last
persistent cachwDesc;
% wavelet filters functions definition prefix
wprefix = 'wfilt_';
% output structure definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w.filts = {};
w.h = {};
w.g = {};
w.a = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% return empty struct if no argument was passed
if nargin<1
return;
end;
% process other parameters
definput.import = {'fwtcommon'};
% Contents of the arg_fwtcommon:
% definput.flags.ansy = {'ana','syn'};
% definput.keyvals.a = [];
[flags,kv]=ltfatarghelper({},definput,varargin);
% was the function called before with the same parameters?
% if yes, return the chached one
if(isequal(cachwDesc,{wavname,kv.a}))
w = cachw;
w = updateTransDirect(flags.do_ana,w);
return;
else
cachwDesc = {wavname,kv.a};
end
if ischar(wavname)
try
wname = parseNameValPair(wavname,wprefix);
catch err
% If failed, clean the cache.
cachwDesc = [];
cachw = [];
error(err.message);
end
elseif iscell(wavname)
wname = wavname;
if ~ischar(wname{1})
if iscell(wname{1})
if(length(wname)==1)
w = formatFilters(wname,flags.do_ana,[],w);
w = updateTransDirect(flags.do_ana,w);
cachw = w;
return;
else
error('%s: Unrecognizer format of the filterbank definition.',upper(mfilename));
end
elseif isnumeric(wname{1})
w = formatFilters(wname,flags.do_ana,[],w);
w = updateTransDirect(flags.do_ana,w);
else
error('%s: Unrecognizer format of the filterbank definition.',upper(mfilename));
end
w.a = formata(max([length(w.h),length(w.g)]),kv.a,w.a);
cachw = w;
return;
end
elseif isstruct(wavname)
if(isequal(fieldnames(wavname),fieldnames(w)))
w = wavname;
w = updateTransDirect(flags.do_ana,w);
w.a = formata(length(w.h),kv.a,w.a);
cachw = w;
return;
else
error('%s: Passed structure has different fields.',upper(mfilename));
end
else
error('%s: First argument must be a string, cell or struct.',upper(mfilename));
end;
% Search for m-file containing string wname
wfiltFiles = dir(fullfile(ltfatbasepath,sprintf('wavelets/%s%s.m',wprefix,lower(wname{1}))));
if(isempty(wfiltFiles))
error('%s: Unknown wavelet type: %s',upper(mfilename),name);
else
% if found, crop '.m' from the filename
tmpFile = wfiltFiles.name(1:end-2);
end
tmpFile
% Synthesis filters delay
d = [];
wfiltNargout = nargout(tmpFile);
if(nargin(tmpFile)~=numel(wname)-1)
error('%s: Incorrect number of parameters to be passed to the wfilt_ func.',upper(mfilename));
end
if(wfiltNargout==3)
[tmph, tmpg, w.a] = feval(tmpFile,wname{2:end});
elseif(wfiltNargout==4)
[tmph, tmpg, w.a, d] = feval(tmpFile,wname{2:end});
else
error('%s: Function %s does not return 3 or 4 arguments.',upper(mfilename),upper(tmpFile));
end
w.a = formata(length(w.h),kv.a,w.a);
w = formatFilters(tmph,1,d,w);
w = formatFilters(tmpg,0,d,w);
w = updateTransDirect(flags.do_ana,w);
cachw = w;
function w = formatFilters(cellh,do_ana,d,w)
noFilts = numel(cellh);
if(isempty(d))
d = findFiltDelays(cellh,do_ana,'half');
end
if(do_ana)
w.h = cell(noFilts,1);
for ff=1:noFilts
w.h{ff} = wfiltstruct('FIR');
w.h{ff}.h = cellh{ff};
w.h{ff}.d = d(ff);
end
else
w.g = cell(noFilts,1);
for ff=1:noFilts
w.g{ff} = wfiltstruct('FIR');
w.g{ff}.h = cellh{ff};
w.g{ff}.d = d(ff);
end
end
%END FORMATFILTERS
function w = updateTransDirect(do_ana,w)
if(do_ana)
w.filts = w.h;
else
w.filts = w.g;
end
function a = formata(filtsNo,a,origa)
if(~isempty(a))
if(length(a)==1)
atmp = ones(filtsNo,1);
a = atmp*a;
end
else
if(~isempty(origa))
a = origa;
else
atmp = ones(filtsNo,1);
a = atmp*filtsNo;
end
end
function wcell = parseNameValPair(wchar,wprefix)
%PARSENAMEVALPAIR
%Parses string in the following format wnameN1:N2... , where wname have to
%be name of the existing function with wfilt_ prefix. N1,N2,... are doubles
%delimited by character ':'.
%The output is cell array {wname,str2double(N1),str2double(N2),...}
wcharNoNum = wchar(1:find(isstrprop(wchar,'digit')~=0,1)-1);
numDelimiter = ':';
wfiltFiles = dir(fullfile(ltfatbasepath,sprintf('wavelets/%s*',wprefix)));
wfiltNames = arrayfun(@(fEl) fEl.name(1+find(fEl.name=='_',1):find(fEl.name=='.',1,'last')-1),wfiltFiles,'UniformOutput',0);
wcharMatch = cellfun(@(nEl) strcmpi(wcharNoNum,nEl),wfiltNames);
wcharMatchIdx = find(wcharMatch~=0);
if(isempty(wcharMatchIdx))
error('%s: Unknown wavelet filter definition string.',upper(mfilename));
end
if(numel(wcharMatchIdx)>1)
error('%s: Ambiguous wavelet filter definition string. Probably bug somewhere.',upper(mfilename));
end
wcell{1} = wfiltNames{wcharMatchIdx};
numString = wchar(numel(wcell{1})+1:end);
if(isempty(numString))
error('%s: No numeric parameter specified in %s.',upper(mfilename),wchar);
end
wcharNum = textscan(numString,'%f','Delimiter',numDelimiter);
if(~isnumeric(wcharNum{1})||any(isnan(wcharNum{1})))
error('%s: Incorrect numeric part of the wavelet filter definition string.',upper(mfilename));
end
wcell = [wcell, num2cell(wcharNum{1}).'];
%END PARSENAMEVALPAIR
function d = findFiltDelays(cellh,do_ana,type)
filtNo = numel(cellh);
d = ones(filtNo,1);
for ff=1:filtNo
if(strcmp(type,'half'))
if(do_ana)
d(ff) = ceil((length(cellh{ff})+1)/2);
else
d(ff) = floor((length(cellh{ff})+1)/2);
end
elseif(strcmp(type,'zeroana'))
if(do_ana)
d(ff) = 1;
else
d(ff) = length(cellh{ff});
end
elseif(strcmp(type,'zerosyn'))
if(do_ana)
d(ff) = length(cellh{ff});
else
d(ff) = 1;
end
elseif(strcmp(type,'energycent'))
tmph =cellh{ff};
tmphLen = length(tmph);
ecent = sum((1:tmphLen-1).*tmph(2:end).^2)/sum(tmph.^2);
if(do_ana)
d(ff) = round(ecent)+1;
if(rem(abs(d(ff)-d(1)),2)~=0)
d(ff)=d(ff)+1;
end
else
anad = round(ecent)+1;
d(ff) = tmphLen-anad;
if(rem(abs(d(ff)-d(1)),2)~=0)
d(ff)=d(ff)-1;
end
end
else
error('fail');
end
end