function [Magnitude, Frequency, Ts] = computeSpectrum(sysc, specVar, wspec, NominalRPM, IsPS)
% Compute power spectrum for a cell array of LTIs.
% NominalRPM: used if "sysc" stores order spectrum.
% IsPS: is the computation for power spectrum (IsPS=true) or for order
% spectrum (IsPS=false). For order spectrum, the frequencies are orders
% themselves and no scaling is required. If NominalRPM is specified and
% IsPS is true, that means we are trying to put order spectrum data on a
% power spectrum plot. This is not user facing right now.

%  Copyright 1986-2018 The MathWorks, Inc.

if ~iscell(sysc), sysc = {sysc}; end
Magnitude = cell(1,numel(sysc)); Frequency = cell(1,numel(sysc));
for ct = 1:numel(sysc)
   if isa(sysc{1},'idfrd')
      sys = sysc{ct};
      if isempty(sys) || isa(sys,'missing')
          Frequency{ct} = NaN;
          Magnitude{ct} = NaN;
          Ts = 1;
          continue;
      end
      sysdata = getPlotLTIData(sys);
      if ~isa(sysdata, 'idpack.frddata') && isa(sysdata,'ltipack.ltidata')
          sysdata = idss(sysdata);
      elseif isa(sysdata, 'ltipack.frddata') && ~isnan(NominalRPM) && IsPS
          sysdata.Frequency = sysdata.Frequency*NominalRPM/60*2*pi;
      end
%       [~,nu] = iosize(sysdata); NOT Needed anymore
      
      % Default behavior: regenerate data on appropriate grid based on input
      % arguments
      %g2329971:  We do not need outputSpectrum which uses ResponseData
      %because it is not supported. We are interested in Spectrum only (We
      % are only supporting the case where ResponseData is empty => 
      %(y = Gu + He) : we are ignoring 'Gu' here, even if it is present in the
      %data). 
      
      [ps,w] = noiseSpectrumSpec(sysdata,3,wspec,true);
      
      % store in rad/s units
      Frequency{ct} = w*funitconv('rad/TimeUnit','rad/s',sys.TimeUnit);
      Magnitude{ct} = abs(ps);
      Ts = sys.Ts*tunitconv(sys.TimeUnit,'seconds');
   else
      tbldata = sysc{ct};
      ps = tbldata.(specVar.Variables(1).Name);
      w = tbldata.(specVar.IndependentVariables(1).Name);
      freqUnit = specVar.IndependentVariables.Unit;
      if specVar.Type == predmaint.internal.pmdata.var.enum.VariableType.OrderSpectrum
          wFactor = 2*pi;
      elseif specVar.IndependentVariables.IVType == predmaint.internal.pmdata.var.enum.IndependentVariableType.Time
          if startsWith(freqUnit, 'rad/')
              wFactor = funitconv('rad/TimeUnit', 'rad/s', specVar.IVUnit);
          elseif startsWith(freqUnit, 'cycles/')
              wFactor = funitconv('cycles/TimeUnit', 'rad/s', specVar.IVUnit);
          else
              wFactor = funitconv(freqUnit, 'rad/s');
          end
      else
          % specVar's independent variable's unit must be 'cycles/XX' or
          % 'rad/XX' for non-time-based
          wFactor = 1;
          if startsWith(freqUnit,"cycles/")
              wFactor = 2*pi;
          end
      end

      % store in rad/s units
      Frequency{ct} = wFactor*w;
      Magnitude{ct} = abs(ps);

      % convert Ts to seconds
      FsFactor = wFactor/(2*pi);
      if ~isempty(tbldata.Properties.CustomProperties.SampleFrequency)
          Ts = 1/(FsFactor*tbldata.Properties.CustomProperties.SampleFrequency);
      else
          % If somebody import a table of spectrum data, without sampling
          % frequency, need to guess a sampling frequency. But the
          % following guess can be very wrong, need to investigate more on
          % this situation.
          Ts = 1/(FsFactor*2*w(end));
      end
   end
end
