数字滤波器之高通滤波器设计

Posted 2023-04-02 Ivy9808

要求：设计一款高通滤波器，用在音频信号处理过程中，滤掉100Hz以下的信号。

实现方法：通过在Z平面放置零极点的来设计数字滤波器。
设计要点：单位圆上对应于需要加强频率的点附近放置极点，即w=π处附近落极点，需要拉低的频点附近放置零点，即w=0处附近落零点。为了滤波器的稳定，所有极点必须在单位圆内，而零点可以在任何位置。为了使滤波器系数是实数，所有的复数值的极点和零点必须以共轭的形式出现。
设计过程中需解决的问题：

1. 模拟信号频率变量F（或Ω）和离散时间信号的频率变量f（或ω）之间的关系为：f = F/Fs ; 其中Fs为采样频率。因此若设计截止频率为100Hz的高通滤波器，采样频率16KHz时，则其在Z平面的对应值为π/80≈0.04。

2. 简单高通滤波器的零极点位置有以下几种形式:(1)单极点单零点（2）双极点单零点（3）多极点多零点。为了降低滤波器频率特性对误差的敏感程度，不采用单极点、单零点的方式，因此选用方式3进行高通滤波器的设计。
   
   
   

3. 滤波器在Z域由零极点表示的表达式如下：
   
   其中，z=e^jw
   高通滤波器的幅频曲线中，w=0处，H(w)=0；w=π处，H(w)=1;
   初步假设高通滤波器所具有三个极点和两个零点：
   z1=0；z2=1；p1=-√3/2；p2=-√3/2+j1/2;p2=-√3/2-j1/2。
   将各零极点值带入到H（z）表达式中，并根据w=π时H（w）=1，可得出b0的值，再将H（z）写成多项式的形式，如下：
   

此时，该多项式对应的时域差分方程为：

至此，该高通滤波器设计初步完成，但是并未使用到截止频率的相关信息，并且需使用matlab进行频域信息的验证，后续再进行探究。待验证完成后，便可添加到音频处理过程中，进行低频信号的过滤。

[Matlab]双线性变换法设计数字高通滤波器

测试代码：

%%****bin_hp.m*******************%%
%% 使用双线性变换法设计高通滤波器
%% 2018年6月13日 14:27:37
%% author:Alimy

close all;
clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%代码正文
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%给定数字滤波器指标
f_p     =   200     ;   %阻带上截止频率
f_st    =   210     ;   %通带下截止频率
R_p     =   3       ;   %通带允许的最大衰减
R_st    =   30      ;   %阻带允许的最小衰减
f_s     =   1000    ;   %采样频率
T_s     =   1 / f_s ;   %采样间隔
%1.将数字高通滤波器的频率参数变换为归一化的数字角频率参数
omega_p  = 2 * pi * f_p  / f_s;
omega_st = 2 * pi * f_st / f_s;
%2.预畸变处理,将归一化数字角频率参数变换成模拟高通滤波器的角频率参数
C = 2 / T_s ;
Omega_p   = C * tan( omega_p  / 2 );
Omega_st  = C * tan( omega_st / 2 ); 
%3.对模拟高通滤波器的角频率参数做归一化处理
lamda_p  = Omega_p  / Omega_p;
lamda_st = Omega_st / Omega_p;
%4.设计归一化模拟滤波器,得到归一化模拟高通滤波器的转移函数
[ N , Wn ] = buttord( Omega_p , Omega_st , R_p , R_st , ‘s‘ ); %选择模拟巴特沃斯滤波器的最小阶数
[ z , p , k ] = buttap(N); %创建巴特沃斯模拟低通滤波器
[ Bp , Ap ] = zp2tf(z,p,k); %由零点、极点、增益确定传输函数的分子与分母系数
%5.利用归一化模拟高通滤波器的转移函数确定模拟高通滤波器的转移函数
[ b , a ] = lp2hp(Bp,Ap,Wn);
%6.利用模拟高通滤波器的转移函数确定IIR数字滤波器的转移函数 
[ bz , az ] = bilinear(b,a,f_s);

figure(1);
freqz(bz,az);
title(‘高通滤波器幅度谱和相位谱特性‘);

%滤波效果测试 
N = 1000
t = [ 0 : N - 1 ] * T_s ;
f1 = 20;
f2 = 100;
f3 = 300;
x1 = 2*1*sin( 2 * pi * f1 * t );
x2 = 2*2*sin( 2 * pi * f2 * t );
x3 = 2*1*sin( 2 * pi * f3 * t );
x  = x1 + x2 + x3;

fft_x = fft( x );
X_mag = fftshift( abs   ( fft_x ) ) / N ;
X_ang = fftshift( angle ( fft_x ) );
delta_f = f_s/N;
f = ( -N / 2 : N / 2 - 1 )*delta_f;

figure(2);
subplot(3,1,1);
plot(t,x);
title(‘原信号时域波形‘);
xlabel(‘t(s)‘);
subplot(3,1,2);
plot( f , X_mag );
title(‘原信号幅度谱‘);
xlabel(‘f(hz)‘);
subplot(3,1,3);
plot( f , X_ang );
title(‘原信号相位谱‘);
xlabel(‘f(hz)‘);


%滤波
lp_x = filter( bz , az , x );
lp_fft_x = fft( lp_x );
lp_X_mag = fftshift( abs   ( lp_fft_x ) ) / N ;
lp_X_ang = fftshift( angle ( lp_fft_x ) );
figure(3);
subplot(3,1,1);
plot(t,lp_x);
title(‘滤波后信号时域波形‘);
xlabel(‘t(s)‘);
subplot(3,1,2);
plot( f , lp_X_mag );
title(‘滤波后信号幅度谱‘);
xlabel(‘f(hz)‘);
subplot(3,1,3);
plot( f , lp_X_ang );
title(‘滤波后信号相位谱‘);
xlabel(‘f(hz)‘);

　　

　效果：

滤波器特性：

待滤波的信号：

滤波后信号：