数字信号处理公式变程序（四）——巴特沃斯滤波器（下）
2024-05-29 13:15:07
之前已经讲过巴特沃斯滤波器的基础知识和数字滤波器求系统函数的代码实现，本节讲如何使用数字滤波器的系统函数实现对信号的滤波。
注：可能会有不足或者理解偏差的地方，路过的高人请不吝赐教。
OK，开始！
====================================================
一、实现filter方法对信号滤波
在理论讲解部分已经介绍过有关filter对信号滤波的相关知识，抄录如下：
滤波过程就是解常系数线性差分方程的过程，形式如：
其中，x(n)序列为滤波前的信号序列，ak, bm为H(z)系统函数分母与分子的系统数组，求出的y(n)即为滤波后的信号序列。
注：x(n)与y(n)的长度要相等，且a0=1。
公式的化简工程如下：
默认条件，当k<0时，x(k), y(k)都为0。例如n=0时，y(0)=b0*x0+b1*x(-1)+...+bM*x(-M)-a1*y(-1)-...-aN*y(-N)=b0*x(0)。经过迭代，可以求出y(n)序列的所有值。
filter的流程图如下所示
至此，滤波器设计的实现已经全部讲完，为了查看滤波效果，下面进行简单测试。
二、数字滤波的滤波测试
1.先来看看低通情况
选择的参数数值和波形点的值如下所示
//------------------------------------------------------------------------------------------// 1.低通滤波器//// 通带截止频率：passF = 2000Hz，阻带截止频率：stopF = 2500Hz，抽样频率：fs = 10000Hz// 通带衰减：rp = 2dB，阻带衰减：rs = 20dB//------------------------------------------------------------------------------------------float passF1 = 2000, stopF1 = 2500, fs1 = 10000, rp1 = 2, rs1 = 20; //数字滤波器性能float xVector1[100],yVector1[100];for(int i = 0; i < 100; i++){xVector1[i] = sin(2*M_PI*1000*i/10000);xVector1[i] += sin(2*M_PI*3000*i/10000);}代入设计的程序得到的结果为
将得到的数据复制到matlab中，绘制出滤波器的系统曲线和滤波前后的信号曲线，如下图所示：
2.来看看高通情况
测试过程与低通一致，这里就直接贴图了。参数、信号如下
//------------------------------------------------------------------------------------------// 2.高通滤波器//// 通带截止频率：passF = 3000Hz，阻带截止频率：stopF = 2800Hz，抽样频率：fs = 10000Hz// 通带衰减：rp = 3dB，阻带衰减：rs = 10dB//------------------------------------------------------------------------------------------float passF2 = 3000, stopF2 = 2800, fs2 = 10000, rp2 = 3, rs2 = 10; //数字滤波器性能float xVector2[100],yVector2[100];for(int i = 0; i < 100; i++){xVector2[i] = sin(2*M_PI*1000*i/1230);xVector2[i] += sin(2*M_PI*4000*i/1230);}输出结果
数据导入matlab绘图
3.来看看带通情况
直接贴图了。参数、信号如下
//------------------------------------------------------------------------------------------// 3.带通滤波器//// 通带截止频率：passF = [15000 40000]Hz，阻带截止频率：stopF = [10000 48000]Hz，抽样频率：fs = 100000Hz// 通带衰减：rp = 2dB，阻带衰减：rs = 30dB//------------------------------------------------------------------------------------------float passF3[2] = {15000.0, 40000.0}, stopF3[2] = {10000.0, 48000.0}, fs3 = 100000, rp3 = 2, rs3 = 30; //数字滤波器性能//float passF3[2] = {20000.0, 30000.0}, stopF3[2] = {10000.0, 45000.0}, fs3 = 100000, rp3 = 2, rs3 = 20; //数字滤波器性能float xVector3[100],yVector3[100];for(int i = 0; i < 100; i++){xVector3[i] = sin(2*M_PI*10000*i/12300);xVector3[i] += sin(2*M_PI*30000*i/12300);}输出结果
数据导入matlab绘图
4.来看看带阻情况
直接贴图了。参数、信号如下
//------------------------------------------------------------------------------------------// 4.带阻滤波器//// 通带截止频率：passF = [10000 45000]Hz，阻带截止频率：stopF = [20000 30000]Hz，抽样频率：fs = 100000Hz// 通带衰减：rp = 1dB，阻带衰减：rs = 20dB//------------------------------------------------------------------------------------------float passF4[2] = {10000.0, 45000.0}, stopF4[2] = {20000.0, 30000.0}, fs4 = 100000, rp4 = 1, rs4 = 20; //数字滤波器性能float xVector4[100],yVector4[100];for(int i = 0; i < 100; i++){xVector4[i] = sin(2*M_PI*10000*i/12300);xVector4[i] += sin(2*M_PI*25000*i/12300);}输出结果
数据导入matlab绘图
三、数字滤波的算法分析
1.从上图中的滤波后的信号对比可以看出iOS程序的滤波效果与MATLAB的滤波效果几乎一样，这说明功能已经实现；
2.但是！为毛系统函数曲线和相位曲线会有差异？原因是①iOS计算出的所有参数我保存的有效位数为8位，也就是精度与matlab有差异，导致了上述问题；②计算过程中中间值取值不一样。
四、测试用的MATLAB程序
%% 滤波器测试
close all;
clear,clc;
%% 低通滤波器
% 原始信号
t1 = 0:1/10000:99/10000;
signal1 = sin(2*pi*1000*t1) + sin(2*pi*3000*t1);% Matlab设计及滤波
Wp1 = 2000/5000; Ws1 = 2500/5000;Rp1 = 2; Rs1 = 20;
[n1,Wn1] = buttord(Wp1,Ws1,Rp1,Rs1);
[b11,a11] = butter(n1,Wn1,'low');
fprintf('\n\n');
disp(['阶数 n = ' num2str(n1) ' Wn = ' num2str(Wn1)]);
disp(['分子 b = ' num2str(b11)]);
disp(['分母 a = ' num2str(a11)]);
figure;
freqz(b11,a11,128,10000);
title('MATLAB设计低通');
signal11 = filter(b11, a11, signal1);% ios设计及滤波
% 分子、分母系数
b12 = [0.00149779 0.01348012 0.05392048 0.12581446 0.18872169 0.18872169 0.12581446 0.05392048 0.01348012 0.00149779];
a12 = [1.00000000 -1.44011279 2.06045886 -1.54849047 0.99939943 -0.41504301 0.13498033 -0.02782488 0.00372147 -0.00021985];
figure;
freqz(b12,a12,128,10000);
title('IOS设计低通');% IOS滤波后信号
signal12 = [0.000000 0.002305 0.024607 0.119104 0.344789 0.663375 0.891053 0.844175 0.507787 -0.007554 -0.564512 -0.962252 -0.976236 -0.567524 0.060922 0.641034 0.966314 0.923191 0.538677 -0.036528 -0.609295 -0.969491 -0.947205 -0.546920 0.050043 0.621007 0.964918 0.936902 0.545752 -0.043208 -0.617246 -0.968174 -0.940973 -0.545176 0.046354 0.618183 0.966335 0.939519 0.545864 -0.045024 -0.618104 -0.967247 -0.939966 -0.545405 0.045539 0.617997 0.966835 0.939866 0.545656 -0.045360 -0.618099 -0.967006 -0.939866 -0.545534 0.045413 0.618035 0.966942 0.939884 0.545588 -0.045402 -0.618069 -0.966963 -0.939870 -0.545566 0.045401 0.618053 0.966957 0.939878 0.545574 -0.045404 -0.618060 -0.966958 -0.939874 -0.545572 0.045402 0.618057 0.966959 0.939876 0.545572 -0.045403 -0.618058 -0.966958 -0.939875 -0.545572 0.045403 0.618058 0.966958 0.939876 0.545572 -0.045403 -0.618058 -0.966958 -0.939875 -0.545572 0.045403 0.618058 0.966958 0.939875 0.545572 -0.045403];% 滤波后信号绘图比较
figure;
subplot(211);
plot(t1,signal1);
title('低通原始信号');subplot(223);
plot(t1,signal11);
title('Matlab滤波后效果');subplot(224);
plot(t1,signal12);
title('IOS滤波后效果');%% 高通滤波器
% 原始信号
t2 = 0:1/1230:99/1230;
signal2 = sin(2*pi*1000*t2) + sin(2*pi*4000*t2);% Matlab设计及滤波
Wp2 = 3000/5000; Ws2 = 2800/5000;Rp2 = 3; Rs2 = 10;
[n2,Wn2] = buttord(Wp2,Ws2,Rp2,Rs2);
[b21,a21] = butter(n2,Wn2, 'high');
fprintf('\n\n');
disp(['阶数 n = ' num2str(n2) ' Wn = ' num2str(Wn2)]);
disp(['分子 b = ' num2str(b21)]);
disp(['分母 a = ' num2str(a21)]);
figure;
freqz(b21,a21,128,10000);
title('MATLAB设计高通');
signal21 = filter(b21, a21, signal2);% ios设计及滤波
% 分子、分母系数
b22 = [0.00111091 -0.00999820 0.03999278 -0.09331649 0.13997473 -0.13997473 0.09331649 -0.03999278 0.00999820 -0.00111091];
a22 = [1.00000000 1.74937261 2.46983593 2.04364540 1.31994556 0.58247111 0.19027367 0.04094845 0.00550463 0.00033600];
figure;
freqz(b22,a22,128,10000);
title('IOS设计高通');% IOS滤波后信号
signal22 = [0.000000 0.000086 -0.001742 0.012600 -0.047186 0.100503 -0.111185 0.013206 0.122525 -0.111056 -0.061493 0.143739 -0.001911 -0.130231 0.039947 0.098974 -0.051564 -0.070078 0.046668 0.052642 -0.036133 -0.047343 0.027154 0.050315 -0.022482 -0.056516 0.021789 0.061909 -0.023545 -0.064402 0.026367 0.063836 -0.029499 -0.061347 0.032642 0.058502 -0.035587 -0.056500 0.038019 0.055732 -0.039617 -0.055780 0.040302 0.055808 -0.040397 -0.055097 0.040538 0.053452 -0.041354 -0.051267 0.043111 0.049250 -0.045546 -0.047970 0.047998 0.047518 -0.049783 -0.047452 0.050592 0.047070 -0.050663 -0.045839 0.050650 0.043716 -0.051230 -0.041174 0.052706 0.038925 -0.054830 -0.037490 0.056949 0.036892 -0.058393 -0.036640 0.058874 0.036023 -0.058652 -0.034531 0.058387 0.032161 -0.058738 -0.029420 0.059975 0.027023 -0.061819 -0.025469 0.063604 0.024745 -0.064677 -0.024333 0.064783 0.023520 -0.064219 -0.021824 0.063655 0.019276 -0.063736 -0.016410 0.064699 0.013939];% 滤波后信号绘图比较
figure;
subplot(211);
plot(t1,signal2);
title('高通原始信号');subplot(223);
plot(t1,signal21);
title('Matlab滤波后效果');subplot(224);
plot(t1,signal22);
title('IOS滤波后效果');%% 带通滤波器
% 原始信号
t3 = 0:1/12300:99/12300;
signal3 = sin(2*pi*10000*t3) + sin(2*pi*30000*t3);% Matlab设计及滤波
Wp3 = [15000 40000]/50000; Ws3 = [10000 48000]/50000;Rp3 = 2; Rs3 = 30;
[n3,Wn3] = buttord(Wp3,Ws3,Rp3,Rs3);
[b31,a31] = butter(n3,Wn3);
fprintf('\n\n');
disp(['阶数 n = ' num2str(n3) ' Wn = ' num2str(Wn3)]);
disp(['分子 b = ' num2str(b31)]);
disp(['分母 a = ' num2str(a31)]);
figure;
freqz(b31,a31,128,100000);
title('Matlab设计带通');
signal31 = filter(b31, a31, signal3);% ios设计及滤波
% 分子、分母系数
b32 = [0.02089926 0.00000000 -0.14629485 0.00000000 0.43888454 0.00000000 -0.73147424 0.00000000 0.73147424 0.00000000 -0.43888454 0.00000000 0.14629485 0.00000000 -0.02089926];
a32 = [1.00000000 1.48232097 0.68685736 0.40477466 1.26650163 1.07464274 0.26942426 0.12877156 0.27672881 0.12952094 0.00562433 0.00391479 0.00953371 0.00144346 -0.00026785];
figure;
freqz(b32,a32,128,100000);
title('IOS设计带通');% IOS滤波后信号
signal32 = [0.000000 -0.011470 -0.012362 0.133377 0.020989 -0.504277 0.060238 0.856009 -0.034678 -0.646678 -0.533126 0.070799 1.098486 0.403407 -0.617302 -0.907387 -0.360178 0.992291 0.840365 -0.115532 -0.951133 -0.767798 0.556551 0.952045 0.361045 -0.738074 -0.935849 0.104443 0.893475 0.676192 -0.438216 -0.993646 -0.306054 0.734025 0.892950 -0.050160 -0.953524 -0.659203 0.424352 0.992014 0.365742 -0.757231 -0.890260 0.009727 0.939750 0.700832 -0.414879 -0.975872 -0.399614 0.733296 0.906831 0.005338 -0.918347 -0.721237 0.390575 0.980315 0.413980 -0.709504 -0.919479 -0.030340 0.913009 0.733389 -0.361708 -0.984191 -0.436449 0.692539 0.925505 0.060260 -0.905874 -0.748485 0.336629 0.982200 0.461034 -0.675117 -0.932113 -0.086892 0.895190 0.764759 -0.312205 -0.980479 -0.483469 0.655561 0.939733 0.112786 -0.884287 -0.779758 0.286378 0.979200 0.505085 -0.635462 -0.946374 -0.139267 0.873351 0.793849 -0.260238 -0.977091 -0.526652 0.615244 0.952188 0.165655];% 滤波后信号绘图比较
figure;
subplot(211);
plot(t1,signal3);
title('带通原始信号');subplot(223);
plot(t1,signal31);
title('Matlab滤波后效果');subplot(224);
plot(t1,signal32);
title('IOS滤波后效果');%% 带阻滤波器
% 原始信号
t4 = 0:1/12300:99/12300;
signal4 = sin(2*pi*10000*t4) + sin(2*pi*25000*t4);% Matlab设计及滤波
Wp4 = [10000 45000]/50000; Ws4 = [20000 30000]/50000;Rp4 = 1; Rs4 = 20;
[n4,Wn4] = buttord(Wp4,Ws4,Rp4,Rs4);
[b41,a41] = butter(n4,Wn4,'stop');
% b41 = [0.2691 0.0000 0.8074 0.0000 0.8074 0.0000 0.2691];
% a41 = [1.0000 0.0000 0.6451 0.0000 0.4439 0.0000 0.0640];
fprintf('\n\n');
disp(['阶数 n = ' num2str(n4) ' Wn = ' num2str(Wn4)]);
disp(['分子 b = ' num2str(b41)]);
disp(['分母 a = ' num2str(a41)]);
figure;
freqz(b41,a41,128,100000);
title('Matlab设计带阻');
signal41 = filter(b41, a41, signal4);% ios设计及滤波
% 分子、分母系数
b42 = [0.26912293 -0.00000000 0.80736878 -0.00000000 0.80736878 -0.00000000 0.26912293];
a42 = [1.00000000 -0.00000000 0.64508705 -0.00000000 0.44390629 -0.00000000 0.06399006];
figure;
freqz(b42,a42,128,100000);
title('IOS设计带阻');% IOS滤波后信号
signal42 = [0.000000 -0.193699 -0.084565 -0.200709 0.266234 0.737165 1.074895 1.223653 0.967087 0.713045 0.768339 0.963786 1.068077 0.873396 0.440608 0.125843 0.083439 0.113238 -0.046823 -0.435940 -0.805283 -0.902246 -0.788876 -0.746856 -0.920383 -1.144710 -1.143013 -0.866163 -0.558718 -0.459143 -0.512245 -0.446050 -0.113495 0.316975 0.557749 0.546198 0.511826 0.685090 1.006668 1.193786 1.082902 0.828061 0.710817 0.795579 0.853842 0.650082 0.237963 -0.097386 -0.180531 -0.151449 -0.283443 -0.637455 -0.977900 -1.053837 -0.900350 -0.790095 -0.894778 -1.072668 -1.043293 -0.734639 -0.376138 -0.222987 -0.250289 -0.192509 0.120599 0.544769 0.788483 0.761613 0.677457 0.783149 1.050018 1.204605 1.062408 0.754359 0.567767 0.597082 0.633452 0.428540 0.008133 -0.350135 -0.447985 -0.400058 -0.483280 -0.787519 -1.097171 -1.149493 -0.952550 -0.772293 -0.804780 -0.936327 -0.886269 -0.556772 -0.159678 0.035764 0.024635 0.062605 0.341711 0.744464 0.978831 0.928630];% 滤波后信号绘图比较
figure;
subplot(211);
plot(t1,signal4);
title('带阻原始信号');subplot(223);
plot(t1,signal41);
title('Matlab滤波后效果');subplot(224);
plot(t1,signal42);
title('IOS滤波后效果');==========================================================
OK，滤波器暂告一段落。
数字信号处理公式变程序（四）——巴特沃斯滤波器（下）相关推荐
- MATLAB_数字信号处理_模拟滤波器_设计巴特沃斯滤波器
简介 巴特沃斯滤波器-百度百科 巴特沃斯滤波器-维基百科 简介:巴特沃斯滤波器是一种模拟滤波器,它在频率响应方面具有特殊的属性.它被设计为具有均匀的幅度响应,即在通带内,它对所有频率的增益是相等的,而 ...
- 数字信号处理——巴特沃斯滤波器设计
设计思路 这里采用间接法设计数字滤波器(先设计模拟滤波器再设计数字滤波器) 滤波器理解: 1.数字滤波器可以用H(z),h(n)or系统差分方程来表示,对应的就是一个系统,信号输入该系统即可改变其所含 ...
- Matlab语音信号去噪程序，使用低通巴特沃斯滤波器
Matlab语音信号去噪程序,使用低通巴特沃斯滤波器. 1.读取一段歌曲的信号,绘制时域频域图,并播放. 2.添加正弦噪声: 3.设计巴特沃斯低通滤波器: 4.使用滤波器去除噪声,并画出时域频域图,播 ...
- 巴特沃斯滤波器应用场合_巴特沃斯数字低通滤波器设计及应用
龙源期刊网 http://www.qikan.com.cn 巴特沃斯数字低通滤波器设计及应用 作者:汪其锐 王桂华 王永军 来源:<山东工业技术> 2016 年第 24 期 摘 要:现实生 ...
- VLSI数字信号处理系统——第十四章冗余运算
VLSI数字信号处理系统--第十四章冗余运算 作者:夏风喃喃 参考: (1) VLSI数字信号处理系统:设计与实现 (美)Keshab K.Parhi/著 (2) socvista https://w ...
- 信号处理之巴特沃斯滤波器的理解----2022/11/30
巴特沃斯理解 问题展示 巴特沃斯滤波器概念理解 实现以及效果 结果展示 问题展示 在处理信号的时候发现对巴特沃斯理解得不到位,因此设计不出想要的滤波器.特此学习记录. 巴特沃斯滤波器概念理解 学习来源 ...
- 设计一个三阶巴特沃斯滤波器_巴特沃斯滤波器matlab实现
巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有起伏,而在阻频带则逐渐下降为零. 在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大. 一 ...
- matlab模拟巴特沃斯滤波器设计,巴特沃斯滤波器matlab实现
描述 巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有起伏,而在阻频带则逐渐下降为零. 在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大 ...
- 巴特沃斯滤波器 python_巴特沃斯、切比雪夫、贝塞尔滤波器的区别
巴特沃斯滤波器.切比雪夫滤波器.贝塞尔滤波器均包括模拟滤波器和数字滤波器两种形式. 数字滤波器是指完成信号滤波处理功能的,用有限精度算法实现的离散时间线性非时变系统,其输入是一组数字量,其输出是经过变 ...
最新文章
- linux下修改/dev/shm tmpfs文件系统大小
- cgi硬盘安装器_简简单单，玩转虚拟硬盘装多系统
- 理解模板引擎Razor 的原理
- 保存到redis的字符串类型出现斜杆_深入浅出Redis：这次从Redis底层数据结构开始...
- 实现机器学习的循序渐进指南I——KNN
- php取mysql某列的值,php – 获取MYSQL中某些列为null的表中的值
- java自学难点_学习JAVA遇到的难点总结
- Flume-NG源码阅读之Interceptor(原创)
- button3 电脑上mouse,鼠标侧键设置工具(X-Mouse Button Control)
- python实践日记二
- Google搜索的基本语法
- 快速删除excel中的空行
- centos xfs硬盘扩容
- Simulink基于level 2的s-function C语言编写
- 爬虫 第六讲 Scrapy框架
- 用chatgpt做ppt
- BZOJ 小约翰的游戏John 反尼姆博弈
- 搜索引擎技术 —— 检索模型
- iphone分辨率终极指南(含iphone6/6+)
- 一种字符编码猜测工具的实现方法