view signals/@sigresample/construct.m @ 1:289445d368a7

import.

function s2=construct(sig)

	fprintf('\nResampling %s to %g Hz...\n',tostring(sig.source),sig.rate);
	f1=rate(sig.source);
	f2=sig.rate;
	c=channels(sig);

	[p,q]=rat(f2/f1,1e-12);
	sysobj=design(sig,p,q);


	delay=(length(sysobj.Numerator)-1)/2; 
	outdelay = ceil(delay/q);
	indelay  = ceil(delay/p);
	m=max(1,floor(sig.opts.bs/q)); % try to approximate requested input block size
	fprintf('Input/output delays are %d and %d.\n',indelay,outdelay);
	fprintf('Read block size is %d.\n\n',m*q);

	s1=construct(sig.source & sigarray(zeros(channels(sig.source),indelay)));
	r1=s1.reader(m*q); % one and only input reader
	chunk=uint32(m*p); % filter output block size
	CHUNK=1:chunk;
	queue=[];

	s2.start   = s1.start;
	s2.stop    = s1.stop;
	s2.dispose = @dispose;
	s2.reader  = @reader;

	% drop some samples to account for filter delay
	s2.start();
	sigreadn(s2,outdelay);
	s2.stop();

	function dispose, s1.dispose(); release(sysobj); end 
	function r2=reader(n)
		outbuf=zeros(c,n);
		r2=@next;

		function [x,rem]=next 
			% transfer up to n queued samples to outbuf
			pos=uint32(size(queue,2));
			if pos>=n % enough samples already waiting
				outbuf=queue(:,1:n);
				queue=queue(:,n+1:end);
				rem=0;
			else
				if pos==0, toread=n;
				else % use up queue
					outbuf(:,1:pos)=queue; queue=[];
					toread=n-pos;
				end

				% transfer complete chunks
				while toread>=chunk
					[outbuf(:,pos+CHUNK),rem]=filter_next;
					toread=toread-chunk;
					pos=pos+chunk;
					if rem>0 % we ran out of samples
						rem=rem+toread; % account for rest of missing samples (not just this chunk) 
						toread=0; % causes immediate exit after this
					end
				end

				% transfer partial chunk if necessary
				if toread>0
					[y,rem]=filter_next;
					outbuf(:,pos+1:n)=y(:,1:toread);
					queue=y(:,toread+1:chunk-rem);
					rem=max(0,toread-(chunk-rem));
				end
			end
			x=outbuf;
		end

		% returns the next block of m*p output samples 
		function [x,rem]=filter_next
			[y,rem1]=r1();
			if rem1>0,
				y(:,end-rem1+1:end)=zeros(c,rem1); % pad with zeros
				rem=uint32(ceil(rem1*p/q)); % round down the number of valid samples
			else
				rem=rem1;
			end
			x=step(sysobj,y')';
		end
	end
end