数字带阻滤波器，满足如下技术指标要求:对模拟信号在1.9KHz、4.7kHz处的衰减大于40 dB，在1.2 kHz、9.2 kHz 处的衰减小于 1dB。请分别设计R、FIR(要求具有线性相位)类型

本文链接：https://blog.youkuaiyun.com/m0_73811793/article/details/144654925

clc;
close all;
clear;

% 设计数字带阻滤波器
fs = 20000;            % 采样频率20kHz

% 定义关键频率点
f_check = [1200 1900 4700 9200];  % Hz
w_norm = f_check/(fs/2);          % 转换为归一化频率

% 定义通带和阻带的边界频率（归一化频率）
wp1 = w_norm(1);      % 第一通带频率 (1.2kHz)
ws1 = w_norm(2);      % 第一阻带频率 (1.9kHz)
ws2 = w_norm(3);      % 第二阻带频率 (4.7kHz)
wp2 = w_norm(4);      % 第二通带频率 (9.2kHz)

% 计算过渡带宽度和中心频率
deltaw = min((ws1-wp1), (wp2-ws2));
wc1 = (wp1+ws1)/2;
wc2 = (ws2+wp2)/2;

% 计算滤波器阶数
N0 = ceil(20/deltaw);  
N = N0 + mod(N0+1, 2);
beta = 8;             % Kaiser窗参数

% 设计滤波器
h = fir1(N-1, [wc1, wc2], 'stop', kaiser(N, beta));

% 1. 绘制零极点图
figure(1)
zplane(h, 1)
title('零极点分布图')
grid on

% 2. 绘制单位脉冲响应
figure(2)
stem(0:length(h)-1, h, 'filled')
title('单位脉冲响应')
xlabel('n')
ylabel('h(n)')
grid on

% 3. 计算频率响应
[h_resp, w, phi] = freqz(h, 1, 1024, fs);
mag_db = 20*log10(abs(h_resp));
phase = unwrap(angle(h_resp));  % 解卷绕相位
gd = -diff(phase)./diff(w);     % 计算群延时

% 绘制完整的频率响应（幅度和相位）
figure(4)
subplot(3,1,1)
plot(w, mag_db)
grid on
title('幅度响应')
xlabel('频率 (Hz)')
ylabel('幅度 (dB)')

subplot(3,1,2)
plot(w, phase)
grid on
title('相位响应')
xlabel('频率 (Hz)')
ylabel('相位 (弧度)')

subplot(3,1,3)
plot(w(1:end-1), gd)
grid on
title('群延时')
xlabel('频率 (Hz)')
ylabel('群延时 (样本)')

% 4. 计算关键频率点的响应
[h_resp_check, w_check] = freqz(h, 1, w_norm*pi);
mag_db_check = 20*log10(abs(h_resp_check));

% 输出系统特性分析
fprintf('\n系统特性分析：\n');
fprintf('1. 滤波器阶数：%d\n', length(h)-1);
fprintf('2. 滤波器类型：FIR带阻滤波器\n');
fprintf('3. 系统特点：\n');
fprintf('   - 线性相位（FIR对称设计）\n');
fprintf('   - 因为是FIR，所以系统稳定\n');
fprintf('   - 所有极点都在原点\n\n');

% 输出相位特性分析
fprintf('\n相位特性分析：\n');
fprintf('1. 相位类型：线性相位（I型FIR，长度为奇数，系数对称）\n');
fprintf('2. 滤波器长度：%d（奇数）\n', length(h));
fprintf('3. 群延时：%.2f 个样本\n', mean(gd));
fprintf('4. 相位特点：\n');
fprintf('   - 严格线性相位\n');
fprintf('   - 恒定群延时\n');
fprintf('   - 中心对称的系数\n\n');

% 验证系数对称性
fprintf('系数对称性验证：\n');
max_diff = max(abs(h - flip(h)));
fprintf('最大对称差异：%.2e\n', max_diff);

% 检查FIR类型
N = length(h);
center = ceil(N/2);
if abs(h(center)) > 1e-10
    fprintf('FIR类型：I型（奇数长度，偶对称）\n');
else
    fprintf('FIR类型：III型（奇数长度，奇对称）\n');
end

% 输出关键频率点的衰减值
fprintf('\n关键频率点的衰减值：\n');
fprintf('1.2 kHz: %.2f dB (要求 < 1 dB)\n', -mag_db_check(1));
fprintf('1.9 kHz: %.2f dB (要求 > 40 dB)\n', -mag_db_check(2));
fprintf('4.7 kHz: %.2f dB (要求 > 40 dB)\n', -mag_db_check(3));
fprintf('9.2 kHz: %.2f dB (要求 < 1 dB)\n', -mag_db_check(4));

% 保存滤波器系数（如果需要）
% save('bandstop_filter.mat', 'h');

clc;
close all;
clear;

% 采样频率
fs = 20000;

% 滤波器设计参数
N = 401;  % 滤波器长度（奇数以确保Type-I线性相位）
f = [0:N-1]*fs/N;  % 频率点

% 设计理想频率响应
H = ones(1,N);
% 定义阻带（归一化后的频率索引）
stop_band1 = round([1900 4700]*N/fs);
% 设置阻带衰减
H(stop_band1(1):stop_band1(2)) = 0;
% 对称填充频谱
H(N-stop_band1(2)+1:N-stop_band1(1)+1) = 0;

% 使用IFFT获取时域冲激响应
h = real(ifft(H));
% 移动到因果系统
h = fftshift(h);

% 应用Hamming窗来减少吉布斯现象
w = hamming(N)';
h = h.*w;

% 1. 绘制零极点图
figure(1)
zplane(h,1)
title('零极点分布图')
grid on

% 2. 绘制单位脉冲响应
figure(2)
stem(0:length(h)-1,h,'filled')
title('单位脉冲响应')
xlabel('n')
ylabel('h(n)')
grid on

% 3. 计算频率响应
[h_resp,w,phi] = freqz(h,1,1024,fs);
mag_db = 20*log10(abs(h_resp));
phase = unwrap(angle(h_resp));
gd = -diff(phase)./diff(w);  % 群延时

% 绘制频率响应
figure(3)
subplot(3,1,1)
plot(w,mag_db)
grid on
title('幅度响应')
xlabel('频率 (Hz)')
ylabel('幅度 (dB)')
yline(-40,'r--')  % 阻带衰减要求
yline(-1,'r--')   % 通带纹波要求

subplot(3,1,2)
plot(w,phase)
grid on
title('相位响应')
xlabel('频率 (Hz)')
ylabel('相位 (弧度)')

subplot(3,1,3)
plot(w(1:end-1),gd)
grid on
title('群延时')
xlabel('频率 (Hz)')
ylabel('群延时 (样本)')

% 检查关键频点的响应
f_check = [1200 1900 4700 9200];
w_norm_check = f_check/(fs/2);
[h_resp_check,w_check] = freqz(h,1,w_norm_check*pi);
mag_db_check = 20*log10(abs(h_resp_check));

% 输出系统分析
fprintf('\n系统特性分析：\n');
fprintf('1. 滤波器阶数：%d\n', length(h)-1);
fprintf('2. 滤波器类型：频率采样法设计的FIR带阻滤波器\n');
fprintf('3. 系统特点：\n');
fprintf(' - 使用Hamming窗函数优化\n');
fprintf(' - 线性相位（Type-I FIR）\n');
fprintf(' - 系统稳定\n\n');

% 输出关键频率点的衰减值
fprintf('\n关键频率点的衰减值：\n');
fprintf('1.2 kHz: %.2f dB (通带)\n', -mag_db_check(1));
fprintf('1.9 kHz: %.2f dB (阻带)\n', -mag_db_check(2));
fprintf('4.7 kHz: %.2f dB (阻带)\n', -mag_db_check(3));
fprintf('9.2 kHz: %.2f dB (通带)\n', -mag_db_check(4));

% 验证线性相位
is_symmetric = all(abs(h - flip(h)) < 1e-10);
fprintf('\n线性相位验证：\n');
if is_symmetric
    fprintf('系数对称性：对称\n');
else
    fprintf('系数对称性：不对称\n');
end
fprintf('FIR类型：Type-I (奇数长度，偶对称)\n');
fprintf('平均群延时：%.2f 样本\n', mean(gd));

clc;
close all;
clear;

% 设计数字IIR带阻滤波器（切比雪夫 Type II）
fs = 20000; % 采样频率20kHz

% 定义关键频率点
f_check = [1200 1900 4700 9200]; % Hz
w_norm = f_check/(fs/2); % 转换为归一化频率

% 定义通带和阻带的边界频率
wp1 = w_norm(1); % 第一通带频率 (1.2kHz)
ws1 = w_norm(2); % 第一阻带频率 (1.9kHz)
ws2 = w_norm(3); % 第二阻带频率 (4.7kHz)
wp2 = w_norm(4); % 第二通带频率 (9.2kHz)

% 设计规格
Rp = 1;     % 通带波纹(dB)
Rs = 45;    % 阻带衰减(dB)，略微提高以确保性能

% 使用cheb2ord函数计算滤波器阶数
[N, Wn] = cheb2ord([ws1 ws2], [wp1 wp2], Rp, Rs);
[b, a] = cheby2(N, Rs, [ws1 ws2], 'stop');

% 1. 绘制零极点图
figure(1)
zplane(b, a)
title('零极点分布图')
grid on

% 2. 计算和绘制频率响应
[h_resp, w] = freqz(b, a, 1024, fs);
mag_db = 20*log10(abs(h_resp));
phase = unwrap(angle(h_resp));
gd = -diff(phase)./diff(w); % 计算群延时

% 绘制完整的频率响应（幅度、相位和群延时）
figure(2)
subplot(3,1,1)
plot(w, mag_db)
grid on
title('幅度响应')
xlabel('频率 (Hz)')
ylabel('幅度 (dB)')
ylim([-60 5])

subplot(3,1,2)
plot(w, phase)
grid on
title('相位响应')
xlabel('频率 (Hz)')
ylabel('相位 (弧度)')

subplot(3,1,3)
plot(w(1:end-1), gd)
grid on
title('群延时')
xlabel('频率 (Hz)')
ylabel('群延时 (样本)')

% 3. 计算关键频率点的响应
[h_resp_check, w_check] = freqz(b, a, w_norm*pi);
mag_db_check = 20*log10(abs(h_resp_check));

% 输出系统特性分析
fprintf('\n系统特性分析：\n');
fprintf('1. 滤波器阶数：%d\n', N);
fprintf('2. 滤波器类型：IIR切比雪夫Type II带阻滤波器\n');
fprintf('3. 系统特点：\n');
fprintf(' - 等纹波阻带响应\n');
fprintf(' - 非线性相位\n');
fprintf(' - 递归结构\n');
fprintf(' - 需要检查稳定性\n\n');

% 检查稳定性
poles = roots(a);
isStable = all(abs(poles) < 1);
fprintf('稳定性分析：\n');
if isStable
    fprintf('系统稳定（所有极点在单位圆内）\n');
else
    fprintf('系统不稳定（存在极点在单位圆外）\n');
end
fprintf('最大极点模值：%.4f\n\n', max(abs(poles)));

% 输出关键频率点的衰减值
fprintf('关键频率点的衰减值：\n');
fprintf('1.2 kHz: %.2f dB (要求 < 1 dB)\n', -mag_db_check(1));
fprintf('1.9 kHz: %.2f dB (要求 > 40 dB)\n', -mag_db_check(2));
fprintf('4.7 kHz: %.2f dB (要求 > 40 dB)\n', -mag_db_check(3));
fprintf('9.2 kHz: %.2f dB (要求 < 1 dB)\n', -mag_db_check(4));

% 计算并显示滤波器的具体特性
fprintf('\n滤波器性能指标：\n');
fprintf('1. 通带波纹：%.2f dB\n', Rp);
fprintf('2. 阻带衰减：%.2f dB\n', Rs);
fprintf('3. 滤波器阶数：%d\n', N);
fprintf('4. 极点数量：%d\n', length(poles));
fprintf('5. 零点数量：%d\n', length(roots(b)));

% 输出系数
fprintf('\n滤波器系数：\n');
fprintf('分子系数 b：\n');
disp(b);
fprintf('分母系数 a：\n');
disp(a);

% 保存滤波器系数
% save('iir_bandstop_filter.mat', 'b', 'a');

clc;
close all;
clear;

% 设计数字IIR带阻滤波器（椭圆）
fs = 20000; % 采样频率20kHz

% 定义关键频率点
f_check = [1200 1900 4700 9200]; % Hz
w_norm = f_check/(fs/2); % 转换为归一化频率

% 定义通带和阻带的边界频率
wp1 = w_norm(1); % 第一通带频率 (1.2kHz)
ws1 = w_norm(2); % 第一阻带频率 (1.9kHz)
ws2 = w_norm(3); % 第二阻带频率 (4.7kHz)
wp2 = w_norm(4); % 第二通带频率 (9.2kHz)

% 修改设计规格以提高性能
Rp = 0.5;    % 降低通带波纹以获得更好的整体性能
Rs = 60;     % 提高阻带衰减要求以确保达到40dB

% 使用ellipord函数计算滤波器阶数
[N, Wn] = ellipord([ws1 ws2], [wp1 wp2], Rp, Rs);
[b, a] = ellip(N, Rp, Rs, [ws1 ws2], 'stop');

% 增加滤波器阶数以提高性能
N = N + 1;  % 增加一阶
[b, a] = ellip(N, Rp, Rs, [ws1 ws2], 'stop');

% [剩余代码与之前相同，包括绘图和分析部分]


% 1. 绘制零极点图
figure(1)
zplane(b, a)
title('零极点分布图')
grid on

% 2. 计算和绘制频率响应
[h_resp, w] = freqz(b, a, 1024, fs);
mag_db = 20*log10(abs(h_resp));
phase = unwrap(angle(h_resp));
gd = -diff(phase)./diff(w); % 计算群延时

% 绘制完整的频率响应（幅度、相位和群延时）
figure(2)
subplot(3,1,1)
plot(w, mag_db)
grid on
title('幅度响应')
xlabel('频率 (Hz)')
ylabel('幅度 (dB)')
ylim([-60 5])

subplot(3,1,2)
plot(w, phase)
grid on
title('相位响应')
xlabel('频率 (Hz)')
ylabel('相位 (弧度)')

subplot(3,1,3)
plot(w(1:end-1), gd)
grid on
title('群延时')
xlabel('频率 (Hz)')
ylabel('群延时 (样本)')

% 3. 计算关键频率点的响应
[h_resp_check, w_check] = freqz(b, a, w_norm*pi);
mag_db_check = 20*log10(abs(h_resp_check));

% 输出系统特性分析
fprintf('\n系统特性分析：\n');
fprintf('1. 滤波器阶数：%d\n', N);
fprintf('2. 滤波器类型：椭圆带阻滤波器\n');
fprintf('3. 系统特点：\n');
fprintf(' - 通带等波纹响应\n');
fprintf(' - 阻带单调衰减\n');
fprintf(' - 非线性相位\n');
fprintf(' - 递归结构\n');
fprintf(' - 需要检查稳定性\n\n');

% 检查稳定性
poles = roots(a);
isStable = all(abs(poles) < 1);
fprintf('稳定性分析：\n');
if isStable
    fprintf('系统稳定（所有极点在单位圆内）\n');
else
    fprintf('系统不稳定（存在极点在单位圆外）\n');
end
fprintf('最大极点模值：%.4f\n\n', max(abs(poles)));

% 输出关键频率点的衰减值
fprintf('关键频率点的衰减值：\n');
fprintf('1.2 kHz: %.2f dB (要求 < 1 dB)\n', -mag_db_check(1));
fprintf('1.9 kHz: %.2f dB (要求 > 40 dB)\n', -mag_db_check(2));
fprintf('4.7 kHz: %.2f dB (要求 > 40 dB)\n', -mag_db_check(3));
fprintf('9.2 kHz: %.2f dB (要求 < 1 dB)\n', -mag_db_check(4));

% 计算并显示滤波器的具体特性
fprintf('\n滤波器性能指标：\n');
fprintf('1. 通带波纹：%.2f dB\n', Rp);
fprintf('2. 阻带衰减：%.2f dB\n', Rs);
fprintf('3. 滤波器阶数：%d\n', N);
fprintf('4. 极点数量：%d\n', length(poles));
fprintf('5. 零点数量：%d\n', length(roots(b)));

% 输出系数
fprintf('\n滤波器系数：\n');
fprintf('分子系数 b：\n');
disp(b);
fprintf('分母系数 a：\n');
disp(a);

% 保存滤波器系数
% save('iir_bandstop_filter.mat', 'b', 'a');