X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=octave_packages%2Fsignal-1.1.3%2Fpulstran.m;fp=octave_packages%2Fsignal-1.1.3%2Fpulstran.m;h=bd0426caabad5d71fa0c0bc979b6a8004c1a5577;hb=f5f7a74bd8a4900f0b797da6783be80e11a68d86;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/signal-1.1.3/pulstran.m b/octave_packages/signal-1.1.3/pulstran.m new file mode 100644 index 0000000..bd0426c --- /dev/null +++ b/octave_packages/signal-1.1.3/pulstran.m @@ -0,0 +1,144 @@ +## Copyright (C) 2000 Paul Kienzle +## +## 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 . + +## usage: y=pulstran(t,d,'func',...) +## y=pulstran(t,d,p,Fs,'interp') +## +## Generate the signal y=sum(func(t+d,...)) for each d. If d is a +## matrix of two columns, the first column is the delay d and the second +## column is the amplitude a, and y=sum(a*func(t+d)) for each d,a. +## Clearly, func must be a function which accepts a vector of times. +## Any extra arguments needed for the function must be tagged on the end. +## +## Example +## fs = 11025; # arbitrary sample rate +## f0 = 100; # pulse train sample rate +## w = 0.001; # pulse width of 1 millisecond +## auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, 'rectpuls', w), fs); +## +## If instead of a function name you supply a pulse shape sampled at +## frequency Fs (default 1 Hz), an interpolated version of the pulse +## is added at each delay d. The interpolation stays within the the +## time range of the delayed pulse. The interpolation method defaults +## to linear, but it can be any interpolation method accepted by the +## function interp1. +## +## Example +## fs = 11025; # arbitrary sample rate +## f0 = 100; # pulse train sample rate +## w = boxcar(10); # pulse width of 1 millisecond at 10 kHz +## auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, w, 10000), fs); + +## TODO: Make it faster. It is currently unusable for anything real. +## TODO: It may not be possible to speed it up with the present interface. +## TODO: See speech/voice.m for a better way. + +## Note that pulstran can be used for some pretty strange things such +## as simple band-limited interpolation: +## xf = 0:0.05:10; yf = sin(2*pi*xf/5); +## xp = 0:10; yp = sin(2*pi*xp/5); # .2 Hz sine sampled every second +## s = pulstran(xf, [xp, yp],'sinc'); +## plot(f, yf, ";original;", xf, s, ";sinc;",xp,yp,"*;;"); +## You wouldn't want to do this in practice since it is expensive, and +## since it works much better with a windowed sinc function, at least +## for short samples. + +function y = pulstran(t, d, pulse, varargin) + + if nargin<3 || (!ischar(pulse) && nargin>5) + print_usage; + endif + y = zeros(size(t)); + if isempty(y), return; endif + if rows(d) == 1, d=d'; endif + if columns(d) == 2, + a=d(:,2); + else + a=ones(rows(d),1); + endif + if ischar(pulse) + ## apply function t+d for all d + for i=1:rows(d) + y = y+a(i)*feval(pulse,t-d(i,1),varargin{:}); + endfor + else + ## interpolate each pulse at the specified times + Fs = 1; method = 'linear'; + if nargin==4 + arg=varargin{1}; + if ischar(arg), + method=arg; + else + Fs = arg; + endif + elseif nargin==5 + Fs = varargin{1}; + method = varargin{2}; + endif + span = (length(pulse)-1)/Fs; + t_pulse = (0:length(pulse)-1)/Fs; + for i=1:rows(d) + dt = t-d(i,1); + idx = find(dt>=0 & dt<=span); + y(idx) = y(idx) + a(i)*interp1(t_pulse, pulse, dt(idx), method); + endfor + endif +endfunction + +%!error pulstran +%!error pulstran(1,2,3,4,5,6) + +%!## parameter size and shape checking +%!shared t,d +%! t = 0:0.01:1; d=0:0.1:1; +%!assert (isempty(pulstran([], d, 'sin'))); +%!assert (pulstran(t, [], 'sin'), zeros(size(t))); +%!assert (isempty(pulstran([], d, boxcar(5)))); +%!assert (pulstran(t, [], boxcar(5)), zeros(size(t))); +%!assert (size(pulstran(t,d,'sin')), size(t)); +%!assert (size(pulstran(t,d','sin')), size(t)); +%!assert (size(pulstran(t',d,'sin')), size(t')); +%!assert (size(pulstran(t,d','sin')), size(t)); + +%!demo +%! fs = 11025; # arbitrary sample rate +%! f0 = 100; # pulse train sample rate +%! w = 0.003; # pulse width of 3 milliseconds +%! t = 0:1/fs:0.1; d=0:1/f0:0.1; # define sample times and pulse times +%! a = hanning(length(d)); # define pulse amplitudes +%! +%! subplot(221); title("rectpuls"); +%! auplot(pulstran(t', d', 'rectpuls', w), fs); +%! hold on; plot(d*1000,ones(size(d)),'g*;pulse;'); hold off; +%! +%! subplot(223); title("sinc => band limited interpolation"); +%! auplot(pulstran(f0*t, [f0*d', a], 'sinc'), fs); +%! hold on; plot(d*1000,a,'g*;pulse;'); hold off; +%! +%! subplot(222); title("interpolated boxcar"); +%! pulse = boxcar(30); # pulse width of 3 ms at 10 kHz +%! auplot(pulstran(t, d', pulse, 10000), fs); +%! hold on; plot(d*1000,ones(size(d)),'g*;pulse;'); hold off; +%! +%! subplot(224); title("interpolated asymmetric sin"); +%! pulse = sin(2*pi*[0:0.0001:w]/w).*[w:-0.0001:0]; +%! auplot(pulstran(t', [d', a], pulse', 10000), fs); +%! hold on; plot(d*1000,a*w,'g*;pulse;'); hold off; title(""); +%! +%! %---------------------------------------------------------- +%! % Should see (1) rectangular pulses centered on *, +%! % (2) rectangular pulses to the right of *, +%! % (3) smooth interpolation between the *'s, and +%! % (4) asymetric sines to the right of *