Mercurial > forge
changeset 10420:bfc8b33828a5 octave-forge
signal: Adding new functions. clustersegment is old and can be improved a lot.
| author | jpicarbajal |
|---|---|
| date | Sat, 09 Jun 2012 00:52:28 +0000 |
| parents | 66ec9da43f53 |
| children | 4b25237504f0 |
| files | main/signal/INDEX main/signal/NEWS main/signal/inst/clustersegment.m main/signal/inst/movingrms.m main/signal/inst/schtrig.m |
| diffstat | 5 files changed, 268 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/main/signal/INDEX Fri Jun 08 19:20:19 2012 +0000 +++ b/main/signal/INDEX Sat Jun 09 00:52:28 2012 +0000 @@ -27,6 +27,7 @@ sgolayfilt sosfilt medfilt1 + movingrms Filter analysis freqs freqs_plot grpdelay @@ -119,7 +120,8 @@ kaiser nuttallwin triang - gaussian gausswin + gaussian + gausswin tukeywin rectwin welchwin @@ -160,3 +162,5 @@ wrev zerocrossing sampled2continuous + schtrig + clustersegment
--- a/main/signal/NEWS Fri Jun 08 19:20:19 2012 +0000 +++ b/main/signal/NEWS Sat Jun 09 00:52:28 2012 +0000 @@ -4,6 +4,9 @@ signal-1.1.4 Release Date: 2012-XX-XX Release Manager: =============================================================================== + ** The following functions are new: + movingrms schtrig clustersegment + ** signal is no longer dependent on the image package. ** signal is now dependent on the general package.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/clustersegment.m Sat Jun 09 00:52:28 2012 +0000 @@ -0,0 +1,69 @@ +## Copyright (c) 2010 Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, see <http://www.gnu.org/licenses/>. + +%% contRange = clustersegment(xhi) +% The function calculates the initial and end index values of the sequences of +% 1's in the rows of xhi. The result is returned in a cell of size 1xNp, with Np +% the numer of rows in xhi. Each element of the cell has two rows; the first row +% is the inital index of a sequence of 1's and the second row is the end index +% of that sequence. + +function contRange = clustersegment(xhi) + + % Find discontinuities + bool_discon = diff(xhi,1,2); + [Np Na] = size(xhi); + contRange = cell(1,Np); + + for i = 1:Np + idxUp = find(bool_discon(i,:)>0)+1; + idxDwn = find(bool_discon(i,:)<0); + tLen = length(idxUp) + length(idxDwn); + + if xhi(i,1)==1 + % first event was down + contRange{i}(1) = 1; + contRange{i}(2:2:tLen+1) = idxDwn; + contRange{i}(3:2:tLen+1) = idxUp; + else + % first event was up + contRange{i}(1:2:tLen) = idxUp; + contRange{i}(2:2:tLen) = idxDwn; + end + + if xhi(i,end)==1 + % last event was up + contRange{i}(end+1) = Na; + end + + tLen = length(contRange{i}); + if tLen ~=0 + contRange{i}=reshape(contRange{i},2,tLen/2); + end + + end + +endfunction + +%!demo +%! xhi = [0 0 1 1 1 0 0 1 0 0 0 1 1]; +%! ranges = intervalSegment(xhi) +%! +%! % The first sequence of 1's in xhi is +%! xhi(ranges{1}(1,:)) + +%!demo +%! xhi = rand(3,10)>0.4 +%! ranges = intervalSegment(xhi)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/movingrms.m Sat Jun 09 00:52:28 2012 +0000 @@ -0,0 +1,96 @@ +## Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{rmsx},@var{w}] =} movingrms (@var{x},@var{w},@var{rc},@var{Fs}=1) +## Calculates moving RMS value of the signal in @var{x}. +## +## The signla is convoluted against a sigmoid window of width @var{w} and risetime +## @var{rc}. The units of these to parameters are relative ot the vlaue of the sampling +## frequency given in @var{Fs} (Default value = 1). +## +## Run @code{demo movingrms} to see an example. +## +## @seealso{sigmoid_train} +## @end deftypefn + +function [rmsx w]= movingrms (x,width, risetime, Fs=1) + + [N nc] = size (x); + if width*Fs > N/2 + idx = [1 N]; + w = ones(N,1); + else + idx = round ((N + width*Fs*[-1 1])/2); + w = sigmoid_train ((1:N)', idx, risetime*Fs); + end + fw = fft (w.^2); + fx = fft (x.^2); + rmsx = real(ifft (fx.*fw)/(N-1)); + + rmsx (rmsx < eps*max(rmsx(:))) = 0; + + rmsx = circshift (sqrt (rmsx), round(mean(idx))); +% w = circshift (w, -idx(1)); + +endfunction + +%!demo +%! N = 128; +%! t = linspace(0,1,N)'; +%! x = sigmoid_train (t,[0.4 inf],1e-2).*(2*rand(size(t))-1); +%! +%! Fs = 1/diff(t(1:2)); +%! width = 0.05; +%! rc = 5e-3; +%! [wx w] = movingrms (zscore (x),width,rc,Fs); +%! +%! close all +%! figure () +%! +%! area (t,wx,'facecolor',[0.85 0.85 1],'edgecolor','b','linewidth',2); +%! hold on; +%! h = plot (t,x,'r-;Data;',t,w,'g-;Window;'); +%! set (h, 'linewidth', 2); +%! hold off +%! +%! % --------------------------------------------------------------------------- +%! % The shaded plot shows the local RMS of the Data: white noise with onset at +%! % aprox. t== 0.4. +%! % The observation window is also shown. + +%!demo +%! N = 128; +%! t = linspace(0,1,N)'; +%! x = exp(-((t-0.5)/0.1).^2) + 0.1*rand(N,1); +%! +%! Fs = 1/diff(t(1:2)); +%! width = 0.1; +%! rc = 2e-3; +%! [wx w] = movingrms (zscore (x),width,rc,Fs); +%! +%! close all +%! figure () +%! +%! area (t,wx,'facecolor',[0.85 0.85 1],'edgecolor','b','linewidth',2); +%! hold on; +%! h = plot (t,x,'r-;Data;',t,w,'g-;Window;'); +%! set (h, 'linewidth', 2); +%! hold off +%! +%! % --------------------------------------------------------------------------- +%! % The shaded plot shows the local RMS of the Data: Gausian with centered at +%! % aprox. t== 0.5. +%! % The observation window is also shown.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/schtrig.m Sat Jun 09 00:52:28 2012 +0000 @@ -0,0 +1,95 @@ +## Copyright (c) 2012 Juan Pablo Carbajal <carbajal@ifi.uzh.ch> +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{rmsx},@var{w}] =} schtrig (@var{x},@var{lvl},@var{rst}=1) +## Implements a multisignal Schmitt trigger with levels @var{lvl}. +## +## The triger works along the first dimension of the array @var{x}. When @code{@var{rst}==1} +## the state of the trigger for all signals is set to the low state (i.e. 0). +## +## Run @code{demo schtrig} to see an example. +## +## @seealso{clustersegment} +## @end deftypefn + +function v = schtrig (x, lvl, rst = 1) + + persistent st0; + + if length(lvl) == 1 + lvl = abs (lvl)*[1 -1]; + else + lvl = sort(lvl,'descend'); + end + + [nT nc] = size(x); + + v = NA (nT, nc); + + if rst || isempty(st0) + st0 = zeros(1,nc); + printf ("Trigger initialized!\n"); + flush (stdout); + end + + v(1,:) = st0; + + % Signal is above up level + up = x > lvl(1); + v(up) = 1; + + % Signal is below down level + dw = x < lvl(2); + v(dw) = 0; + + % Resolve intermediate states + % Find data between the levels + idx = isnan (v); + ranges = clustersegment (idx'); + + for i=1:nc + % Record the state at the begining of the interval between levels + if !isempty (ranges{i}) + prev = ranges{i}(1,:)-1; + prev(prev<1) = 1; + st0 = v(prev,i); + + % Copy the initial state to the interval + ini_idx = ranges{i}(1,:); + end_idx = ranges{i}(2,:); + for j =1:length(ini_idx) + v(ini_idx(j):end_idx(j),i) = st0(j); + end + end + end + + st0 = v(end,:); + +endfunction + +%!demo +%! t = linspace(0,1,100)'; +%! x = sin (2*pi*2*t) + sin (2*pi*5*t).*[0.8 0.3]; +%! +%! lvl = [0.8 0.25]'; +%! v = schtrig (x,lvl); +%! +%! h = plot(t,x,t,v); +%! set (h([1 3]),'color','b'); +%! set (h([2 4]),'color',[0 1 0.5]); +%! set (h,'linewidth',2); +%! line([0; 1],lvl([1; 1]),'color','r'); +%! line([0;1],lvl([2;2]),'color','b')