数字信号处理公式变程序（四）——巴特沃斯滤波器（下）
2015-07-08 14:50
645 查看
之前已经讲过巴特沃斯滤波器的基础知识和数字滤波器求系统函数的代码实现，本节讲如何使用数字滤波器的系统函数实现对信号的滤波。
注：我比较喜欢使用matlab（也喜欢自己修改matlab的源码做测试，所以重装了好几次了，囧。。。），[u]有部分参考matlab的算法进行（或者直接翻译为OC），算法的运行结果经常会跟matlab运算结果进行比较，算法只做学习用，转载请注明。[/u]
另外可能会有不足或者理解偏差的地方，路过的高人请不吝赐教。
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.先来看看低通情况
选择的参数数值和波形点的值如下所示
代入设计的程序得到的结果为
将得到的数据复制到matlab中，绘制出滤波器的系统曲线和滤波前后的信号曲线，如下图所示：
2.来看看高通情况
测试过程与低通一致，这里就直接贴图了。参数、信号如下
输出结果
数据导入matlab绘图
3.来看看带通情况
直接贴图了。参数、信号如下
数据导入matlab绘图
4.来看看带阻情况
直接贴图了。参数、信号如下
数据导入matlab绘图
三、数字滤波的算法分析
1.从上图中的滤波后的信号对比可以看出iOS程序的滤波效果与MATLAB的滤波效果几乎一样，这说明功能已经实现；
2.但是！为毛系统函数曲线和相位曲线会有差异？原因是①iOS计算出的所有参数我保存的有效位数为8位，也就是精度与matlab有差异，导致了上述问题；②计算过程中中间值取值不一样。
四、测试用的MATLAB程序
==========================================================
OK，滤波器暂告一段落。
注：我比较喜欢使用matlab（也喜欢自己修改matlab的源码做测试，所以重装了好几次了，囧。。。），[u]有部分参考matlab的算法进行（或者直接翻译为OC），算法的运行结果经常会跟matlab运算结果进行比较，算法只做学习用，转载请注明。[/u]
另外可能会有不足或者理解偏差的地方，路过的高人请不吝赐教。
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，滤波器暂告一段落。
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 软件开发过程学习笔记（六）之测试报告模板
- windows 2012 r2下安装sharepoint 2013错误解决
- 软件开发过程学习笔记（六）之测试报告模板 分类： 开发过程 2015-07-08 14:50 8人阅读 评论(0) 收藏
- 判断一个数是否为2的若干次幂
- 实现编辑和删除
- Linux下重启mysql数据库的方法
- ipcs 使用
- C++顺序容器vector、deque、list
- 纯python实现的web: tornado性能测试以及实际使用解析
- linux中的sleep
- .NET 验证码/验证图片
- 计算一个整数中含1的个数
- 通过抽象工厂模式再理解java的上转型
- 大数据平台架构及数据分析平台架构
- eclipse自动补全的设置
- TCP Retransmission
- dp与px
- 高通8x12平台 mipi屏 调试
- 软件开发过程学习笔记（五）之测试方案模板
- C++管理指针成员