滤波器设计与信号处理、数字模拟滤波器的转换（完整整理版）

最新推荐文章于 2024-06-08 18:36:41 发布

擦擦擦大侠

最新推荐文章于 2024-06-08 18:36:41 发布

阅读量3.4k

点赞数 3

分类专栏：信号处理

本文链接：https://blog.csdn.net/chenshiming1995/article/details/109130940

版权

信号处理专栏收录该内容

8 篇文章 6 订阅

订阅专栏

滤波器设计与信号处理、数字模拟滤波器的转换

文章目录

滤波器设计与信号处理、数字模拟滤波器的转换

本文主要讨论数字信号处理、滤波器设计、模拟滤波器与数字滤波器的互相转换几个部分。

1 信号处理

1.1 双线性变换设计递归滤波器

如何设计数字滤波器¹：

试着写出双线性变换法设计IIR数字高通滤波器的主要步骤
- 将数字高通的频率指标转换为模拟高通的频率指标（其中将高通截止频率通过预畸转换为模拟高通的截止频率）
- 将模拟高通技术指标转换成为模拟低通技术指标
- 设计模拟低通原型滤波器
- 将模拟低通原型滤波器通过双线形映射为数字低通原型滤波器（下面一小节主要介绍这方面的工作）

（理解主要就是，将频率滤波的带宽从数字转换到模拟信号，再将相应的技术指标从模拟高通到模拟低通，然后设计一个相应的低通的模拟滤波器，再使用双线性映射到数字低通滤波器）

矩形窗、加窗对于滤波器性能的影响：

简述用窗函数法设计FIR滤波器的时，对理想低通滤波器加矩形窗处理后的影响；为了改善FIR的性能，窗函数要尽可能的满足什么样的条件？
- 幅频特性的陡直边沿被加宽，形成一个过渡带，其带宽取决于窗函数频响的主瓣宽度
- 过渡带的两侧附近产生起伏的尖峰和波纹，这是由于窗函数旁瓣引起的
- 增加截取长度N，将缩小窗函数主瓣宽度，但不能减小旁瓣的宽度值；能减小过渡带的宽度；能减小过渡带宽度，但不能改善带通内的平稳性和阻带衰减
- 为改善滤波器性能，窗函数需要满足：
  - 主瓣宽度尽量窄，以获得较陡的过渡带
  - 旁瓣值尽可能的小，以改善通带的平稳度和增大阻带的衰减。

1.2 相关数学内容

这里²详细介绍了矩形函数、三角函数与高斯函数的傅立叶变换的推导过程，作为参考。

下面是信号处理相关的一些记录。

可以用傅立叶变化的方法去得到滤波器或其他系统的幅频相应，如下。

H = fft(b,N)./fft(a,N);

1.3 使用matlab设计滤波器

1）设计过程

Matlab实现IIR数字滤波器设计
傅立叶变换、拉普拉斯变换、Z变换之间的联系？

20200628160951

2）给滤波器进行加窗

加窗过程是对FIR滤波器系数b矩阵进行计算：

w = kaiser(N,beta);
h = b.*w';% 设计的滤波器
mag = freqz(h,[1],omega);

b矩阵相当于“系统的单位冲激响应（时域上）”，而利用freqz函数可以转换到频域上观察幅度谱与相位谱，可以清晰看到滤波器的抖动有所改善。FIR的系数常常以sinc函数的形式呈现，其傅立叶变换之后在频域上的信号便为类似一个矩阵的形状。

2 模拟电路滤波器到数字滤波器

2.1 双线性变换原理部分³

下面解释了双线性变换方法的推导。

2.2 选取示例

1）双线性滤波器转换³

一个二阶低通滤波器可以通过以下无源电路实现

20200628164931

模拟系统的传递函数为：
$H(s)=\frac{Y(s)}{X(s)}=\frac{1}{L C s^{2}+R C s+1}$
双线性变换之后：
$H_{d}(z)=\frac{T^{2} z^{-2}+2 T^{2} z^{-1}+T^{2}}{\left(4 L C-2 T R C+T^{2}\right) z^{-2}+\left(-8 L C+2 T^{2}\right) z^{-1}+\left(4 L C-2 T R C+T^{2}\right)}$
可以看到模拟电路滤波器与数字电路的滤波器有如下的对应关系。此外，当采样频率越高（如下采样率为64100Hz时），其滤波器与模拟电路滤波器越接近。（但是这里的结果相反）

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-O51a1cpi-1602906396011)(/Users/chenshiming/Library/Application%20Support/typora-user-images/image-20200628195110471.png)]

下面给出matlab代码：

clear all;
close all;

fs = 44100;  %# sampling frequency
fc = 1000;  %# corner frequency of the lowpass 低通频率

% # coefficients of analog lowpass filter 模拟低通滤波器系数
Qinf = 0.8;
sinf = 2*pi*fc ;
C = 1e-6;
L = 1/(sinf^2*C);
R = sinf*L/Qinf;%% 这个应该是涉及到无源函数的性质
B = [0, 0, 1];
A = [L*C, R*C, 1];
% # cofficients of digital filter
T = 1/fs;
b = [T^2, 2*T^2, T^2];
a = [(4*L*C+2*T*R*C+T^2), (-8*L*C+2*T^2), (4*L*C-2*T*R*C+T^2)];

% # compute frequency responses
[ Hd , Om ] = freqz(b, a, 1024,fs);%#离散的
[ H , tmp ] = freqs(B, A);%# 连续的模拟信号
figure;
semilogx(tmp/pi/2 , 20*log10(H),'LineWidth',2);
hold on; semilogx(Om , 20*log10(Hd),'LineWidth',2);
set(gca,'Fontname','Times New Roman','fontsize',16);
xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')

2）matlab的双线性变化函数表示

使用双线性变换将连续滤波器变成离散滤波器。（这个非常重要）

%% 这是一个非常重要的实例
fs = 1000;
T = 1/fs;
b = [1 2 3];
a = [2 1 3];
figure;bode(b,a);
[dnum,dden]= bilinear(b,a,fs);
[ H , f ] = freqz(dnum,dden);
[h1 , Om] = freqs(b,a);
figure;
semilogx(f*fs/2/pi, 20*log10(abs(H)),'LineWidth',1);
hold on;
semilogx(Om/2/pi , 20*log10(abs(h1)),'LineWidth',2);
legend 数字滤波 模拟滤波

20200730104904

或者：（这与freqz函数的设定相关）

%% 这是一个非常重要的实例
fs = 1000;
T = 1/fs;
b = [1 2 3];
a = [2 1 3];
figure;bode(b,a);
[dnum,dden]= bilinear(b,a,fs);
[ H , f ] = freqz(dnum,dden,1024,fs);
[h1 , Om] = freqs(b,a);
figure;
semilogx(f, 20*log10(abs(H)),'LineWidth',1);
hold on;
semilogx(Om/2/pi , 20*log10(abs(h1)),'LineWidth',2);
legend 数字滤波 模拟滤波

3）freqz函数与freqs函数

所有的freqz的直接输出的f都是从0～pi的。freqz直接画出来的图都是norminialized。

freqz(b,a)

20200730111004

下面给出了freqs函数与传递函数之间的关系。
$H(s)=\frac{B(s)}{A(s)}=\frac{b(1) s^{n}+b(2) s^{n-1}+\cdots+b(n+1)}{a(1) s^{m}+a(2) s^{m-1}+\cdots+a(m+1)}$

下面给出了freqz函数与传递函数之间的关系。

            jw               -jw              -jmw
     jw  B(e)    b(1) + b(2)e + .... + b(m+1)e
  H(e) = ---- = ------------------------------------
            jw               -jw              -jnw
         A(e)    a(1) + a(2)e + .... + a(n+1)e

这里有一个问题就是：这个结果与bode图的结果并不一致。bode与freqs的横坐标不需要除以2*pi的。

下面解释了这个结果：

clear;
%一阶模拟抗混滤波器
omiga1 = 10^5 / 2^17;
f  = 0.001 : 0.0001 :100;  % 模拟频率Hz
Omega =  2*pi.*f;         % 模拟角频率 rad/sec
H = omiga1 ./ (1j.*(Omega) + omiga1);
figure
semilogx( f , 20*log10(H) );
xlabel('模拟频率 (Hz)');
ylabel('dB');            % 截止频率为 0.76/2/pi = 0.121Hz
title('一阶抗混滤波器');
grid on
hold on;
[h,f] = freqs(omiga1,[1,omiga1]);
semilogx(f/pi/2 , 20*log10(h),'LineWidth',2);
grid on;

那么Freqz函数是如何计算滤波器的频率响应的呢？下面给出了解释：

psi = 0:0.001:900;% 不能超过fs，其实就是频率
fs = 1000;
wd = 0.001;
Rz = (1-wd)*(1- exp(-j*2*pi*psi/fs))./( 1-(1-wd)*exp(-j*2*pi*psi/fs) );
figure;semilogx(psi,20*log10(Rz),'LineWidth',1);
[h,f] = freqz(b,a,1e5,fs);%1e5越大，那么画出来的分辨率越高
hold on;semilogx(f,20*log10(h),'--r','LineWidth',2)

20200923161740

完全重合，说明这种方式画图的准确性

2.3 反向差分法

反向差分法：参考连续系统离散化方法.pdf（还有其他的变换方法）。离散滤波器的过程特性及频率特性同原连续滤波器比较有一定的失真，需要较小的采样周期T。

$s=\frac{1-z^{-1}}{T}$

2.4 讨论滤波器的设计

1） fdatool设计

fdatool:

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-PiXhcXGw-1602906396015)(pic/image-20200914100014027.png)]

b = Num;
a = [1];
fs = 48000;
[h,f] = freqz(b,a,1024,fs);
figure;plot(f/1000, 20*log10(abs(h)),'LineWidth',1);

2) 系统频率响应的画图

b = [1,2,3];
a = [2,3,4];
fs = 1000;
N = 512;
H = fft(b,N)./fft(a,N);   % H矩阵
mag = 20*log10(abs(H));    % get magnitude of spectrum in dB 幅值
phase = angle(H)*2*pi;     % get phase in deg.相位
faxis = fs/2*linspace(0,1,N/2);  % the axis of frequency

figure,
subplot(2,1,1),plot(faxis,mag(1:N/2))
xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
grid on
subplot(2,1,2),plot(faxis,phase(1:N/2),'r')
xlabel('Frequency (Hz)'),ylabel('Phase (deg.)')
grid on
figure;freqz(b,a)

参考文献

博客 - 数据处理与滤波器性质 ↩︎
https://blog.csdn.net/qq_22943397/article/details/80301398 ↩︎
参考开源项目digital-signal-processing-lecture-master ↩︎ ↩︎

擦擦擦大侠

关注

3
点赞
踩
24

收藏

觉得还不错? 一键收藏
打赏
0
评论
滤波器设计与信号处理、数字模拟滤波器的转换（完整整理版）

滤波器设计与信号处理、数字模拟滤波器的转换文章目录滤波器设计与信号处理、数字模拟滤波器的转换1 信号处理1.1 双线性变换设计递归滤波器1.2 相关数学内容1.3 使用matlab设计滤波器1）设计过程2）给滤波器进行加窗2 模拟电路滤波器到数字滤波器2.1 双线性变换原理部分[^3]2.2 选取示例1）双线性滤波器转换[^3]2）matlab的双线性变化函数表示3）freqz函数与freqs函数2.3 反向差分法2.4 讨论滤波器的设计1） fdatool设计2) 系统频率响应的画图参考文献本文主要讨
复制链接

扫一扫