octave_packages/signal-1.1.3/pulstran.m

   1 ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## usage: y=pulstran(t,d,'func',...)
  17 ##        y=pulstran(t,d,p,Fs,'interp')
  18 ##
  19 ## Generate the signal y=sum(func(t+d,...)) for each d.  If d is a
  20 ## matrix of two columns, the first column is the delay d and the second
  21 ## column is the amplitude a, and y=sum(a*func(t+d)) for each d,a.
  22 ## Clearly, func must be a function which accepts a vector of times.
  23 ## Any extra arguments needed for the function must be tagged on the end.
  24 ##
  25 ## Example
  26 ##   fs = 11025;  # arbitrary sample rate
  27 ##   f0 = 100;    # pulse train sample rate
  28 ##   w = 0.001;   # pulse width of 1 millisecond
  29 ##   auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, 'rectpuls', w), fs);
  30 ##
  31 ## If instead of a function name you supply a pulse shape sampled at
  32 ## frequency Fs (default 1 Hz),  an interpolated version of the pulse
  33 ## is added at each delay d.  The interpolation stays within the the
  34 ## time range of the delayed pulse.  The interpolation method defaults
  35 ## to linear, but it can be any interpolation method accepted by the
  36 ## function interp1.
  37 ##
  38 ## Example
  39 ##   fs = 11025;  # arbitrary sample rate
  40 ##   f0 = 100;    # pulse train sample rate
  41 ##   w = boxcar(10);  # pulse width of 1 millisecond at 10 kHz
  42 ##   auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, w, 10000), fs);
  43
  44 ## TODO: Make it faster.  It is currently unusable for anything real.
  45 ## TODO: It may not be possible to speed it up with the present interface.
  46 ## TODO: See speech/voice.m for a better way.
  47
  48 ## Note that pulstran can be used for some pretty strange things such
  49 ## as simple band-limited interpolation:
  50 ##     xf = 0:0.05:10; yf = sin(2*pi*xf/5);
  51 ##     xp = 0:10; yp = sin(2*pi*xp/5); # .2 Hz sine sampled every second
  52 ##     s = pulstran(xf, [xp, yp],'sinc');
  53 ##     plot(f, yf, ";original;", xf, s, ";sinc;",xp,yp,"*;;");
  54 ## You wouldn't want to do this in practice since it is expensive, and
  55 ## since it works much better with a windowed sinc function, at least
  56 ## for short samples.
  57
  58 function y = pulstran(t, d, pulse, varargin)
  59
  60   if nargin<3 || (!ischar(pulse) && nargin>5)
  61     print_usage;
  62   endif
  63   y = zeros(size(t));
  64   if isempty(y), return; endif
  65   if rows(d) == 1, d=d'; endif
  66   if columns(d) == 2,
  67     a=d(:,2);
  68   else
  69     a=ones(rows(d),1);
  70   endif
  71   if ischar(pulse)
  72     ## apply function t+d for all d
  73     for i=1:rows(d)
  74       y = y+a(i)*feval(pulse,t-d(i,1),varargin{:});
  75     endfor
  76   else
  77     ## interpolate each pulse at the specified times
  78     Fs = 1; method = 'linear';
  79     if nargin==4
  80       arg=varargin{1};
  81       if ischar(arg),
  82         method=arg;
  83       else
  84         Fs = arg;
  85       endif
  86     elseif nargin==5
  87       Fs = varargin{1};
  88       method = varargin{2};
  89     endif
  90     span = (length(pulse)-1)/Fs;
  91     t_pulse = (0:length(pulse)-1)/Fs;
  92     for i=1:rows(d)
  93       dt = t-d(i,1);
  94       idx = find(dt>=0 & dt<=span);
  95       y(idx) = y(idx) + a(i)*interp1(t_pulse, pulse, dt(idx), method);
  96     endfor
  97   endif
  98 endfunction
  99
 100 %!error pulstran
 101 %!error pulstran(1,2,3,4,5,6)
 102
 103 %!## parameter size and shape checking
 104 %!shared t,d
 105 %! t = 0:0.01:1; d=0:0.1:1;
 106 %!assert (isempty(pulstran([], d, 'sin')));
 107 %!assert (pulstran(t, [], 'sin'), zeros(size(t)));
 108 %!assert (isempty(pulstran([], d, boxcar(5))));
 109 %!assert (pulstran(t, [], boxcar(5)), zeros(size(t)));
 110 %!assert (size(pulstran(t,d,'sin')), size(t));
 111 %!assert (size(pulstran(t,d','sin')), size(t));
 112 %!assert (size(pulstran(t',d,'sin')), size(t'));
 113 %!assert (size(pulstran(t,d','sin')), size(t));
 114
 115 %!demo
 116 %! fs = 11025;                   # arbitrary sample rate
 117 %! f0 = 100;                     # pulse train sample rate
 118 %! w = 0.003;                    # pulse width of 3 milliseconds
 119 %! t = 0:1/fs:0.1; d=0:1/f0:0.1; # define sample times and pulse times
 120 %! a = hanning(length(d));       # define pulse amplitudes
 121 %!
 122 %! subplot(221); title("rectpuls");
 123 %! auplot(pulstran(t', d', 'rectpuls', w), fs);
 124 %! hold on; plot(d*1000,ones(size(d)),'g*;pulse;'); hold off;
 125 %!
 126 %! subplot(223); title("sinc => band limited interpolation");
 127 %! auplot(pulstran(f0*t, [f0*d', a], 'sinc'), fs);
 128 %! hold on; plot(d*1000,a,'g*;pulse;'); hold off;
 129 %!
 130 %! subplot(222); title("interpolated boxcar");
 131 %! pulse = boxcar(30);  # pulse width of 3 ms at 10 kHz
 132 %! auplot(pulstran(t, d', pulse, 10000), fs);
 133 %! hold on; plot(d*1000,ones(size(d)),'g*;pulse;'); hold off;
 134 %!
 135 %! subplot(224); title("interpolated asymmetric sin");
 136 %! pulse = sin(2*pi*[0:0.0001:w]/w).*[w:-0.0001:0];
 137 %! auplot(pulstran(t', [d', a], pulse', 10000), fs);
 138 %! hold on; plot(d*1000,a*w,'g*;pulse;'); hold off; title("");
 139 %!
 140 %! %----------------------------------------------------------
 141 %! % Should see (1) rectangular pulses centered on *,
 142 %! %            (2) rectangular pulses to the right of *,
 143 %! %            (3) smooth interpolation between the *'s, and
 144 %! %            (4) asymetric sines to the right of *