function [freq, spec] = SpectrumAnalyzer(signal_in, Fs, Fc, ResBW, PlotColor)
%SpectrumAnalyzer Computes the power spectral density function of a signal
%Inputs:
%    signal_in    -    The complex input signal
%    Fs            -    The sampling frequency (in Hz).
%    Fc            -    The carrier frequency (in Hz). If empty, defaults to 0.
%    ResBW        -    Resolution Bandwidth of the spectrum. If empty, uses a
%                    ResBW = Fs/1024.
%    PlotColor    -    The plot color to use. Same format as plot()
%Outputs: (if no outputs specified, plots on the current axes).
%    freq        -    Frequency x-axis
%    spec        -    Spectrum y-axis.

% Notes:
% By defn, mean of power spectrum must equal mean of input signal. To
% check, use mean(abs(signal_in).^2) versus mean(10.^((spec-30)/10) * 100)

% CVS
% $Revision: 1.2 $ $Date: 2008/09/26 18:59:35 $

% Compute the number of points needed, according to resbw.
iResBWPoints = 1024;
if(exist('ResBW'))
    if(~isempty(ResBW))
        iResBWPoints = ceil(Fs/ResBW);
        % Compensate, always have odd number points for even number bins.
        if(floor(iResBWPoints/2) == iResBWPoints/2)
            iResBWPoints = iResBWPoints + 1;
        end
    end
end

iWindowPoints = iResBWPoints;

% Create my windowing signal
windowsignal = hann(iWindowPoints);

% Compute number of 'blocks', rounded down.
% Change this if needed.
iOverlap = floor(iResBWPoints / 2);

% Use while loop instead of for if overlap is needed.
iIndex = 1;
iBlocks = 0;
signal_in = signal_in(:);

% Initialize accumulation vector
power_mag_accum = zeros(iWindowPoints, 1);

while(iIndex + iWindowPoints < length(signal_in))

    % use chirp-z transform
    block_mag = abs(czt(complex(windowsignal,windowsignal) .* signal_in(iIndex:(iIndex + iWindowPoints - 1)))).^2;
    
    % find current spectral power
    block_mag_mean = mean(block_mag);
    
    % find mean power of pre-windowed signal
    input_mag_mean = mean(abs(signal_in(iIndex:(iIndex + iWindowPoints - 1)).^2));
    
    % normalize, use / 10 since squared already
    norm_block_mag = block_mag * input_mag_mean / block_mag_mean;
    
    % check again (debugging)
    block_mag_mean = mean(norm_block_mag);
    
    % accumulate
    power_mag_accum = power_mag_accum + norm_block_mag;
    
    % Increment
    iIndex = iIndex + iWindowPoints - iOverlap;
    iBlocks = iBlocks + 1;
end

% Averaging all the points
% There is a reason why i multiply by 4.
power_mag_accum = power_mag_accum * 4 / iBlocks / iWindowPoints;
power_dBm = 10 * log10(fftshift(abs(power_mag_accum))/100) + 30;

% return frequency and power spectrum
freq = (Fc-Fs/2):(Fs/(iWindowPoints-1)):(Fc+Fs/2);
spec = power_dBm;

if nargout == 0
    if(exist('PlotColor', 'var'))
        %%figure;plot(freq, spec, PlotColor);
        plot(freq, spec, PlotColor);
        xlabel('frequency (Hz)')
        ylabel('PSD (dBm/Hz)')
    else
        %%figure;plot(freq, spec, '-b');
        figure; plot(freq, spec, '-b');
        xlabel('frequency (Hz)')
        ylabel('PSD (dBm/Hz)')
    end
end

return;
