数字信号处理公式变程序(四)——巴特沃斯滤波器(下)
来源:互联网 发布:aaa云主机怎么绑定域名 编辑:程序博客网 时间:2024/06/06 18:44
之前已经讲过巴特沃斯滤波器的基础知识和数字滤波器求系统函数的代码实现,本节讲如何使用数字滤波器的系统函数实现对信号的滤波。
注:我比较喜欢使用matlab(也喜欢自己修改matlab的源码做测试,所以重装了好几次了,囧。。。),有部分参考matlab的算法进行(或者直接翻译为OC),算法的运行结果经常会跟matlab运算结果进行比较,算法只做学习用,转载请注明。
另外可能会有不足或者理解偏差的地方,路过的高人请不吝赐教。
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的流程图如下所示代码如下:
#pragma mark - 滤波方法/*====================================================================== * 方法名: filter * 方法功能:根据数字滤波器系统方法(系数),对原始信号进行滤波 * * 变量名称: * butterValue - 存放滤波器参数(阶数和截止频率)的结构体变量 * numerator - 系统方法,分子系数数组 * denominator - 系统方法,分母系数数组 * xVector - 输入的原始信号(数组) * length - 原始信号的长度,也是滤波后信号的长度 * yVector - 滤波后的信号(数组) * * 返回值: 设计是否成功,true-成功,false-失败 *=====================================================================*/+ (Boolean)filter:(ButterFilterStruct)butterValue andNumerator:(float *)numerator andDenominator:(float *)denominator andXVector:(float *)xVector andXLength:(int)length andReturnY:(float *)yVector{ Boolean isFilterOK = false; if(!butterValue.isFOK) { NSLog(@"系统方法错误!"); isFilterOK = false; return isFilterOK; } if(butterValue.N > 10) { NSLog(@"失败!滤波器的阶数不能大于10。"); isFilterOK = false; return isFilterOK; } int N = butterValue.length; //系数数组的长度 //返回值初始化 for(int i = 0; i < length; i++) { yVector[i] = 0.0; //后面循环中用到y递归算法,需要提前初始化 } //第一层循环,计算length个y的输出值 for(int i = 0; i < length; i++) { if(i == 0) { yVector[i] = numerator[i]*xVector[i]; } else { yVector[i] = numerator[0]*xVector[i]; //第二层循环,计算每个y的每一项 for(int j = 1; j <= i && j < N; j++) { yVector[i] += numerator[j]*xVector[i-j] - denominator[j]*yVector[i-j]; } } yVector[i] /= denominator[0]; } isFilterOK = true; return isFilterOK;}至此,滤波器设计的实现已经全部讲完,为了查看滤波效果,下面进行简单测试。
二、数字滤波的滤波测试
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的spectrogram函数(STFT)
- 数字信号处理公式变程序(一)——DFT、FFT
- 数字信号处理公式变程序(二)——插值
- 数字信号处理公式变程序(三)——压缩算法
- 数字信号处理——中值滤波器
- [数字信号处理]IIR滤波器的间接设计(C代码)
- [数字信号处理]IIR滤波器的间接设计(C代码)
- [数字信号处理] FIR滤波器基础
- [数字信号处理]IIR滤波器基础
- [数字信号处理] FIR滤波器基础
- [数字信号处理]IIR滤波器基础
- KCF(核化相关滤波)跟踪公式推导笔记(1)——线性情况下滤波器的解
- 数字信号处理集中滤波器的实现
- 数字信号处理——卷积
- 数字信号处理 DSP(一)
- 纯python实现的web: tornado性能测试以及实际使用解析
- uart过来数据不够,组包
- 软件开发过程学习笔记(六)之测试报告模板
- 单例模式
- Oracle高版本dmp导入到低版本
- 数字信号处理公式变程序(四)——巴特沃斯滤波器(下)
- abstract factory(抽象工厂) 对象创建型模式
- 类String的构造函数、析构函数和赋值函数
- 70天复习一次通过信息系统项目管理师考试经验和心得
- 句柄
- 用git更新线上项目代码后回滚到之前的稳定版本
- Elasticsearch scoring detailed explanation
- 论文提要“You Only Look Once: Unified, Real-Time Object Detection”
- 【Oracle】dbconsole exited with retCode 2. Server2003