function showEQNoisePSD(eq_noise, options)
arguments
    eq_noise 
    options.postfilter_taps = NaN
    options.fignum (1,1) double = NaN % Default to NaN if not provided
    options.displayname (1,:) char = '' % Default to an empty string if not provided
    options.color = [];
    options.colormode (1,1) string {mustBeMember(options.colormode, ["diverging", "qualitative", "paired"])} = "diverging";
end

    persistent nextGroupIdx
    if isempty(nextGroupIdx)
        nextGroupIdx = 0;
    end

    % Determine the figure number to use or create a new figure
    if isnan(options.fignum)
        fig = figure; % Create a new figure and get its handle
    else
        fig = figure(options.fignum); % Use the specified figure number
    end
    hold on
    ax = gca;
    linesBefore = findall(ax, 'Type', 'Line');

    useAutomaticColor = isempty(options.color);
    if useAutomaticColor
        options.color = selectEQNoisePSDColor(ax, options.colormode);
    end

    % Ensure the figure is ready before calling spectrum
    eq_noise = eq_noise - mean(eq_noise.signal);
    eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color,"fft_length",4096);

    if ~isnan(options.postfilter_taps)
        % Hold on to the figure for further plotting
        hold on;
    
        % Compute the frequency response of the postfilter
        [h, w] = freqz(1, options.postfilter_taps, length(eq_noise), "whole", eq_noise.fs);
        h = h / max(abs(h)); % Normalize the filter response
    
        % Adjust frequency axis to center at 0
        w_ = (w - eq_noise.fs / 2);
    
        % Plot the inverted postfilter response
        plot(w_ * 1e-9, 20 * log10(fftshift(abs(h))), 'DisplayName', ['Burg Coeffs: ', num2str(round(options.postfilter_taps, 2)), ' '], 'LineWidth', 1,'Color',options.color,'LineStyle','--');
    
        % Ensure a legend is displayed
        legend('show');
    end

    if useAutomaticColor
        nextGroupIdx = nextGroupIdx + 1;
        newLines = getNewLines(ax, linesBefore);
        tagEQNoisePSDLines(newLines, nextGroupIdx);
        recolorEQNoisePSDLines(ax, options.colormode);
    end

    xlim([-eq_noise.fs/2* 1e-9 eq_noise.fs/2* 1e-9]);
    % ylim([-15, 0]);

end

function color = selectEQNoisePSDColor(ax, colormode)
if colormode == "qualitative"
    N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children));
    cmap = linspecer(8);
    color = cmap(mod(N, size(cmap, 1)) + 1, :);
else
    color = [0 0 0];
end
end

function newLines = getNewLines(ax, linesBefore)
linesAfter = findall(ax, 'Type', 'Line');
newLines = linesAfter(~ismember(linesAfter, linesBefore));
end

function tagEQNoisePSDLines(lines, groupIdx)
for lineIdx = 1:numel(lines)
    setappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx', groupIdx);
end
end

function recolorEQNoisePSDLines(ax, colormode)
if colormode == "qualitative"
    return
end

lines = findall(ax, 'Type', 'Line');
groupIdx = NaN(size(lines));

for lineIdx = 1:numel(lines)
    if isappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx')
        groupIdx(lineIdx) = getappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx');
    end
end

lines = lines(~isnan(groupIdx));
groupIdx = groupIdx(~isnan(groupIdx));

if isempty(lines)
    return
end

groups = unique(groupIdx, 'stable');
groups = sort(groups);
cols = getEQNoisePSDColormap(colormode, numel(groups));

for groupColorIdx = 1:numel(groups)
    sameGroup = groupIdx == groups(groupColorIdx);
    set(lines(sameGroup), 'Color', cols(groupColorIdx, :));
end
end

function cols = getEQNoisePSDColormap(colormode, numColors)
switch colormode
    case "paired"
        cols = cbrewer2('paired', numColors);
    otherwise
        cols = getDivergingWithoutBrightCenter(numColors);
end
end

function cols = getDivergingWithoutBrightCenter(numColors)
numBaseColors = max(6, numColors);
centerExcludeFraction = 0.3;

baseCols = cbrewer2('RdYlBu', numBaseColors);
centerIdx = (numBaseColors + 1) / 2;
excludeHalfWidth = centerExcludeFraction * numBaseColors / 2;
keepIdx = find(abs((1:numBaseColors) - centerIdx) > excludeHalfWidth);
sampleIdx = round(linspace(1, numel(keepIdx), numColors));

cols = baseCols(keepIdx(sampleIdx), :);
end