clc,clear;
[freq,H]=freqAndHFunction();
freq=[0,freq];
magnitude=[0,H];
% 采样率
Fs= 16000;
% 进行插值
interp_freq = linspace(0, 8000, 160); % 选择插值点数
interp_magnitude = interp1(freq, magnitude, interp_freq, 'linear');
% PSO选项
options = optimoptions('particleswarm', 'SwarmSize', 150, 'MaxIterations', 2000);
% 参数界限
lb = -1e-7 * ones(1, 162); % 下界
ub = 1e-7 * ones(1, 162); % 上界
% 运行PSO
best_params = particleswarm(@objective_function, 162, lb, ub, options);
b_estimated = best_params(1:41);
a_estimated = best_params(42:82);
% 计算估计滤波器的频率响应
[H_est,~] = freqz(b_estimated, a_estimated, interp_freq, 16000);
% 可视化结果
figure;
plot(interp_freq, 20*log10(abs(H_est)), 'r'); % 估计的响应
hold on;
plot(interp_freq, 20*log10(interp_magnitude), 'b'); %原始响应
legend('Estimated Response', 'Original Response');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
function error = objective_function(params)
try
b = params(1:41); % 分子系数
a = params(42:82); % 分母系数
% 计算当前IIR滤波器的频率响应
[H,~] = freqz(b, a, interp_freq, 16000);
% 计算误差
error = sum((abs(H) - interp_magnitude).^2);
% 检查滤波器的稳定性
if ~isstable(b, a)
error = inf; % 给不稳定的滤波器一个大的错误值
end
catch
error = inf; % 对于任何其他错误，返回一个大的错误值
end
end
function [f, H] = freqAndHFunction()
% Define freq array
f = [1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500, 3600, 3700, 3800, 3900, 4000, 4100, 4200, 4300, 4400, 4500, 4600, 4700, 4800, 4900, 5000, 5100, 5200, 5300, 5400, 5500, 5600, 5700, 5800, 5900, 6000, 6100, 6200, 6300, 6400, 6500, 6600, 6700, 6800, 6900, 7000, 7100, 7200, 7300, 7400, 7500, 7600, 7700, 7800, 7900, 8000];
% Define tvr array
tvr = [128.5163628, 129.051837, 128.0333238, 127.200628, 128.4267629, 129.4437987, 130.9112373, 133.5186641, 135.2522284, 134.405548, 130.9848593, 128.2534974, 134.0918335, 138.610268, 140.9095005, 142.6631096, 143.7293213, 144.7608089, 145.3854305, 145.5075137, 145.554636, 144.7531332, 143.9672297, 143.410465, 142.2522647, 141.169289, 140.4608309, 139.6028821, 138.8584932, 138.4957287, 138.1202031, 137.6812736, 137.4824973, 137.3362914, 137.3640559, 137.3579424, 137.6974913, 138.0771592, 138.0983559, 137.5053153, 136.0987116, 133.593733, 131.4955087, 130.7312329, 131.25016, 131.7451198, 132.140601, 132.1775175, 131.4427206, 130.7975073, 129.799245, 128.9741335, 129.4695932, 130.020389, 129.89581, 129.3831058, 127.7675183, 124.4633058, 118.3829614, 108.3285636, 113.1397294, 115.693913, 115.7284724, 114.169323, 116.192023, 120.3573138];
% Calculate H array
H = 10.^((tvr-120)/20);
end
