function w = psdfreqvec(varargin)
%This undocumented function may be removed in a future release.
%
%PSDFREQVEC Frequency vector
%   PSDFREQVEC('Npts',NPTS) returns a frequency vector in radians based on
%   the number of points specified in NPTS. The vector returned assumes 2pi
%   periodicity.
%
%   PSDFREQVEC('Fs',FS) specifies the sampling frequency FS in hertz. By
%   default Fs is set to empty indicating normalized frequency.
%   
%   PSDFREQVEC('CenterDC',CENTERDC) specifies a boolean value in CENTERDC
%   which indicates if zero hertz should be in the center of the frequency
%   vector. CENTERDC can be one of these two values [ {false} | true ]. 
%
%   PSDFREQVEC('Range',RANGE) specifies the range of frequency in RANGE.
%   RANGE can be one of the two strings [ {whole} | half ]. Assuming
%   CenterDC=false then:
%       'whole' = [0, 2pi)
%       'half'  = [0, pi] for even NPTS or [0, pi) for odd NPTS
%
%   When CenterDC=true then:
%       'whole' = (-pi, pi] for even NPTs or (-pi, pi) for odd NPTs
%       'half'  = [-pi/2, pi/2] for even* NPTS or (-pi/2, pi/2) for odd NPTS
%
%       *When NPTS is not divisible by 4, then the range is (-pi/2, pi/2).
%
%   When Range='half' the frequency vector has length (NPTS/2+1) if NPTS is
%   even**, and (NPTS+1)/2 if NPTS is odd***.
%
%       **If CenterDc=true and the number of points specified is even is
%       not divisible by 4, then the number of points returned is NPTS/2.
%       This is to avoid frequency points outside the range [-pi/2 pi/2]. 
%
%       ***If CenterDC=true and the number of points NPTS specified is odd
%       and (NPTS+1)/2 is even then the length of the frequency vector is
%       (NPTS-1)/2.

%#codegen

%   Copyright 1988-2023 The MathWorks, Inc.
coder.internal.prefer_const(varargin)
[Npts,Fs,Range,CenterDC] = local_pvparse(varargin{:});
w = frequencygrid(Range,Npts(1),Fs,CenterDC);

end

%--------------------------------------------------------------------------

function w = frequencygrid(Range,Npts,Fs,CenterDC)
% Compute the frequency grid.
coder.inline('always');
coder.internal.prefer_const(Range,Npts,CenterDC)
% Compute the whole frequency range, e.g., [0,2pi) to avoid round off errors.
if isempty(Fs) || isnan(Fs)
    Fs1 = cast(2*pi,'like',Fs);
else
    Fs1 = Fs;
end

freq_res = Fs1/Npts;
w1 = freq_res*(0:Npts-1);

% There can still be some minor round off errors in the frequency grid.  
% Fix the known points, i.e., those near pi and 2pi.
Nyq = Fs1/2;
half_res = freq_res/2; 

% Determine if Npts is odd and calculate half and quarter of Npts.
[isNPTSodd,halfNPTS,ishalfNPTSodd,quarterNPTS] = NPTSinfo(Npts);

if isNPTSodd
    % Adjust points on either side of Nyquist.
    w1(halfNPTS)   = Nyq - half_res;
    w1(halfNPTS+1) = Nyq + half_res;
else
    % Make sure we hit Nyquist exactly, i.e., pi or Fs/2 
    w1(halfNPTS) = Nyq;
end
w1(Npts) = Fs1-freq_res;

% Get the right grid based on range, centerdc, etc.
w = finalgrid(w1,Npts,Nyq,Range,CenterDC,isNPTSodd,ishalfNPTSodd,halfNPTS,quarterNPTS);

end

%--------------------------------------------------------------------------

function [isNPTSodd,halfNPTS,ishalfNPTSodd,quarterNPTS] = NPTSinfo(NPTS)
% Determine if we're dealing with even or odd lengths of NPTS, 1/2 NPTS,
% and 1/4 NPTS.

coder.inline('always')

isNPTSodd = false;
if rem(NPTS,2)
    isNPTSodd = true;
end

% Determine half the number of points.
if isNPTSodd
    halfNPTS = (NPTS+1)/2;  
else
    halfNPTS = (NPTS/2)+1; 
end


ishalfNPTSodd = false;     
if rem(halfNPTS,2)
    ishalfNPTSodd = true;  
end


if ishalfNPTSodd
    quarterNPTS = (halfNPTS+1)/2; 
else
    quarterNPTS = (halfNPTS/2)+1;
end

end

%--------------------------------------------------------------------------

function w2 = finalgrid(w1,Npts,Nyq,Range,CenterDC,isNPTSodd,ishalfNPTSodd,halfNPTS,quarterNPTS)
% Calculate the correct grid based on user specified values for range,
% centerdc, etc.

coder.inline('always');

% Calculated by default.% [0, 2pi)

switch lower(Range)
    case 'whole'
        w = coder.nullcopy(zeros(Npts,1,'like',w1));
        
        if CenterDC          % (-pi, pi] even or (-pi, pi) odd
            if isNPTSodd
                negEndPt = halfNPTS;
            else
                negEndPt = halfNPTS-1;
            end
            
            w = [-fliplr(w1(2:negEndPt)), w1(1:halfNPTS)];
        else
            
            w = w1(:);  % Return a column vector.
        end
        
        
    case 'half'
        
        w = w1(1:halfNPTS);      % [0, pi] even or [0, pi) odd
        
        % For even number of points that are not divisible by 4 you get
        % less one point to avoid going outside the [-pi/2 pi/2] range.
        if CenterDC            % [-pi/2, pi/2] even (-pi/2, pi/2) odd
            if ishalfNPTSodd
                negEndPt = quarterNPTS;
            else
                quarterNPTS = quarterNPTS-1; % Avoid going over pi/2
                negEndPt = quarterNPTS;
            end
            
            w = [-fliplr(w1(2:negEndPt)), w1(1:quarterNPTS)];
            
            
            if ~rem(Npts,4)
                % Make sure we hit pi/2 exactly when Npts is divisible
                % by 4! In this case it's due to roundoff.
                w(end) = Nyq/2;
            end
        end
    otherwise
        coder.internal.error('siglib:psdfreqvec:InternalError');
        
end

w2 = w(:);

end

%--------------------------------------------------------------------------
function [Npts,Fs,Range,CenterDC] = local_pvparse(varargin)
coder.internal.prefer_const(varargin)
% Define default parameters
defaultNpts = 1024;
defaultFs = 2*pi;
defaultRange = 'whole';
defaultCenterDC = false;

if coder.target('MATLAB')
    p = inputParser;
    p.addParameter('npts',defaultNpts);
    p.addParameter('Fs',defaultFs);
    p.addParameter('Range',defaultRange);
    p.addParameter('CenterDC',defaultCenterDC);
    p.parse(varargin{:});
    Npts = p.Results.npts;
    Fs = p.Results.Fs;
    Range = p.Results.Range;
    CenterDC = p.Results.CenterDC;
else
    parms = struct(...
        'npts',uint32(0),...
        'Fs',uint32(0),...
        'Range',uint32(0),...
        'CenterDC',uint32(0));
    pstruct = coder.internal.parseParameterInputs(parms,[], varargin{:});
    Npts = coder.internal.getParameterValue(pstruct.npts, defaultNpts, varargin{:});
    Fs = coder.internal.getParameterValue(pstruct.Fs, defaultFs, varargin{:});
    Range = coder.internal.getParameterValue(pstruct.Range, defaultRange, varargin{:});
    CenterDC = coder.internal.getParameterValue(pstruct.CenterDC, defaultCenterDC, varargin{:});
end

end

% LocalWords:  Npts NPTS Fs FS roundoff siglib npts
