%% 参数设置Fs = 44100; % 采样频率 (Hz)fcut_list = [10, 21; % 滤波器1截止频率 [低, 高] (Hz) 21, 42; 42, 85; 85, 170; 170, 500; 500, 1000; 1000, 2000; % 滤波器2截止频率 2000, 4000; 4000, 8000 ]; % 滤波器3截止频率N = 4096; % 滤波器阶数windowType = 'hamming'; % 窗函数类型%windowType = 'hann'; % 窗函数类型%windowType = 'kaiser'; % 窗函数类型%windowType = 'blackman'; % 窗函数类型N_fft = 8192; % 频响分析点数%% 窗函数参数表(主瓣系数和旁瓣衰减)windowParams = struct(... 'blackman', struct('C',12, 'attenuation',-74),... 'hamming', struct('C',4, 'attenuation',-41),... 'hann', struct('C',4, 'attenuation',-31),... 'kaiser', struct('C',4.5,'attenuation',-53));%% 参数校验if ~isfield(windowParams, windowType) error('不支持的窗类型,可选: blackman/hamming/hann/kaiser');end%% 生成窗函数if strcmpi(windowType, 'kaiser') beta = 5; % Kaiser窗参数 win = kaiser(N+1, beta);else win = feval(lower(windowType), N+1); end%% 设计三个并联滤波器num_filters = size(fcut_list,1);b = cell(num_filters, 1);% for k = 1:num_filters% if any(fcut_list(k,:) >= Fs/2)% error('滤波器%d截止频率超过Nyquist频率(%.1fHz)', k, Fs/2);% end% Wn = fcut_list(k,:)/(Fs/2);% b{k} = fir1(N, Wn, 'bandpass', win); % 设计带通滤波器% endmin_transition = 43.07; % 最小过渡带宽度(Hz)transition_ratio = 0.1; % 过渡带占带宽比例for k = 1:num_filters % 1. 安全检查 fLow = fcut_list(k,1); fHigh = fcut_list(k,2); % 截止频率有效性验证 if fHigh >= Fs/2 error('滤波器%d截止频率超过Nyquist频率(%.1fHz)', k, Fs/2); end if fLow >= fHigh error('滤波器%d截止频率范围无效(fLow >= fHigh)', k); end % 2. 计算带宽和过渡带 BW = fHigh - fLow; % 通带带宽 delta_f = max(transition_ratio * BW, min_transition); % 动态过渡带 disp(delta_f); % 3. 自动计算阶数 (基于Hamming窗公式) N = ceil(4 * Fs / delta_f); % 主瓣宽度约束 N = N + mod(N,2); % 强制偶阶数 disp(N); % 4. 带通滤波器设计 Wn = [fLow fHigh]/(Fs/2); % 归一化截止频率 b{k} = fir1(N, Wn, 'bandpass', hamming(N+1)); % 动态阶数设计 % 5. 设计参数输出 fprintf('滤波器%d: BW=%.1fHz, N=%d, Δf=%.1fHz\n',... k, BW, N, delta_f);end%% 幅频响应分析与绘图figure('Position',[100,100,1400,900],'Color','w');colors = {'#E74C3C','#2ECC71','#3498DB'}; % 红/绿/蓝配色H_total = zeros(N_fft,1);lineStyles = {'--', ':', '-.'}; % 不同滤波器的虚线样式% 绘制单个滤波器响应for k = 1:num_filters [H, f] = freqz(b{k}, 1, N_fft, Fs); H_mag = 20*log10(abs(H)); % 计算3dB带宽 try [bw, f_low, f_high, threshold] = calculate_3dB_bandwidth(f, H_mag, true); catch ME warning('滤波器%d带宽计算失败: %s', k, ME.message); continue; end % 绘制幅频响应 semilogx(f, H_mag, 'Color', colors{k}, 'LineWidth', 1.5,... 'DisplayName', sprintf('Filter %d',k)); hold on; % 标记3dB点坐标 plot(f_low, threshold, 'o', 'MarkerSize',6,... 'MarkerFaceColor',colors{k}, 'MarkerEdgeColor','k'); text(f_low, threshold-4, sprintf('%.1fHz\\n%.1fdB',f_low,threshold),... 'Color',colors{k}, 'FontSize',8, 'VerticalAlignment','top'); plot(f_high, threshold, 'o', 'MarkerSize',6,... 'MarkerFaceColor',colors{k}, 'MarkerEdgeColor','k'); text(f_high, threshold-4, sprintf('%.1fHz\\n%.1fdB',f_high,threshold),... 'Color',colors{k}, 'FontSize',8, 'VerticalAlignment','top'); % 标记设计截止频率 design_fc = fcut_list(k,:); line([design_fc(1) design_fc(1)], ylim,... 'LineStyle',':', 'Color',[0.4 0.4 0.4], 'LineWidth',1.2); text(design_fc(1), max(ylim)-5, sprintf('设计低端\\n%.1fHz',design_fc(1)),... 'Color',[0.4 0.4 0.4], 'FontSize',8, 'HorizontalAlignment','center'); line([design_fc(2) design_fc(2)], ylim,... 'LineStyle',':', 'Color',[0.4 0.4 0.4], 'LineWidth',1.2); text(design_fc(2), max(ylim)-5, sprintf('设计高端\\n%.1fHz',design_fc(2)),... 'Color',[0.4 0.4 0.4], 'FontSize',8, 'HorizontalAlignment','center'); % 标注带宽值 text(mean([f_low, f_high]), threshold-8,... sprintf('BW_{3dB} = %.1f Hz', bw),... 'Color', colors{k}, 'HorizontalAlignment', 'center',... 'FontSize', 9, 'BackgroundColor','w'); H_total = H_total + H; end%% 绘制总响应H_total_abs = abs(H_total);H_total_db = 20*log10(H_total_abs);semilogx(f, H_total_db, 'k', 'LineWidth', 2,... 'DisplayName', 'Combined Response');%% 图形美化xlabel('Frequency (Hz)');ylabel('Magnitude (dB)');title(['Parallel Filter Bank Response (', windowType, ' Window)']);set(gca, 'XScale','log', 'XMinorGrid','on', 'XMinorTick','on',... 'XTick',[20 50 100 200 500 1000 2000 5000 10000 20000],... 'XTickLabel',{'20','50','100','200','500','1k','2k','5k','10k','20k'});grid on;legend('Location','best');xlim([20, 20000]);ylim([-120, 5]);%% 改进的3dB带宽计算函数function [bw_3db, f_low, f_high, threshold] = calculate_3dB_bandwidth(f, H_mag, is_dB)% 输入参数校验validateattributes(f, {'numeric'}, {'vector'});validateattributes(H_mag, {'numeric'}, {'vector','size',size(f)});if nargin < 3, is_dB = false; end% 单位转换if ~is_dB H_mag = 20*log10(abs(H_mag) + eps); end% 峰值定位与插值[~, idx] = max(H_mag);if idx <= 1 || idx >= numel(H_mag) error('峰值位于端点');end% 二次多项式插值delta = 0.5*(H_mag(idx-1)-H_mag(idx+1))/(2*H_mag(idx)-H_mag(idx-1)-H_mag(idx+1));peak_refined = H_mag(idx) - 0.25*(H_mag(idx-1)-H_mag(idx+1))*delta;threshold = peak_refined - 3;% 低频交点插值idx_low = find(diff(sign(H_mag(1:idx) - threshold)) > 0, 1, 'last');if isempty(idx_low) error('未找到低频交点');endx_low = f(idx_low:idx_low+1);y_low = H_mag(idx_low:idx_low+1) - threshold;f_low = interp1(y_low, x_low, 0, 'linear');% 高频交点插值idx_high = find(diff(sign(H_mag(idx:end) - threshold)) < 0, 1, 'first') + idx -1;if isempty(idx_high) error('未找到高频交点');endx_high = f(idx_high:idx_high+1);y_high = H_mag(idx_high:idx_high+1) - threshold;f_high = interp1(y_high, x_high, 0, 'linear');% 带宽计算bw_3db = f_high - f_low;% 结果校验if bw_3db <=0 || bw_3db > f(end)-f(1) error('无效带宽值');endend,上面代码报错。> 位置:bandpassfiltertest2 (第 273 行) 索引超过数组元素的数量。索引不能超过 3。出错 bandpassfiltertest2 (第 278 行) semilogx(f, H_mag, 'Color', colors{k}, 'LineWidth', 1.5,...,