function y = chirp(t,varargin)
%CHIRP Swept-frequency cosine generator.
%   Y = CHIRP(T) generates samples of a swept-frequency signal at the time 
%   instances defined in array T. By default, the instantaneous frequency
%   at time 0 is 0, and the instantaneous frequency one second later is 100
%   Hz.    
%
%   Y = CHIRP(T,F0,T1,F1) generates samples of a linear swept-frequency
%   signal at the time instances defined in array T. The instantaneous 
%   frequency at time 0 is F0 Hz. The instantaneous frequency F1 is 
%   achieved at time T1.
%
%   Y = CHIRP(T,F0,T1,F1,method) specifies alternate sweep methods.
%   Available methods are 'linear', 'quadratic', and 'logarithmic'; the
%   default is 'linear'. For a logarithmic-sweep, F0 must be greater than 
%   or equal to 1e-6 Hz, which is also the default value.
%
%   Y = CHIRP(T,F0,T1,F1,method,PHI) specifies an initial phase PHI in 
%   degrees. By default, PHI = 0.
%
%   Y = CHIRP(T,FO,T1,F1,'quadratic',PHI,'concave') generates samples of a
%   quadratic swept-frequency chirp signal whose spectrogram is a parabola 
%   with its concavity in the positive frequency axis.
%
%   Y = CHIRP(T,FO,T1,F1,'quadratic',PHI,'convex') generates samples of a
%   quadratic swept-frequency signal whose spectrogram is a parabola with
%   its convexity in the positive frequency axis.
%   
%   Y = CHIRP(...,sigtype) specifies sigtype as 'real' or 'complex'; the 
%   default is 'real'. When sigtype is set to 'real', CHIRP generates real  
%   chirp signals. When set to 'complex', CHIRP generates complex chirp 
%   signals. 
%
%   EXAMPLE 1: Compute the spectrogram of a linear chirp.
%     t = 0:0.001:2;                    % 2 s at 1 kHz sample rate
%     y = chirp(t,0,1,150);             % Start at DC, cross 150 Hz at 
%                                       % t = 1 s
%     spectrogram(y,256,250,256,1E3);   % Display the spectrogram
%
%   EXAMPLE 2: Compute the spectrogram of a quadratic chirp.
%     t = -2:0.001:2;                   % +/-2 s at 1 kHz sample rate
%     y = chirp(t,100,1,200,'q');       % Start at 100 Hz, cross 200 Hz at
%                                       % t = 1 s
%     spectrogram(y,128,120,128,1E3);   % Display the spectrogram
%
%   EXAMPLE 3: Compute the spectrogram of a convex quadratic chirp
%     t = 0:0.001:1;                    % 1 s at 1 kHz sample rate
%     fo = 25;                          % Start at 25 Hz,
%     f1 = 100;                         % go up to 100 Hz
%     y = chirp(t,fo,1,f1,'q',[],'convex');
%     spectrogram(y,256,200,256,1000);  % Display the spectrogram.
%
%   EXAMPLE 4: Compute the spectrogram of a concave quadratic chirp
%     t = 0:0.001:1;                    % 1 s at 1 kHz sample rate
%     fo = 100;                         % Start at 100 Hz,
%     f1 = 25;                          % go down to 25 Hz
%     y = chirp(t,fo,1,f1,'q',[],'concave');
%     spectrogram(y,256,200,256,1000);  % Display the spectrogram.
%
%   EXAMPLE 5: Compute the spectrogram of a logarithmic chirp
%     t = 0:0.001:10;                   % 10 s at 1 kHz sample rate
%     fo = 10;                          % Start at 10 Hz,
%     f1 = 400;                         % go up to 400 Hz
%     y = chirp(t,fo,10,f1,'logarithmic');
%     spectrogram(y,256,200,256,1000);  % Display the spectrogram
%
%   EXAMPLE 6: Compute a complex-valued linear chirp
%     t = 0:0.001:2;                    % 2 s at 1 kHz sample rate
%     fo = -50;                         % Start at -50 Hz,
%     t1 = 1;
%     f1 = 200;                         % cross 200 Hz at t = 1 s
%     y = chirp(t,fo,t1,f1,'complex'); 
%     spectrogram(y,256,200,256,1000,'centered'); % Display the spectrogram
%
%   See also GAUSPULS, SAWTOOTH, SINC, SQUARE.

%   Copyright 1988-2023 The MathWorks, Inc.

%   References:
%   [1] Agilent 33220A 20 MHz Function/Arbitrary Waveform Generator
%       Users Guide, Agilent Technologies, March 2002, pg 298

%#codegen

narginchk(1,8);

validateattributes(t,{'numeric'},{'real','finite'},'chirp','t',1);
t = double(t);

% Parse inputs, and substitute for defaults:
[f0,t1,f1,method,phi,quadtype,isComplex] = parseOptArgs(varargin{:});

% Initializing output variable
y = coder.nullcopy(t);

% Computing the signal according to method
switch method
    case 'polynomial'
        % Polynomial chirp
        y = cos(2*pi*polyval(polyint(f0),t));
        
    case 'linear'
        % Linear chirp
        y = calculateChirp(f0,f1,t1,1,t,phi,isComplex);
        
    case 'quadratic'
        % Compute the quadratic chirp output based on quadtype
        
        % Compute the quadratic chirp (upsweep or downsweep) for
        % complex or concave modes.
        
        % For the default 'concave-upsweep' and 'convex-downsweep' modes
        % call calculateChirp without any changes to the input parameters.
        % For the forced 'convex-upsweep' and 'concave-downsweep' call
        % calculateChirp with f0 and f1 swapped and t = fliplr(-t)
        
        % For 'convex-upsweep' and 'concave-downsweep' modes
        if ((f0(1) < f1) && strcmp(quadtype,'convex') || ((f0(1) > f1) && strcmp(quadtype,'concave')))
            t = fliplr(-t);
            f0temp = f1;
            f1temp = f0(1);
        else
            f0temp = f0(1);
            f1temp = f1;
        end
        y = calculateChirp(f0temp,f1temp,t1,2,t,phi,isComplex);
        
    case 'logarithmic'
        % Logarithmic chirp
        if f0 < 1e-6
            coder.internal.error('siglib:chirp:InvalidF0','F0');
        else
            tempVector = repmat(f1/f0,size(t,1),size(t,2)).^(t./t1);
            instPhi = (t1/log(f1/f0)*f0)*(tempVector-1);
            if (isComplex)
                y = exp(1i*2*pi*phi/360).*complex(cos(2*pi*(instPhi)), -cos(2*pi*(instPhi+90/360)));
            else
                y = cos(2*pi*(instPhi+phi/360));
            end
        end
end
end
%---------------------------------------------------------------------------
function yvalue = calculateChirp(f0,f1,t1,p,t,phi,isComplex)
% General function to compute beta and y for both
% linear and quadratic modes.
% p is the polynomial order (1 for linear and 2 for quadratic)
coder.internal.prefer_const(isComplex);
beta = (f1-f0).*(t1.^(-p));
if (isComplex)
    yvalue = exp(1i*2*pi*phi/360).*complex(cos(2*pi*(beta./(1+p)*(t.^(1+p))+f0*t)),-cos(2*pi*(beta./(1+p)*(t.^(1+p))+f0*t+90/360)));
else
    yvalue = cos(2*pi*(beta./(1+p)*(t.^(1+p))+f0*t+phi/360));
end
end
%---------------------------------------------------------------------------

function [f0,t1,f1,method,phi,quadtype,isComplex] = parseOptArgs(varargin)
% PARSEOPTARGS   Parse optional input arguments.
% Up to 6 optional input arguments,
% f0 - Instantaneous frequency at time 0
% t1 - Reference time
% f1 - Instantaneous frequency at time t1
% method - Sweep method
% phi - Initial phase
% quadtype - Spectrogram shape of quadratic chirp

% Flag to keep track of coder target
isCodegen = ~coder.target('MATLAB');

% Flag to keep track of output signal type
isComplex = false;
TFcomprl = false;

if ( ~ isempty(varargin))
    lastArg = varargin{nargin};
    if ( ~ isnumeric(lastArg)) && (strcmpi(lastArg,'complex'))
        TFrc = strcmpi(lastArg,'complex') || strcmpi(lastArg,'real');
        isComplex = true;
        TFcomprl = true;
        % checks if the output signal type is coder.const()
        coder.internal.errorIf((~coder.internal.isConst(lastArg) && ~TFrc),...
            'siglib:chirp:CoderConstantExpected');

    elseif ( ~ isnumeric(lastArg)) && (strcmpi(lastArg,'real'))
        TFrc = strcmpi(lastArg,'complex') || strcmpi(lastArg,'real');
        isComplex = false;
        TFcomprl = true;
        coder.internal.errorIf((~coder.internal.isConst(lastArg) && ~TFrc),...
            'siglib:chirp:CoderConstantExpected');
    end
end

% Default value depends upon method
f0 = 0;
% Flag to check input status of f0
isF0Set = false;

% Initialize to default values
t1 = 1;
f1 = 100;
phi = 0;

if TFcomprl
    processArgs = nargin-1;
else
    processArgs = nargin;
end

if nargin >= 1 && processArgs >= 1 && ~isempty(varargin{1})   
    f0inp = varargin{1};
    validateattributes(f0inp,{'numeric'},{'vector','finite','real'},'chirp','F0');
    if length(f0inp)>1
        coder.internal.errorIf(isCodegen && length(f0inp) > 1, 'siglib:chirp:PolynomialNotSupported');
        coder.internal.warning('siglib:chirp:syntaxObsolete');
        method = 'polynomial';
    else
        method = 'linear';
    end
    f0 = double(f0inp);
    isF0Set = true;
else
    method = 'linear';
end

if nargin >= 2 && processArgs >= 2 && ~isempty(varargin{2})
    validateattributes(varargin{2},{'numeric'},{'scalar','real','positive','finite'},'chirp','T1');
    t1Inp = varargin{2};
    t1 = double(t1Inp(1));
end

if nargin >= 3 && processArgs >= 3 && ~isempty(varargin{3})
    validateattributes(varargin{3},{'numeric'},{'scalar','real','finite'},'chirp','F1');
    f1Inp = varargin{3};
    f1 = double(f1Inp(1));
end

if nargin >= 4 && processArgs >= 4 && ~isempty(varargin{4})
    methodStr = varargin{4};
    coder.internal.errorIf(~ischar(methodStr) && ~isStringScalar(methodStr), 'siglib:chirp:InvalidInputType','Method','Chirp','character array',class(methodStr));
    if ~strcmp(methodStr, '')
        optionStrings = {'linear','quadratic','logarithmic'};
        method = validatestring(methodStr,optionStrings,'chirp');
    end
end

if nargin >= 5 && processArgs >= 5 && ~isempty(varargin{5})
    validateattributes(varargin{5},{'numeric'},{'scalar','finite','real'},'chirp','PHI');
    phiInp = varargin{5};
    phi = double(phiInp(1));
end

if nargin >= 6 && processArgs >= 6 && ~isempty(varargin{6})
    methodStr = varargin{6};
    coder.internal.errorIf(~ischar(methodStr) && ~isStringScalar(methodStr), 'siglib:chirp:InvalidInputType','Method','Chirp','character array',class(methodStr));
    if ~strcmp(methodStr, '')
        optionStrings = {'concave', 'convex'};
        quadtypeInp = validatestring(methodStr,optionStrings,'chirp');
    else
        quadtypeInp = '';
    end
else
    quadtypeInp = '';
end


% Assigning default value to f0 based on method
if ~isF0Set && strcmp(method,'logarithmic')
    f0 = 1e-6; % Minimum output frequency in Agilent Function generator is 1 microHz
end

% Quadtype is valid only with quadratic method
if strcmp(method,'quadratic') 
    if ~isempty(quadtypeInp)
        quadtype = quadtypeInp;        
    else
        if (f1 > f0)
            quadtype = 'concave'; % Default for upsweep
        else
            quadtype = 'convex'; % Default for downsweep.
        end
    end
else
    coder.internal.errorIf(~isempty(quadtypeInp), ...
        'siglib:chirp:InvalidMode',quadtypeInp, method);
    quadtype = '';
end


end


% [EOF] chirp.m

% LocalWords:  FO fo Agilent quadtype upsweep downsweep PARSEOPTARGS siglib
