[ Pobierz całość w formacie PDF ]
It is in principle possible to use a similar approach for longer filters, but ironically, computing the initial
condition itself requires filtering. To implement noncausal filtering and filtering with other boundary
handling, it is usually fastest to pad the input signal, applyfilter, and then truncate the result.
Upsampling and Downsampling
Upsamplex(insert zeros) by factorp:
y = zeros(length(x)*p-p+1,1); % For trailing zeros, use y = zeros(length(x)*p,1);
y(1:p:length(x)*p) = x;
Downsamplexby factorp, where 1 d"qd"p:
y = x(q:p:end);
27
Haar Wavelet
This code performsKstages of the Haar wavelet transform on the input dataxto produce wavelet
coefficientsy. The input arrayxmust have length divisible by 2K.
Forward transform: Inverse transform:
y = x(:); N = length(y); x = y(:); N = length(x)*pow2(-K);
for k = 1:K for k = 1:K
tmp = y(1:2:N) + y(2:2:N); N = N*2;
y([1:N/2,N/2+1:N]) = ... tmp = x(N/2+1:N) + 0.5*x(1:N/2);
[tmp;y(1:2:N)- 0.5*tmp]/sqrt(2); x([1:2:N,2:2:N]) = ...
N = N/2; [tmp;x(1:N/2)- tmp]*sqrt(2);
end end
Daubechies 4-Tap Wavelet
This code implements the Daubechies 4-tap wavelet in lifting scheme form [2]. The input arrayxmust
have length divisible by 2K. Filtering is done with symmetric boundary handling.
Forward transform: Inverse transform:
U = 0.25*[2-sqrt(3),-sqrt(3)]; U = 0.25*[2-sqrt(3),-sqrt(3)];
ScaleS = (sqrt(3)- 1)/sqrt(2); ScaleS = (sqrt(3)- 1)/sqrt(2);
ScaleD = (sqrt(3) + 1)/sqrt(2); ScaleD = (sqrt(3) + 1)/sqrt(2);
y = x(:); N = length(y); x = y(:); N = length(x)*pow2(-K);
for k = 1:K for k = 1:K
N = N/2; y1 = x(N+1:2*N)/ScaleD;
y1 = y(1:2:2*N); y2 = x(1:N)/ScaleS + y1([min(2,N),1:N-1]);
y2 = y(2:2:2*N) + sqrt(3)*y1; y1 = y1- filter(U,1,...
y1 = y1 + filter(U,1,... y2([2:N,max(N-1,1)]),y2(1)*U(2));
y2([2:N,max(N-1,1)]),y2(1)*U(2)); x([1:2:2*N,2:2:2*N]) = [y1;y2- sqrt(3)*y1];
y(1:2*N) = [ScaleS*... N = 2*N;
(y2- y1([min(2,N),1:N-1]));ScaleD*y1]; end
end
To use periodic boundary handling rather than symmetric boundary handling, change appearances of
[2:N,max(N-1,1)]to[2:N,1]and[min(2,N),1:N-1]to[N,1:N-1].
28
Miscellaneous
11
Clip a value without using if statements
To clip (or clamp, bound, or saturate) a value to within a range, the straightforward implementation is
if x
x = lowerBound;
elseif x> upperBound
x = upperBound;
end
Unfortunately, this is slow. A faster method is to use theminandmaxfunctions:
x = max(x,lowerBound); % Clip elements from below, force x>= lowerBound
x = min(x,upperBound); % Clip elements from above, force x
Moreover, it works per-element ifxa matrix of any size.
Convert any array into a column vector
It is often useful to force an array to be a column vector, for example, when writing a function expecting
a column vector as an input. This simple trick will convert the input array of any shape and number
of dimensions to a column vector, even if the input is already a column vector.
x = x(:); % Convert x to a column vector
Similarly,x = x(:,:)reduces an N-d array to 2D, where the number of rows stays the same.
Find the min/max of a matrix or N-d array
Given a matrix input (with no other inputs), theminandmaxfunctions operate along the columns,
finding the extreme element in each column. To find the extreme element over the entire matrix, one
way for a two-dimensional matrix ismin(min(A)). A better method that works regardless of the number
of dimensions and determines the extreme element s location is
[MinValue, MinIndex] = min(A(:)); % Find the minimum element in A
% The minimum value is MinValue, the index is MinIndex
MinSub = ind2sub(size(A), MinIndex); % Convert MinIndex to subscripts
The minimum element isA(MinIndex)orA(MinSub(1), MinSub(2), . . .)as a subscript reference.
(Similarly, replaceminwithmaxfor maximum value.)
29
Flood filling
Flood filling, like the bucket tool in image editors, can be elegantly written as a recursive function:
function I = flood1(I,c,x,y)
% Flood fills image I from point (x,y) with color c.
c2 = I(y,x);
I(y,x) = c;
if x> 1 & I(y,x-1) == c2 & I(y,x-1) Ü= c, I = flood1(I,c,x-1,y); end
if x
if y> 1 & I(y-1,x) == c2 & I(y-1,x) Ü= c, I = flood1(I,c,x,y-1); end
if y
Being a highly recursive function, this is inefficient in Matlab. The following code is faster:
function I = flood2(I,c,x,y)
[ Pobierz całość w formacie PDF ]