m8ta
You are not authenticated, login. |
|
{496} | ||
| ||
{585} | ||
LMS-based adaptive decorrelator, xn is the noise, xs is the signal, len is the length of the signal, delay is the delay beyond which the autocorrelation function of the signal is zero but the acf of the noise is non-zero. The filter is very simple, and should be easy to implement in a DSP. function [y,e,h] = lms_test(xn, xs, len, delay) h = zeros(len, 1); x = xn + xs; for k = 1:length(x)-len-delay y(k) = x(k+delay:k+len-1+delay) * h ; e(k) = x(k) - y(k); h = h + 0.0004 * e(k) * x(k+delay:k+len-1+delay)'; endIt works well if the noise source is predictable & stable: (black = sinusoidal noise, red = output, green = error in output) Now, what if the amplitude of the corrupting sinusoid changes (e.g. due to varying electrode properties during movement), and the changes per cycle are larger than the amplitude of the signal? The signal will be swamped! The solution to this is to adapt the decorrelating filter slowly, by adding an extra (multiplicative, nonlinear) gain term to track the error in terms of the absolute values of the signals (another nonlinearity). So, if the input signal is on average larger than the output, the gain goes up and vice-versa. See the code. function [y,e,h,g] = lms_test(xn, xs, len, delay) h = zeros(len, 1); x = xn + xs; gain = 1; e = zeros(size(x)); e2 = zeros(size(x)); for k = 1:length(x)-len-delay y(k) = x(k+delay:k+len-1+delay) * h; e(k) = (x(k) - y(k)); h = h + 0.0002 * e(k) * x(k+delay:k+len-1+delay)'; % slow adaptation. y2(k) = y(k) * gain; e2(k) = abs(x(k)) - abs(y2(k)); gain = gain + 1 * e2(k) ; gain = abs(gain); if (gain > 3) gain = 3; end g(k) = gain; end If, like me, you are interested in only the abstract features of the signal, and not an accurate reconstruction of the waveform, then the gain signal (g above) reflects the signal in question (once the predictive filter has adapted). In my experiments with a length 16 filter delayed 16 samples, extracting the gain signal and filtering out out-of-band information yielded about +45db improvement in SNR. This was with a signal 1/100th the size of the disturbing amplitude-modulated noise. This is about twice as good as the human ear/auditory system in my tests.
It doesn't look like much, but it is just perfect for EMG signals corrupted by time-varying 60hz noise. | ||
{519} | ||
devices:
devices that can be turned off & on to save power (e.g. actually disconnected from power through a P channel MOSFET. must be careful to tristate all outputs before disabling, otherwise we'll get current through the ESD protection diodes )
| ||
{760} |
ref: -0
tags: LDA myopen linear discriminant analysis classification
date: 01-03-2012 02:36 gmt
revision:2
[1] [0] [head]
|
|
How does LDA (Linear discriminant analysis) work? It works by projecting data points onto a series of planes, one per class of output, and then deciding based which projection plane is the largest. Below, to the left is a top-view of this projection with 9 different classes of 2D data each in a different color. Right is a size 3D view of the projection - note the surfaces seem to form a parabola. Here is the matlab code that computes the LDA (from myopen's ceven % TrainData and TrainClass are inputs, column major here. % (observations on columns) N = size(TrainData,1); Ptrain = size(TrainData,2); Ptest = size(TestData,2); % add a bit of interpolating noise to the data. sc = std(TrainData(:)); TrainData = TrainData + sc./1000.*randn(size(TrainData)); K = max(TrainClass); % number of classes. %%-- Compute the means and the pooled covariance matrix --%% C = zeros(N,N); for l = 1:K; idx = find(TrainClass==l); % measure the mean per class Mi(:,l) = mean(TrainData(:,idx)')'; % sum all covariance matrices per class C = C + cov((TrainData(:,idx)-Mi(:,l)*ones(1,length(idx)))'); end C = C./K; % turn sum into average covariance matrix Pphi = 1/K; Cinv = inv(C); %%-- Compute the LDA weights --%% for i = 1:K Wg(:,i) = Cinv*Mi(:,i); % this is the slope of the plane Cg(:,i) = -1/2*Mi(:,i)'*Cinv*Mi(:,i) + log(Pphi)'; % and this, the origin-intersect. end %%-- Compute the decision functions --%% Atr = TrainData'*Wg + ones(Ptrain,1)*Cg; % see - just a simple linear function! Ate = TestData'*Wg + ones(Ptest,1)*Cg; errtr = 0; AAtr = compet(Atr'); % this compet function returns a sparse matrix with a 1 % in the position of the largest element per row. % convert to indices with vec2ind, below. TrainPredict = vec2ind(AAtr); errtr = errtr + sum(sum(abs(AAtr-ind2vec(TrainClass))))/2; netr = errtr/Ptrain; PeTrain = 1-netr; | ||
{594} | ||
{586} | ||
Myopen amplifiers & analog/digital filters & NLMS are working properly! Below, a recording from my deltiod as I held my arm up: (only one EMG channel active, ground was my knee)) Yellow traces are raw inputs from ADC, blue are the output from the IIR / adaptive filters; hence, you only see 8 of the 16 channels. Read from bottom to top (need a -1 in some opengl matrix somewhere...) Below, the system with no input except for free wires attached to one channel (and picking up ambient noise). For this channel, NLMS could not remove the square wave - too many harmonics - but for all other channels the algorthim properly removes 60hz interference :) Now, let me clean this EEG paste off my shoulder & leg ;) |