function [sp,ce] = ampModulate(ss, fc, BW, modIndex, fs)
% ampModulate - amplitude modulate a baseband signal
%
% Inputs:
%   ss - input smaples (maybe complex)
%   fc - carrier frequency (in Hertz)
%   BW - baseband bandwidth (in Hertz); ss will be fitered to enforce BW
%   modIndex - modulation index (as a fraction, e.g., 0.5 for 50%)
%   fs - sampling rate (in Hertz)
%
% Outputs:
%   sp - passband, modulated signal
%   ce - complex envelope
%
% Usage:
%   sp = ampModulate(ss, fc, BW, modIndex, fs)

%% check inputs
if fc+BW >= fs/2
    error('Carrier frequency (%f) plus bandwidth (%d) must be less than the Nyquist frequency (%f)', fc, BW, fs/2);
end

if modIndex == 0
    error('Modulation index can''t be zero! Set modIndex = inf to suppress carrier.')
end

%% low pass filter the input signal to ensure BW
f2 = min(BW*1.1,fs/2*0.95); % end of transition band, aim for 10% above BW

% make a lowpass filter using firpm
hh = firpm(50, [0 BW f2 fs/2]/(fs/2), [1 1 0 0]);

% apply lowpass filter 
ss = conv(hh, ss); % could also use filter(hh, 1, ss);

%% normalize input signal
% find scaling to make min(ss) == -1
scale = abs(min(real(ss)));

if ~isreal(ss)
    % if the signal is complex, check imaginary part of signal for scaling
    scaleI = abs(min(imag(ss)));
    
    scale = max(scale, scaleI);
end

ce = ss/scale;

%% add carrier
if ~isinf(modIndex)
    % figure out carrier amplitude
    Ac = 1/modIndex;
    
    % add carrier to signal
    ce = ce + Ac;
end

%% upconversion
tt = (0:length(ss)-1)/fs; % time axis
% convert complex envelope to passband signal
sp = real(ce .* exp(2j*pi*fc*tt));