% Define center frequencies
center_freqs = [31.5, 63, 125, 250, 500, 1000, 2000, 4000, 8000, 16000];
num_filters = length(center_freqs);

% Set the sample frequency
Fs = 16000;

% Define the quality factor (Q)
Q = 5;

% Calculate normalized cutoff frequencies
normalized_cutoffs = center_freqs / (Fs / 2 * Q);

% Initialize the figure
figure;

% Loop through each center frequency and create filters
for i = 1:num_filters
    % Calculate cutoff frequencies
    cutoff_low = normalized_cutoffs(i) / sqrt(2);
    cutoff_high = normalized_cutoffs(i) * sqrt(2);

    % Design filter using fir1
    b = fir1(30, [cutoff_low, cutoff_high], 'bandpass');

    % Get frequency response using freqz
    [H, Freq] = freqz(b, 1, 1024, Fs);

    % Plot the magnitude components of the filter
    subplot(2, 1, 1);
    plot(Freq, abs(H), 'DisplayName', sprintf('Center Frequency: %g Hz', center_freqs(i)));
    hold on;

    % Plot the filter coefficients
    subplot(2, 1, 2);
    plot(b, 'DisplayName', sprintf('Center Frequency: %g Hz', center_freqs(i)));
    hold on;
end

% Configure subplot 1
subplot(2, 1, 1);
title('Magnitude Components of Filters');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
legend;
grid on;

% Configure subplot 2
subplot(2, 1, 2);
title('Filter Coefficients');
xlabel('Coefficient Index');
ylabel('Coefficient Value');
legend;
grid on;

% Initialize a new figure for the sum of frequency responses
figure;
sum_response = zeros(size(H));

% Loop through each center frequency and add frequency responses
for i = 1:num_filters
    % Design filter using fir1
    b = fir1(30, [normalized_cutoffs(i) / sqrt(2), normalized_cutoffs(i) * sqrt(2)], 'bandpass');

    % Get frequency response using freqz
    [H, Freq] = freqz(b, 1, 1024, Fs);

    % Add frequency response to the sum
    sum_response = sum_response + abs(H);
end

% Plot the sum of frequency responses
plot(Freq, sum_response, 'LineWidth', 2);
title('Sum of Frequency Responses');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;