Dramatically faster Matlab with QLA (beta)
Filed under: Algorithms, Linear Algebra, Numerical Analysis
We recently released the beta version of QLA: The Quick Linear Algebra Library.
QLA provides fast linear algebra in MATLAB, requiring only a regular CPU. For a minuscule sacrifice in accuracy, you get a mammoth increase in speed.
* How fast? Speedups over built-in MATLAB functions range from 10x to 1,000x+ on benchmarks.
* How easy? Just add a ‘q’ in front of your normal function calls – e.g., ‘qsvd(A)’ instead of ’svd(A)’.
* Includes SVD, linear systems, least squares, PCA, and more.
QLA is currently in free public beta. Download and find more details here:
http://massiveanalytics.com
You can also join the QLA discussion here:
http://www.accelereyes.com/forums/viewforum.php?f=15
Enjoy!
Michael Holmes / Massive Analytics
John Melonakos / AccelerEyes
Popularity: 1% [?]
GUI Examples #11: Explore slider and editbox interaction
An editbox is used to display and manipulate the current position of a slider.
function [] = GUI_11() % Demonstrate how to display & change a slider's position with an edit box. % Slide the slider and it's position will be shown in the editbox. % Enter a valid number in the editbox and the slider will be moved to that % position. If the number entered is outside the range of the slider, % the number will be reset. % % % Author: Matt Fig % Date: 7/15/2009 S.fh = figure('units','pixels',... 'position',[300 300 300 100],... 'menubar','none',... 'name','GUI_11',... 'numbertitle','off',... 'resize','off'); S.sl = uicontrol('style','slide',... 'unit','pix',... 'position',[20 10 260 30],... 'min',1,'max',100,'val',50); S.ed = uicontrol('style','edit',... 'unit','pix',... 'position',[20 50 260 30],... 'fontsize',16,... 'string','50'); set([S.ed,S.sl],'call',{@ed_call,S}); % Shared Callback. function [] = ed_call(varargin) % Callback for the edit box and slider. [h,S] = varargin{[1,3]}; % Get calling handle and structure. switch h % Who called? case S.ed L = get(S.sl,{'min','max','value'}); % Get the slider's info. E = str2double(get(h,'string')); % Numerical edit string. if E >= L{1} && E <= L{2} set(S.sl,'value',E) % E falls within range of slider. else set(h,'string',L{3}) % User tried to set slider out of range. end case S.sl set(S.ed,'string',get(h,'value')) % Set edit to current slider. otherwise % Do nothing, or whatever. end
Popularity: 1% [?]
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% [?]
Translator
Follow Us
Video Posts
