## [a9068d]: gabor / gabphasegrad.m Maximize Restore History

 ``` 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``` ```function [tgrad,fgrad,c]=gabphasegrad(method,varargin) %GABPHASEGRAD Phase gradient of the DGT % Usage: [tgrad,fgrad,c] = gabphasegrad('dgt',f,g,a,M); % [tgrad,fgrad] = gabphasegrad('phase',cphase,a); % [tgrad,fgrad] = gabphasegrad('abs',s,g,a); % % `[tgrad,fgrad]=gabphasegrad(method,...)` computes the time-frequency % gradient of the phase of the |dgt| of a signal. The derivative in time % *tgrad* is the instantaneous frequency while the frequency derivative % *fgrad* is the local group delay. % % *tgrad* and *fgrad* measure the deviation from the current time and % frequency, so a value of zero means that the instantaneous frequency is % equal to the center frequency of the considered channel. % % *tgrad* is scaled such that distances are measured in samples. Similarly, % *fgrad* is scaled such that the Nyquest frequency (this highest possible % frequency) corresponds to a value of L/2. % % The computation of *tgrad* and *fgrad* is inaccurate when the absolute % value of the Gabor coefficients is low. This is due to the fact the the % phase of complex numbers close to the machine precision is almost % random. Therefore, *tgrad* and *fgrad* may attain very large random values % when `abs(c)` is close to zero. % % The computation can be done using three different methods. % % 'dgt' Directly from the signal. This is the default method. % % 'phase' From the phase of a DGT of the signal. This is the % classic method used in the phase vocoder. % % 'abs' From the absolute value of the DGT. Currently this % method works only for Gaussian windows. % % `[tgrad,fgrad]=gabphasegrad('dgt',f,g,a,M)` computes the time-frequency % gradient using a DGT of the signal *f*. The DGT is computed using the % window *g* on the lattice specified by the time shift *a* and the number % of channels *M*. The algorithm used to perform this calculation computes % several DGTs, and therefore this routine takes the exact same input % parameters as |dgt|. % % The window *g* may be specified as in |dgt|. If the window used is % 'gauss', the computation will be done by a faster algorithm. % % `[tgrad,fgrad,c]=gabphasegrad('dgt',f,g,a,M)` additionally returns the % Gabor coefficients *c*, as they are always computed as a byproduct of the % algorithm. % % `[tgrad,fgrad]=gabphasegrad('phase',cphase,a)` computes the phase % gradient from the phase *cphase* of a DGT of the signal. The original DGT % from which the phase is obtained must have been computed using a % time-shift of *a*. % % `[tgrad,fgrad]=gabphasegrad('abs',s,g,a)` computes the phase gradient % from the spectrogram *s*. The spectrogram must have been computed using % the window *g* and time-shift *a*. % % `[tgrad,fgrad]=gabphasegrad('abs',s,g,a,difforder)` uses a centered finite % diffence scheme of order *difforder* to perform the needed numerical % differentiation. Default is to use a 4'th order scheme. % % Currently the `'abs'` method only works if the window *g* is a Gaussian % window specified as a string or cell array. % % See also: resgram, gabreassign, dgt % % References: aufl95 cmdaaufl97 fl65 % AUTHOR: Peter L. S��ndergaard, 2008. %error(nargchk(4,6,nargin)); if ~ischar(method) error(['First argument must be the method name, "dgt", "phase" or ' ... '"abs".']); end; switch(lower(method)) case 'dgt' % --------------------------- DGT method ------------------------ [f,g,a,M]=deal(varargin{1:4}); definput.keyvals.L=[]; definput.keyvals.minlvl=eps; definput.keyvals.lt=[0 1]; [flags,kv,L,minlvl]=ltfatarghelper({'L','minlvl'},definput,varargin(5:end)); %% ----- step 1 : Verify f and determine its length ------- % Change f to correct shape. [f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,upper(mfilename),0); %% ------ step 2: Verify a, M and L if isempty(L) % ----- step 2b : Verify a, M and get L from the signal length f---------- L=dgtlength(Ls,a,M,kv.lt); else % ----- step 2a : Verify a, M and get L Luser=dgtlength(L,a,M,kv.lt); if Luser~=L error(['%s: Incorrect transform length L=%i specified. Next valid length ' ... 'is L=%i. See the help of DGTLENGTH for the requirements.'],... upper(mfilename),L,Luser); end; end; %% ----- step 3 : Determine the window [g,info]=gabwin(g,a,M,L,kv.lt,'callfun',upper(mfilename)); if L3 difforder=varargin{4}; else difforder=4; end; if ~(all(s(:)>=0)) error('First input argument must be positive or zero.'); end; [M,N,W]=size(s); L=N*a; tfr=1; g [g,info]=gabwin(g,a,M,L,'callfun','GABPHASEGRAD'); info if ~info.gauss error(['The window must be a Gaussian window (specified as a string or ' ... 'as a cell arrray).']); end; L=N*a; b=L/M; % We must avoid taking the log of zero. % Therefore we add the smallest possible % number logs=log(s+realmin); % XXX REMOVE Add a small constant to limit the dynamic range. This should % lessen the problem of errors in the differentation for points close to % (but not exactly) zeros points. maxmax=max(logs(:)); tt=-11; logs(logs