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% [?]
How to create a Kalman filter into Matlab
Download Now
KALMANF – updates a system state vector estimate based upon an observation, using a discrete Kalman filter.
Version 1.0, June 30, 2004
This tutorial function was written by Michael C. Kleder
INTRODUCTION
Many people have heard of Kalman filtering, but regard the topic as mysterious. While it’s true that deriving the Kalman filter and proving mathematically that it is “optimal” under a variety of circumstances can be rather intense, applying the filter to a basic linear system is actually very easy. This Matlab file is
intended to demonstrate that.
An excellent paper on Kalman filtering at the introductory level, without detailing the mathematical underpinnings, is:
“An Introduction to the Kalman Filter” Greg Welch and Gary Bishop, University of North Carolina
http://www.cs.unc.edu/~welch/kalman/kalmanIntro.html
PURPOSE:
The purpose of each iteration of a Kalman filter is to update the estimate of the state vector of a system (and the covariance of that vector) based upon the information in a new observation.
The version of the Kalman filter in this function assumes that observations occur at fixed discrete time intervals. Also, this function assumes a linear system, meaning that the time evolution of the state vector can be calculated by means of a state transition matrix.
USAGE:
s = kalmanf(s)
“s” is a “system” struct containing various fields used as input and output. The state estimate “x” and its covariance “P” are updated by the function. The other fields describe the mechanics of the system and are left unchanged. A calling routine may change these other fields as needed if state dynamics are time-dependent; otherwise, they should be left alone after initial values are set.
The exceptions are the observation vectro “z” and the input control (or forcing function) “u.” If there is an input function, then “u” should be set to some nonzero value by the calling routine.
SYSTEM DYNAMICS:
The system evolves according to the following difference equations, where quantities are further defined below:
- x = Ax + Bu + w meaning the state vector x evolves during one time step by premultiplying by the “state transition matrix” A. There is optionally (if nonzero) an input vector u which affects the state linearly, and this linear effect on the state is represented by premultiplying by the “input matrix” B. There is also gaussian process noise w.
- z = Hx + v meaning the observation vector z is a linear function of the state vector, and this linear relationship is represented by premultiplication by “observation matrix” H. There is also gaussian measurement noise v.
- where w ~ N(0,Q) meaning w is gaussian noise with covariance
- Q v ~ N(0,R) meaning v is gaussian noise with covariance R
VECTOR VARIABLES:
s.x = state vector estimate. In the input struct, this is the “a priori” state estimate (prior to the addition of the information from the new observation). In the output struct, this is the “a posteriori” state estimate (after the new
measurement information is included).
s.z = observation vector
s.u = input control vector, optional (defaults to zero).
MATRIX VARIABLES:
- s.A = state transition matrix (defaults to identity).
- s.P = covariance of the state vector estimate. In the input struct,
this is “a priori,” and in the output it is “a posteriori.”
(required unless autoinitializing as described below). - s.B = input matrix, optional (defaults to zero).
- s.Q = process noise covariance (defaults to zero).
- s.R = measurement noise covariance (required).
- s.H = observation matrix (defaults to identity).
NORMAL OPERATION:
- define all state definition fields: A,B,H,Q,R
- define intial state estimate: x,P
- obtain observation and control vectors: z,u
- call the filter to obtain updated state estimate: x,P
- return to step (3) and repeat
INITIALIZATION:
If an initial state estimate is unavailable, it can be obtained from the first observation as follows, provided that there are the same number of observable variables as state variables. This “auto- intitialization” is done automatically if s.x is absent or NaN.
x = inv(H)*z P = inv(H)*R*inv(H')
This is mathematically equivalent to setting the initial state estimate
covariance to infinity.
SCALAR EXAMPLE (Automobile Voltimeter):
% Define the system as a constant of 12 volts: clear s s.x = 12; s.A = 1; % Define a process noise (stdev) of 2 volts as the car operates: s.Q = 2^2; % variance, hence stdev^2 % Define the voltimeter to measure the voltage itself: s.H = 1; % Define a measurement error (stdev) of 2 volts: s.R = 2^2; % variance, hence stdev^2 % Do not define any system input (control) functions: s.B = 0; s.u = 0; % Do not specify an initial state: s.x = nan; s.P = nan; % Generate random voltages and watch the filter operate. tru=[]; % truth voltage for t=1:20 tru(end+1) = randn*2+12; s(end).z = tru(end) + randn*2; % create a measurement s(end+1)=kalmanf(s(end)); % perform a Kalman filter iteration end figure hold on grid on % plot measurement data: hz=plot([s(1:end-1).z],'r.'); % plot a-posteriori state estimates: hk=plot([s(2:end).x],'b-'); ht=plot(tru,'g-'); legend([hz hk ht],'observations','Kalman output','true voltage',0) title('Automobile Voltimeter Example') hold off
Download Now
Popularity: 1% [?]
Writing Fast Matlab Code #4: JIT Acceleration
Matlab 6.5 (R13) and later feature the Just-In-Time (JIT) Accelerator for improving the speed of M-functions, particularly with loops. By knowing a few things about the accelerator, you can improve its performance.
The JIT Accelerator is enabled by default. To disable it, type “feature accel off” in the console, and “feature accel on” to enable it again.
As of Matlab R2008b, only a subset of the Matlab language is supported for acceleration. Upon encountering an unsupported feature, acceleration processing falls back to non-accelerated evaluation. Acceleration is most effective when significant contiguous portions of code are supported.
- Data types: Code must use supported data types for acceleration: double (both real and complex), logical, char, int8–32, uint8–32. Some struct, cell, classdef, and function handle usage is supported. Sparse arrays are not accelerated.
- Array shapes: Array shapes of any size with 3 or fewer dimensions are supported. Changing the shape or data type of an array interrupts acceleration. A few limited situations with 4D arrays are accelerated.
- Function calls: Calls to built-in functions and M-functions are accelerated. Calling MEX func- tions and Java interrupts acceleration. (See also page 14 on inlining simple functions.)
- Conditionals and loops: The conditional statements if, elseif, and simple switch statements are supported if the conditional expression evaluates to a scalar. Loops of the form for k=a:b, for k=a:b:c, and while loops are accelerated if all code within the loop is supported.
In-place computation
Introduced in Matlab 7.3 (R2006b), the element-wise operators (+, .*, etc.) and some other functions can be computed in-place. That is, a computation like
x = 5*sqrt(x.ˆ2 + 1);
is handled internally without needing temporary storage for accumulating the result. An M-function can also be computed in-place if its output argument matches one of the input arguments.
x = myfun(x); function x = myfun(x) x = 5*sqrt(x.ˆ2 + 1); return;
To enable in-place computation, the in-place operation must be within an M-function (and for an in- place function, the function itself must be called within an M-function). Currently, there is no support for in-place computation with MEX-functions.
Multithreaded Computation
Matlab 7.4 (R2007a) introduced multithreaded computation for multicore and multiprocessor com- puters. Multithreaded computation accelerates some per-element functions when applied to large arrays (for example, .^, sin, exp) and certain linear algebra functions in the BLAS library. To enable it, select File → Preferences → General → Multithreading and select “Enable multithreaded computation.” Fur- ther control over parallel computation is possible with the Parallel Computing Toolbox using parfor and spmd.
JIT-Accelerated Example
For example, the following loop-heavy code is supported for acceleration:
function B = bilateral(A,sd,sr,R) % The bilateral image denoising filter B = zeros(size(A)); for i = 1:size(A,1) for j = 1:size(A,2) zsum = 0; for m = -R:R if i+m >= 1 && i+m <= size(A,1) for n = -R:R if j+n >= 1 && j+n <= size(A,2) z = exp(-(A(i+m,j+n)-A(i,j))^2/(2*sd^2))... *exp(-(m^2 + n^2)/(2*sr^2)); zsum = zsum + z; B(i,j) = B(i,j) + z*A(i+m,j+n); end end end end end end B(i,j) = B(i,j)/zsum;
For a 128 × 128 input image and R = 3, the run time is 53.3 seconds without acceleration and 0.68 seconds with acceleration.
function B = bilateral(A,sd,sr,R)
% The bilateral image denoising filter
B = zeros(size(A));
for i = 1:size(A,1)
for j = 1:size(A,2)
zsum = 0;
for m = -R:R
if i+m >= 1 && i+m <= size(A,1)
for n = -R:R
if j+n >= 1 && j+n <= size(A,2)
z = exp(-(A(i+m,j+n)-A(i,j))^2/(2*sd^2))…
*exp(-(m^2 + n^2)/(2*sr^2));
zsum = zsum + z;
B(i,j) = B(i,j) + z*A(i+m,j+n);
end
end
end
end
end
end
B(i,j) = B(i,j)/zsum;
Popularity: 1% [?]
Translator
Follow Us
Video Posts
