function [ cc, ss, NCOstate ] = NCO( Ns, fs, NCOstate, multPhase )
% NCO - numerically controlled oscillator
% 
% The numerically controlled oscillator produces a block of Ns samples of a
% cosine and sine signal. The frequency and phase is determined by the
% state maintained in NCOstate. 
% 
% The phase of the cosine and sine signals is continuous between two
% invocations of NCO.
%
% NCOstate is a structure with two fields. 
%   'freq' - frequency of the sinusoids produced by the NCO
%   'firstPhase' - phase angle of the first sample produced by the NCO in a block 
%
% Usage:
%   [ cc, ss, NCOstate ] = NCO( Ns, fs, NCOstate[, multPhase] )

if nargin < 4 || isempty(multPhase)
    multPhase = [1,1];
end

%% figure out the phase function Psi
% first the slope
kPsi = 2*pi*NCOstate.freq/fs;

% then add initial phase
Psi = NCOstate.firstPhase + (0:Ns-1)*kPsi;

%% compute cosine and sine
cc = cos(multPhase(1)*Psi);
ss = sin(multPhase(2)*Psi);

%% update state for next invocation
NCOstate.firstPhase = mod(Psi(end)+kPsi, 2*pi);
end

