数字信号处理实验实验三：数字滤波器的设计与实现

一．实验目的

1、理解滤波器参数的意义。

2、掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的方法。

3、掌握窗函数法设计FIR数字滤波器的方法。

4、掌握利用MATLAB完成各型数字滤波器设计的方法。

5、掌握分析滤波器是否达到性能指标的方法。

二．实验内容

    1、

    2、

    3、

    4、N=33 选择汉明窗

三、实验代码

1、

clc;
clear;

% 滤波器参数
f_p = 200;     % 截止频率
delta_p = 3;   % 最大衰减
f_s = 350;     % 起始频率
delta_s = 35;  % 最小衰减  
T = 0.001;     

w_p = 2 * pi * f_p;
w_s = 2 * pi * f_s;

W_p = w_p / (1/T);
W_s = w_s / (1/T);

[n, Wn] = buttord(W_p, W_s, delta_p, delta_s, 's');
[b, a] = butter(n, Wn, 's');


freq = 0:0.1:f_s;
w=2 * pi * freq *T;
[h,w] = freqs(b, a, w);

% 幅频特性曲线
figure;
plot(freq, 20*log10(abs(h)));
title('幅频特性');
xlabel('(Hz)');
ylabel('(dB)');
grid on;

2、

clc;clear;

fp = 200;   
delta_p = 3;    
fs = 350;   
delta_s = 35;   
T1 = 0.001; 
T2 = 0.0001; 

% 模拟滤波器的通带和阻带截止频率
wp = 2 * pi * fp;  
ws = 2 * pi * fs; 

[n, wn] = buttord(wp, ws, delta_p, delta_s, 's');
[b, a] = butter(n, wn, 's');

%%脉冲响应不变法
% 数字滤波器
[bz1, az1] = impinvar(b, a, 1/T1);
[bz2, az2] = impinvar(b, a, 1/T2);

% 幅频特性
freq_range = 0:0.01:pi;  
w=freq_range/pi;

h1 = freqz(bz1, az1, freq_range);  
h2 = freqz(bz2, az2, freq_range); 
h_db1 = 20 * log10(abs(h1));
h_db2 = 20 * log10(abs(h2));


%%双线性变换法
% 数字滤波器
[bz3, az3] = bilinear(b, a, 1/T1);
[bz4, az4] = bilinear(b, a, 1/T2);
% 幅频特性
 
h3 = freqz(bz3, az3, freq_range);
h4 = freqz(bz4, az4, freq_range);
h_db3 = 20 * log10(abs(h3));
h_db4 = 20 * log10(abs(h4));



% 幅频特性曲线
figure;
subplot(2,2,1);
plot(w, h_db1);
xlabel('(\pi rad)');
ylabel('(dB)');
title('脉冲响应不变法幅频特性T1');
grid on;

subplot(2,2,2);
plot(w, h_db2);
xlabel('(\pi rad)');
ylabel('(dB)');
title('脉冲响应不变法幅频特性T2');
grid on;

subplot(2, 2, 3);
plot(w, h_db3);
xlabel('(\pi rad)');
ylabel('(dB)');
title('双线性变换法幅频特性T1');
grid on;

subplot(2, 2, 4);
plot(w, h_db4);
xlabel('(\pi rad)');
ylabel('(dB)');
title('双线性变换法幅频特性T2');
grid on;

3、

clc;clear;

fs = 40;    
delta_s = 60;   
fp = 150;   
delta_p = 1;  
FT = 20 * fp;   

%通带和阻带频率
wp = fp / (FT/2);   
ws = fs / (FT/2);   

[n, wn] = cheb1ord(wp, ws, delta_p, delta_s);
[b, a] = cheby1(n, delta_p, wn, 'high');

% 幅频特性曲线
[h, w] = freqz(b, a);
h_db = 20 * log10(abs(h));

figure;
subplot(2, 1, 1);
plot(w/pi, h_db);
xlabel('频率 (Hz)');
ylabel('幅度响应 (dB)');
title('幅频特性');
grid on;

% 绘制滤波器的零极点图
subplot(2, 1, 2);
zplane(b, a);
title('零极点图');

4、

clc;clear;
% 滤波器参数
f_p = 1000;  
f_s = 2000;  
delta_s = 50;  
T = 0.1e-3; 
F_T = 1 / T;  
w_p = 2 * pi * f_p / F_T;  

% 窗函数和窗宽
N = 51;
rect_window = rectwin(N);
triang_window = triang(N);
hann_window = hann(N);
hamming_window = hamming(N);
blackman_window = blackman(N);

% 时域窗体图形
figure;
subplot(1, 1, 1);
stem(rect_window, 'b', 'DisplayName', '矩形窗');
hold on;
stem(triang_window, 'r', 'DisplayName', '三角窗');
stem(hann_window, 'g', 'DisplayName', '汉宁窗');
stem(hamming_window, 'm', 'DisplayName', '汉明窗');
stem(blackman_window, 'c', 'DisplayName', '布莱克曼窗');
title('时域窗体图形');
xlabel('样本');
ylabel('幅值');
legend('show');

% 定义滤波器的参数
fp = 1000;   
fs = 2000;   
T = 0.1e-3;  

% 计算归一化的通带和阻带频率
wp = fp * T;
ws = fs * T;
Bt=ws-wp; %过渡带带宽
N0=ceil(3.3/Bt)

% 汉明窗
b = fir1(N0, ws, hamming(N0+1));

% 频率响应
[h, w] = freqz(b, 1, 512);

h_db = 20 * log10(abs(h));
phase = unwrap(angle(h));

% 幅频特性和相频特性曲线
figure;
subplot(2, 1, 1);
plot(w/(2*pi*T), h_db);
xlabel('频率 (Hz)');
ylabel('幅度响应 (dB)');
title('幅频特性');
grid on;

subplot(2, 1, 2);
plot(w/(2*pi*T), phase);
xlabel('频率 (Hz)');
ylabel('相位响应 (弧度)');
title('相频特性');
grid on;