Writing Fast Matlab Code #10:Signal Processing

December 15, 2009 by Pascal Getreuer · Leave a Comment
Filed under: Code Optimization, Signal Processing, Tutorials

Rating: 0 (from 0 votes)

Rating: 0.0/10 (0 votes cast)

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 2^K.

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].

Rating: 0.0/10 (0 votes cast)

Rating: 0 (from 0 votes)

Popularity: 1% [?]

Share and Enjoy:

FFW Fastest Filtering in the West

August 10, 2009 by Luigi Rosa · 2 Comments
Filed under: Signal Processing

Rating: 0 (from 0 votes)

Rating: 0.0/10 (0 votes cast)

The FFW package is an FFT-based algorithm for a fast 2D convolution using the overlap-add method. The overlap-add method is based on the fundamental technique in DSP: decompose the signal into simple components, process each of the components in some useful way, and recombine the processed components into the final signal. This is possible since the convolutional operator is linear. The FFW package works similarly to fftfilt function (Matlab Image Processing Toolbox) but in a deeper way: all possible lengths for vectors are considered and not only lengths which are powers of two. This is highly necessary since the FFTW package ( for more details visit http://www.fftw.org ) includes codelets optimized also for other fixed sizes. Codelets are produced automatically by the FFTW codelet generator: you can add your own codelets and re-calculate the execution times for each FFT. The execution times for:

FFT of real 1D vectors
FFT of complex 1D vectors
IFFT of complex 1D vectors

have been calculated with the script provatempo2.m from length N = 3 up to length N = 2048. A finer determination of such times can be done using PAPI for Matlab ( available here or at http://icl.cs.utk.edu/papi ). A 2D FFT (see Matlab command fft2) is decomposed into several 1D FFTs: the FFT operator for an N-dimensional array can in fact be splitted into several 1-dimensional FFTs of monodimensional arrays. The FFW algorithm automatically selects which is the best choice (first dimension, second dimension and best lengths for overlap-add method) and calculates the 2D convolution.

The FFW package can be easily used to improve speed performances of:

2D convolution (Matlab function conv2)
2D filtering (Matlab function filter2)
2D cross-correlation (Matlab function xcorr2)
Normalized cross-correlation (Matlab function normxcorr2)

Index Terms: Matlab, source, code, conv, conv2, fft, fft2, convolution, correlation, normalized cross-correlation, filtering, filter, filters, fast fourier transform, kernel.

For more details please visit
http://www.advancedsourcecode.com/ffw.asp

Luigi Rosa
Via Centrale 35
67042 Civita Di Bagno
L’Aquila ITALY
mobile +39 3207214179
web site http://www.advancedsourcecode.com