function [y,phi]=block(o,phi,cutoff,f)
% block_sr - Generate block of blit signal data
%
% block_sr :: 
%    blit       ~'blit object'
%    0--1       ~'initial phase (in cycles)'
%    N:natural  ~'number of samples to compute',
%    0--1       ~'normalised frequency (1=sampling freq)',
%    0--0.5     ~'cutoff harmonics above this frequency [0.5]',
% -> [[1,N]]    ~'band-limited impulse train',
%    0--1       ~'initial phase for next block'.

m = max(1,2*floor(max(cutoff./f))-1); % number of harmonics to keep
u = cumsum([phi,f],2);
y = m*f.*diric(2*pi*u(1:end-1),m)-f;
phi = mod(u(end),1);

% If diric is not available then use this
% y=zeros(1,T);
% phase = pi*(phi+f*(0:T-1));
% modcheck=mod(phase,pi);
% I=modcheck>0.001 & modcheck<(pi-0.001);
% y(~I) = f*m;
% y(I)  = f*sin(m*phase(I))./sin(phase(I));


