Writing Fast Matlab Code #10:Signal Processing
Filed under: Code Optimization, Signal Processing, Tutorials
Even without the Signal Processing Toolbox, Matlab is quite capable in signal processing computa- tions. This section lists code snippets to perform several common operations efficiently.
Moving-average filter
To compute an N-sample moving average of x with zero padding:
y = filter(ones(N,1)/N,1,x);
For large N, it is faster to use
y = cumsum(x)/N;
y(N+1:end) = y(N+1:end) − y(1:end−N);
Locating zero-crossings and extrema
To obtain the indices where signal x crosses zero:
i = find(diff(sign(x))); % The kth zero−crossing lies between x(i(k)) and x(i(k)+1)
Linear interpolation can be used for subsample estimates of zero-crossing locations:
i = find(diff(sign(x))); i = i − x(i)./(x(i+1) − x(i)); % Linear interpolation
Since local maximum and minimum points of a signal have zero derivative, their locations can be estimated from the zero-crossings of diff(x), provided the signal is sampled with sufficiently fine resolution. For a coarsely sampled signal, a better estimate is
iMax = find(sign(x(2:end−1)−x(1:end−2)) + sign(x(2:end−1)−x(3:end)) > 0) + 1; iMin = find(sign(x(2:end−1)−x(1:end−2)) + sign(x(2:end−1)−x(3:end)) < 0) + 1;
FFT-based convolution
This line performs FFT-based circular convolution, equivalent to y = filter(b,1,x) except near the boundaries, provided that length(b) < length(x):
y = ifft(fft(b,length(x)).*fft(x));
For FFT-based zero-padded convolution, equivalent to y = filter(b,1,x),
N = length(x)+length(b)−1; y = ifft(fft(b,N).*fft(x,N)); y = y(1:length(x));
In both code snippets above, y will be complex-valued, even if x and b are both real. To force y to be real, follow the computation with y = real(y). If you have the Signal Processing Toolbox, it is faster to use fftfilt for FFT-based, zero-padded filtering.
Noncausal filtering and other boundary extensions
For its intended use, the filter command is limited to causal filters, that is, filters that do not involve “future” values to the right of the current sample. Furthermore, filter is limited to zero-padded boundary extension; data beyond the array boundaries are assumed zero as if the data were padded with zeros.
For two-tap filters, noncausal filtering and other boundary extensions are possible through filter’s fourth initial condition argument. Given a boundary extension value padLeft for x(0), the filter yn = b1 xn + b2 xn−1 may be implemented as
y = filter(b,1,x,padLeft*b(2));
Similarly, given a boundary extension value padRight for x(end+1), the filter yn = b1 xn+1 + b2 xn is implemented as
y = filter(b,1,[x(2:end),padRight],x(1)*b(2));
Choices for padLeft and padRight for various boundary extensions are
| Boundary extension | padLeft | padRight |
| Periodic
Whole-sample symmetric Half-sample symmetric Antisymmetric |
x(end) x(2) x(1)
2*x(1)-x(2) |
x(1) x(end-1) x(end)
2*x(end)-x(end-1) |
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, apply filter, and then truncate the result.
Upsampling and Downsampling
Upsample x (insert zeros) by factor p:
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;
Downsample x by factor p, where 1 ≤ q ≤ p:
y = x(q:p:end);
Haar Wavelet
This code performs K stages of the Haar wavelet transform on the input data x to produce wavelet coefficients y. The input array x must have length divisible by 2K.
Forward transform:
y = x(:); N = length(y); for k = 1:K tmp = y(1:2:N) + y(2:2:N); y([1:N/2,N/2+1:N]) = ... [tmp;y(1:2:N) − 0.5*tmp]/sqrt(2); N = N/2; end
Inverse transform:
x = y(:); N = length(x)*pow2(−K); for k = 1:K N = N*2; tmp = x(N/2+1:N) + 0.5*x(1:N/2); x([1:2:N,2:2:N]) = ... [tmp;x(1:N/2) − tmp]*sqrt(2); end
Daubechies 4-Tap Wavelet
This code implements the Daubechies 4-tap wavelet in lifting scheme form [2]. The input array x must have length divisible by 2K. Filtering is done with symmetric boundary handling.
Forward transform:
U = 0.25*[2−sqrt(3),−sqrt(3)]; ScaleS = (sqrt(3) − 1)/sqrt(2); ScaleD = (sqrt(3) + 1)/sqrt(2); y = x(:); N = length(y); for k = 1:K N = N/2; y1 = y(1:2:2*N); y2 = y(2:2:2*N) + sqrt(3)*y1; y1 = y1 + filter(U,1,... y2([2:N,max(N−1,1)]),y2(1)*U(2)); y(1:2*N) = [ScaleS*... (y2 − y1([min(2,N),1:N−1]));ScaleD*y1]; end
Inverse transform:
U = 0.25*[2−sqrt(3),−sqrt(3)]; ScaleS = (sqrt(3) − 1)/sqrt(2); ScaleD = (sqrt(3) + 1)/sqrt(2); x = y(:); N = length(x)*pow2(−K); for k = 1:K y1 = x(N+1:2*N)/ScaleD; y2 = x(1:N)/ScaleS + y1([min(2,N),1:N−1]); y1 = y1 − filter(U,1,... 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]; N = 2*N; 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].
Popularity: 1% [?]
Writing Fast Matlab Code #9:Numerical Integration
Quadrature formulas are numerical approximation of integrals of the form
where the xk are called the nodes or abscissas and wk are the associated weights. Simpson’s rule is
Simpson’s rule is a quadrature formula with nodes a, a + h, b and node weights h/3,4/3h, h/3 .
Matlab has several functions for quadrature in one dimension:
| quad | adaptive Simpson | Better for low accuracy and nonsmooth integrands |
| quadl | adaptive Gauss-Lobatto | Higher accuracy with smoother integrands |
| quadgk | adaptive Gauss-Kronrod | Oscillatory integrands and high accuracy |
| quadv | adaptive Simpson | Vectorized for multi-valued integrands |
The quad- functions are robust and precise, however, they are not very efficient. They use an adaptive refinement procedure, and while this reduces the number of function calls, they gain little from vectorization and incur significant overhead.
If an application requires approximating an integral whose integrand can be efficiently vectorized, using nonadaptive quadrature may improve speed.
for n = −20:20 % Compute Fourier series coefficients c(n + 21) = quad(@(x)(exp(sin(x).ˆ6).*exp(−1i*x*n)),0,pi,1e−4); end
This code runs in 5.16 seconds. In place of quad, using Simpson’s composite rule with N = 199 nodes yields results with comparable accuracy and allows for vectorized computation. Since the integrals are all over the same interval, the nodes and weights only need to be constructed once.
N = 199; h = pi/(N−1); x = (0:h:pi).'; % Nodes w = ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; % Weights for n = −20:20 c(n + 21) = w * ( exp(sin(x).ˆ6).*exp(−1i*x*n) ); end
This version of the code runs in 0.02 seconds (200 times faster). The quadrature is performed by the dot product multiplication with w. It can be further optimized by replacing the for loop with one vector-matrix multiply:
[n,x] = meshgrid(−20:20, 0:h:pi); c = w * ( exp(sin(x).ˆ6).*exp(−1i*x.*n) );
For this example, quadv can be used on a multi-valued integrand with similar accuracy and speed,
n = −20:20;
c = quadv(@(x)exp(sin(x).ˆ6).*exp(−1i.*x.*n),0,pi,1e−4);
9.1 One-Dimensional Integration
|
|
is approximated by composite Simpson’s rule with
h = (b − a)/(N−1); x = (a:h:b).'; w = ones(1,N); w(2:2:N−1) = 4; w(3:2:N−2) = 2; w = w*h/3; I = w * f(x); % Approximately evaluate the integral
where N is an odd integer specifying the number of nodes.
A good higher-order choice is 4th-order Gauss-Lobatto [4] (as used by quadl), based on
N = max(3*round((N−1)/3),3) + 1; % Adjust N to the closest valid choice h = (b − a)/(N−1); d = (3/sqrt(5) − 1)*h/2; x = (a:h:b).'; x(2:3:N−2) = x(2:3:N−2) − d; x(3:3:N−1) = x(3:3:N−1) + d; w = ones(1,N); w(4:3:N−3) = 2; w([2:3:N−2,3:3:N−1]) = 5; w = w*h/4; I = w * f(x); % Approximately evaluate the integral
The number of nodes N must be such that (N − 1)/3 is an integer. If not, the first line adjusts N to the closest valid choice. It is usually more accurate than Simpson’s rule when f has six continuous derivatives, f ∈ C 6 (a, b).
A disadvantage of this nonadaptive approach is that the accuracy of the result is only indirectly con- trolled by the parameter N. To guarantee a desired accuracy, either use a generously large value for N or, if possible, determine the error bounds 
[6,7].
where h = (b-a)/(N-1). Note that these bounds are valid only when the integrand is sufficiently differentiable: f N −1 , must have four continuous derivatives for the Simpson’s rule error bound, and six continuous derivatives for 4th-order Gauss-Lobatto.
Composite Simpson’s rule is a fast and effective default choice. But depending on the situation, other methods may be faster and more accurate:
- If the integrand is expensive to evaluate, or if high accuracy is required and the error bounds above are too difficult to compute, use one of Matlab’s adaptive methods. Otherwise, consider composite methods.
- Use higher-order methods (like Gauss-Lobatto/quadl) for very smooth integrands and lower-order methods for less smooth integrands.
- Use the substitution u = 1/(1-x) or Gauss-Laguerre quadrature for infinite-domain integrals like f ∞.
9.2 Multidimensional Integration
An approach for evaluating double integrals of the form 
%%% Construct Simpson nodes and weights over x %%% h = (b − a)/(Nx−1); x = (a:h:b).'; wx = ones(1,Nx); wx(2:2:Nx−1) = 4; wx(3:2:Nx−2) = 2; wx = w*h/3; %%% Construct Simpson nodes and weights over y %%% h = (d − c)/(Ny−1); y = (c:h:d).'; wy = ones(1,Ny); wy(2:2:Ny−1) = 4; wy(3:2:Ny−2) = 2; wy = w*h/3; %%% Combine for two−dimensional integration %%% [x,y] = meshgrid(x,y); x = x(:); y = y(:); w = wy.'*wx; w = w(:).'; I = w * f(x,y); % Approximately evaluate the integral
Similarly for three-dimensional integrals, the weights are combined with
[x,y,z] = meshgrid(x,y,z); x = x(:); y = y(:); z = z(:); w = wy.'*wx; w = w(:)*wz; w = w(:).';
|
|
When the integration region is complicated or of high dimension, Monte Carlo integration techniques are appropriate. The disadvantage is that an N -point Monte Carlo quadrature has error on the order O(1/sqrt(n)), so many points are necessary even for moderate accuracy. Suppose that N points,x1 , x2 , . . . , xN , are uniformly randomly selected in a multidimensional region (or volume) Ω. Then
To integrate a complicated region W that is difficult to sample uniformly, find an easier region Ω that contains W and can be sampled [5]. Then
χW (x) is the indicator function of W : χW (x) = 1 when x is within region W and χW (x) = 0 otherwise. Multiplying the integrand by χW sets contributions from outside of W to zero.
%%% Uniformly randomly sample points (x,y) in Ω %%% x = 4*rand(N,1) − 2; y = 4*rand(N,1) − 2; %%% Restrict the points to region W %%% i = (cos(2*sqrt(x.ˆ2 + y.ˆ2)).*x <= y & x.ˆ2 + y.ˆ2 <= 4); x = x(i); y = y(i); %%% Approximately evaluate the integrals %%% area = 4*4; % The area of rectangle Ω M = (area/N) * length(x); Mx = (area/N) * sum(x); My = (area/N) * sum(y);
Region W sampled with N = 1500. The center of mass is ≈ (0.5, 0.7).
More generally, if W is a two-dimensional region contained in the rectangle defined by a ≤ x ≤ b and
c ≤ y ≤ d, the following code approximates 
x = a + (b−a)*rand(N,1); y = c + (d−c)*rand(N,1); i = logical(indicatorW(x,y)); x = x(i); y = y(i); area = (b−a)*(d−c); I = (area/N) * sum(f(x,y)); % Approximately evaluate the integral
where indicatorW(x,y) is the indicator function χW (x, y) for region W .
For refinements and variations of Monte Carlo integration, see for example [1].
Popularity: 1% [?]
Writing Fast Matlab Code #8: Solving a linear system (Ax = b)
Many problems involve solving a linear system, that is, solving the Ax = b problem
The condition number cond(A) describes how sensitive the problem is to small perturbations. Matlab can efficiently approximate cond(A) with condest. Smaller condition number is better, and in this case A is said to be well-conditioned. An infinite condition number implies that A is singular—the problem cannot be solved, though it may be solved approximately with the pseudoinverse :
x = pinv(A)*b
Small- to moderately-sized problems are solved efficiently by the backslash operator:
x = A\b; % Solves A*x = bWe shall focus on large problems.
Iterative Methods
For large problems, Matlab is well-equipped with iterative solvers: bicg, bicgstab, cgs, gmres, lsqr,
minres, pcg, symmlq, qmr. The general calling syntax (except for gmres) is
[x,flag,relres,iter,resvec] = method (A,b,tol,maxit,M1,M2,x0)
where the required arguments are the matrix A, the right-hand side b, and the solution x. The minres, pcg, and symmlq methods require A to be symmetric (A = A’) and pcg additionally requires positive definiteness (all eigenvalues are positive). The optional input variables are
- tol accuracy tolerance (default 10−6 )
- maxit maximum number of iterations (default 20)
- M1, M2 the method uses preconditioner M = M1*M2; it solves M −1 Ax = M −1 b
- x0 initial guess
The optional output variables are
- flag convergence flag (0 = success, nonzero values are error codes)
- relres relative residual kb − Axk/kbk
- iter number of iterations
- resvec a vector listing the residual norm kb − Axk at each iteration
The gmres method includes one additional input,
[x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,...)
The method restarts every restart iterations, or if restart=[], the method does not restart.
Which method works best depends on the particular problem. This diagram (adapted from Demmel [3])
provides a reasonable starting point.
Functional Representation
Rather than storing the matrix explicitly, it is often convenient to represent A as a linear function acting on x. For example,
can be expressed as A*x = cumsum(x). In addition to being conceptually convenient, functional repre- sentation avoids the memory overhead of storing the matrix.
All the iterative solvers support functional representation, where a function handle afun is given in place of matrix A. The function afun(x) should accept a vector x and return A*x. The help for lsqr (“help lsqr”) has a functional example where A is the central differencing operator. The bicg, lsqr, and qmr methods also require the transpose: afun(x,’notransp’) should return A*x and afun(x,’transp’) should return A’*x. Similarly, the preconditioner may be given as a function handle mfun.
Special Matrices
For certain special matrix structures, the Ax = b problem can be solved very quickly.
- Circulant matrices. Matrices corresponding to cyclic convolution Ax = h ∗ x are diagonalized in the Fourier domain, and can be solved by x = ifft(fft(b)./fft(h));
- Triangular and banded. Triangular matrices and diagonally-dominant banded matrices are solved efficiently by sparse LU factorization: [L,U] = lu(sparse(A)); x = U\(L\b);
- Poisson problem. If A is a finite difference approximation of the Laplacian, the problem is effi- ciently solved by multigrid methods (e.g., web search for “matlab multigrid”).
Inlined PCG
As discussed in Section 6, it is sometimes advantageous to inline a simple piece of code to tailor it to a specific application. While most iterative methods are too lengthy for inlining, conjugate gradients has a short implementation. The value of inlining pcg is that calls to afun and mfun may be inlined. Beware that in the following code, A must be symmetric positive definite—otherwise, it will not work!
%% Use PCG to solve Ax = b %% tol = 1e−6; % Convergence tolerance maxit = 20; % Maximum number of iterations x = x0; % Set x to initial guess (if no guess is available, use x0 = 0) Compute w = A x r = b − w; % r = residual vector Set p = r, or if using a preconditioner, p = M −1 r d = r'*p; for iter = 1:maxit if norm(p) < tol break; % PCG converged successfully (small change in x) end Compute w = A p alpha = e/(p'*w); x = x + alpha*p; r = r − alpha*w; Set w = r, or if using a preconditioner, w = M −1 r if norm(w) < tol && norm(r) < tol break; % PCG converged successfully (small residual) end end dlast = d; d = r'*w; p = w + (d/dlast)*p;
After the loop, the final solution is in x. The method converged successfully with accuracy tol if it encountered a break statement, otherwise, it failed or needed more iterations. The relative residual may be computed by relres = norm(r)/norm(b) after the loop and resvec by storing the value of norm(r) each iteration.
Popularity: 1% [?]
Translator
Follow Us
Video Posts
