[数字信号处理]IIR滤波器的间接设计(C代码)(转)
PS:转自http://blog.csdn.net/zhoufan900428/article/details/9069475
1.模拟滤波器的设计
1.1巴特沃斯滤波器的次数
1.2巴特沃斯滤波器的传递函数
根据S解开,可以得到极点。这里,为了方便处理,我们分为两种情况去解这个方程。当N为偶数的时候,
这里,使用了欧拉公式。同样的,当N为奇数的时候,
同样的,这里也使用了欧拉公式。归纳以上,极点的解为
上式所求得的极点,是在s平面内,在半径为Ωc的圆上等间距的点,其数量为2N个。为了使得其IIR滤波器稳定,那么,只能选取极点在S平面左半平面的点。选定了稳定的极点之后,其模拟滤波器的传递函数就可由下式求得。
1.3巴特沃斯滤波器的实现(C语言)
其对应的C语言程序为
- N = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /
- log10 (Stopband/Cotoff) ));
然后是极点的选择,这里由于涉及到复数的操作,我们就声明一个复数结构体就可以了。最重要的是,极点的计算含有自然指数函数,这点对于计算机来讲,不是太方便,所以,我们将其替换为三角函数,
这样的话,实部与虚部就还可以分开来计算。其代码实现为
- typedef struct
- {
- double Real_part;
- double Imag_Part;
- } COMPLEX;
- COMPLEX poles[N];
- for(k = 0;k <= ((2*N)-1) ; k++)
- {
- if(Cotoff*cos((k+dk)*(pi/N)) < 0)
- {
- poles[count].Real_part = -Cotoff*cos((k+dk)*(pi/N));
- poles[count].Imag_Part= -Cotoff*sin((k+dk)*(pi/N));
- count++;
- if (count == N) break;
- }
- }
计算出稳定的极点之后,就可以进行传递函数的计算了。传递的函数的计算,就像下式一样
这里,为了得到模拟滤波器的系数,需要将分母乘开。很显然,这里的极点不一定是整数,或者来说,这里的乘开需要做复数运算。其复数的乘法代码如下,
- int Complex_Multiple(COMPLEX a,COMPLEX b,
- double *Res_Real,double *Res_Imag)
- {
- *(Res_Real) = (a.Real_part)*(b.Real_part) - (a.Imag_Part)*(b.Imag_Part);
- *(Res_Imag)= (a.Imag_Part)*(b.Real_part) + (a.Real_part)*(b.Imag_Part);
- return (int)1;
- }
有了乘法代码之后,我们现在简单的情况下,看看其如何计算其滤波器系数。我们做如下假设
这个时候,其传递函数为
将其乘开,其大致的关系就像下图所示一样。
计算的关系一目了然,这样的话,实现就简单多了。高阶的情况下也一样,重复这种计算就可以了。其代码为
- Res[0].Real_part = poles[0].Real_part;
- Res[0].Imag_Part= poles[0].Imag_Part;
- Res[1].Real_part = 1;
- Res[1].Imag_Part= 0;
- for(count_1 = 0;count_1 < N-1;count_1++)
- {
- for(count = 0;count <= count_1 + 2;count++)
- {
- if(0 == count)
- {
- Complex_Multiple(Res[count], poles[count_1+1],
- &(Res_Save[count].Real_part),
- &(Res_Save[count].Imag_Part));
- }
- else if((count_1 + 2) == count)
- {
- Res_Save[count].Real_part += Res[count - 1].Real_part;
- Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;
- }
- else
- {
- Complex_Multiple(Res[count], poles[count_1+1],
- &(Res_Save[count].Real_part),
- &(Res_Save[count].Imag_Part));
- 1 Res_Save[count].Real_part += Res[count - 1].Real_part;
- Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;
- }
- }
- *(b+N) = *(a+N);
到此,我们就可以得到一个模拟滤波器巴特沃斯低通滤波器了。
2.双1次z变换
2.1双1次z变换的原理
2.2双1次z变换的实现(C语言)
我们将其进行双1次z变换,我们可以得到如下式子
可以看出,我们还是需要将式子乘开,进行合并同类项,这个跟之前说的算法相差不大。其代码为。
- for(Count = 0;Count<=N;Count++)
- {
- for(Count_Z = 0;Count_Z <= N;Count_Z++)
- {
- Res[Count_Z] = 0;
- Res_Save[Count_Z] = 0;
- }
- Res_Save [0] = 1;
- for(Count_1 = 0; Count_1 < N-Count;Count_1++)
- {
- for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
- {
- if(Count_2 == 0) Res[Count_2] += Res_Save[Count_2];
- else if((Count_2 == (Count_1+1))&&(Count_1 != 0))
- Res[Count_2] += -Res_Save[Count_2 - 1];
- else Res[Count_2] += Res_Save[Count_2] - Res_Save[Count_2 - 1];
- for(Count_Z = 0;Count_Z<= N;Count_Z++)
- {
- Res_Save[Count_Z] = Res[Count_Z] ;
- Res[Count_Z] = 0;
- }
- }
- for(Count_1 = (N-Count); Count_1 < N;Count_1++)
- {
- for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
- {
- if(Count_2 == 0) Res[Count_2] += Res_Save[Count_2];
- else if((Count_2 == (Count_1+1))&&(Count_1 != 0))
- Res[Count_2] += Res_Save[Count_2 - 1];
- else
- Res[Count_2] += Res_Save[Count_2] + Res_Save[Count_2 - 1];
- }
- for(Count_Z = 0;Count_Z<= N;Count_Z++)
- {
- Res_Save[Count_Z] = Res[Count_Z] ;
- Res[Count_Z] = 0;
- }
- }
- for(Count_Z = 0;Count_Z<= N;Count_Z++)
- {
- *(az+Count_Z) += pow(2,N-Count) * (*(as+Count)) *
- Res_Save[Count_Z];
- *(bz+Count_Z) += (*(bs+Count)) * Res_Save[Count_Z];
- }
- }
到此,我们就已经实现了一个数字滤波器。
3.IIR滤波器的间接设计代码(C语言)
- #include <stdio.h>
- #include <math.h>
- #include <malloc.h>
- #include <string.h>
- #define pi ((double)3.1415926)
- struct DESIGN_SPECIFICATION
- {
- double Cotoff;
- double Stopband;
- double Stopband_attenuation;
- };
- typedef struct
- {
- double Real_part;
- double Imag_Part;
- } COMPLEX;
- int Ceil(double input)
- {
- if(input != (int)input) return ((int)input) +1;
- else return ((int)input);
- }
- int Complex_Multiple(COMPLEX a,COMPLEX b
- ,double *Res_Real,double *Res_Imag)
- {
- *(Res_Real) = (a.Real_part)*(b.Real_part) - (a.Imag_Part)*(b.Imag_Part);
- *(Res_Imag)= (a.Imag_Part)*(b.Real_part) + (a.Real_part)*(b.Imag_Part);
- return (int)1;
- }
- int Buttord(double Cotoff,
- double Stopband,
- double Stopband_attenuation)
- {
- int N;
- printf("Wc = %lf [rad/sec] \n" ,Cotoff);
- printf("Ws = %lf [rad/sec] \n" ,Stopband);
- printf("As = %lf [dB] \n" ,Stopband_attenuation);
- printf("--------------------------------------------------------\n" );
- N = Ceil(0.5*( log10 ( pow (10, Stopband_attenuation/10) - 1) /
- log10 (Stopband/Cotoff) ));
- return (int)N;
- }
- int Butter(int N, double Cotoff,
- double *a,
- double *b)
- {
- double dk = 0;
- int k = 0;
- int count = 0,count_1 = 0;
- COMPLEX poles[N];
- COMPLEX Res[N+1],Res_Save[N+1];
- if((N%2) == 0) dk = 0.5;
- else dk = 0;
- for(k = 0;k <= ((2*N)-1) ; k++)
- {
- if(Cotoff*cos((k+dk)*(pi/N)) < 0)
- {
- poles[count].Real_part = -Cotoff*cos((k+dk)*(pi/N));
- poles[count].Imag_Part= -Cotoff*sin((k+dk)*(pi/N));
- count++;
- if (count == N) break;
- }
- }
- printf("Pk = \n" );
- for(count = 0;count < N ;count++)
- {
- printf("(%lf) + (%lf i) \n" ,-poles[count].Real_part
- ,-poles[count].Imag_Part);
- }
- printf("--------------------------------------------------------\n" );
- Res[0].Real_part = poles[0].Real_part;
- Res[0].Imag_Part= poles[0].Imag_Part;
- Res[1].Real_part = 1;
- Res[1].Imag_Part= 0;
- for(count_1 = 0;count_1 < N-1;count_1++)
- {
- for(count = 0;count <= count_1 + 2;count++)
- {
- if(0 == count)
- {
- Complex_Multiple(Res[count], poles[count_1+1],
- &(Res_Save[count].Real_part),
- &(Res_Save[count].Imag_Part));
- //printf( "Res_Save : (%lf) + (%lf i) \n" ,Res_Save[0].Real_part,Res_Save[0].Imag_Part);
- }
- else if((count_1 + 2) == count)
- {
- Res_Save[count].Real_part += Res[count - 1].Real_part;
- Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;
- }
- else
- {
- Complex_Multiple(Res[count], poles[count_1+1],
- &(Res_Save[count].Real_part),
- &(Res_Save[count].Imag_Part));
- //printf( "Res : (%lf) + (%lf i) \n" ,Res[count - 1].Real_part,Res[count - 1].Imag_Part);
- //printf( "Res_Save : (%lf) + (%lf i) \n" ,Res_Save[count].Real_part,Res_Save[count].Imag_Part);
- Res_Save[count].Real_part += Res[count - 1].Real_part;
- Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;
- //printf( "Res_Save : (%lf) + (%lf i) \n" ,Res_Save[count].Real_part,Res_Save[count].Imag_Part);
- }
- //printf("There \n" );
- }
- for(count = 0;count <= N;count++)
- {
- Res[count].Real_part = Res_Save[count].Real_part;
- Res[count].Imag_Part= Res_Save[count].Imag_Part;
- *(a + N - count) = Res[count].Real_part;
- }
- //printf("There!! \n" );
- }
- *(b+N) = *(a+N);
- //------------------------display---------------------------------//
- printf("bs = [" );
- for(count = 0;count <= N ;count++)
- {
- printf("%lf ", *(b+count));
- }
- printf(" ] \n" );
- printf("as = [" );
- for(count = 0;count <= N ;count++)
- {
- printf("%lf ", *(a+count));
- }
- printf(" ] \n" );
- printf("--------------------------------------------------------\n" );
- return (int) 1;
- }
- int Bilinear(int N,
- double *as,double *bs,
- double *az,double *bz)
- {
- int Count = 0,Count_1 = 0,Count_2 = 0,Count_Z = 0;
- double Res[N+1];
- double Res_Save[N+1];
- for(Count_Z = 0;Count_Z <= N;Count_Z++)
- {
- *(az+Count_Z) = 0;
- *(bz+Count_Z) = 0;
- }
- for(Count = 0;Count<=N;Count++)
- {
- for(Count_Z = 0;Count_Z <= N;Count_Z++)
- {
- Res[Count_Z] = 0;
- Res_Save[Count_Z] = 0;
- }
- Res_Save [0] = 1;
- for(Count_1 = 0; Count_1 < N-Count;Count_1++)
- {
- for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
- {
- if(Count_2 == 0)
- {
- Res[Count_2] += Res_Save[Count_2];
- //printf( "Res[%d] %lf \n" , Count_2 ,Res[Count_2]);
- }
- else if((Count_2 == (Count_1+1))&&(Count_1 != 0))
- {
- Res[Count_2] += -Res_Save[Count_2 - 1];
- //printf( "Res[%d] %lf \n" , Count_2 ,Res[Count_2]);
- }
- else
- {
- Res[Count_2] += Res_Save[Count_2] - Res_Save[Count_2 - 1];
- //printf( "Res[%d] %lf \n" , Count_2 ,Res[Count_2]);
- }
- }
- //printf( "Res : ");
- for(Count_Z = 0;Count_Z<= N;Count_Z++)
- {
- Res_Save[Count_Z] = Res[Count_Z] ;
- Res[Count_Z] = 0;
- //printf( "[%d] %lf " ,Count_Z, Res_Save[Count_Z]);
- }
- //printf(" \n" );
- }
- for(Count_1 = (N-Count); Count_1 < N;Count_1++)
- {
- for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
- {
- if(Count_2 == 0)
- {
- Res[Count_2] += Res_Save[Count_2];
- //printf( "Res[%d] %lf \n" , Count_2 ,Res[Count_2]);
- }
- else if((Count_2 == (Count_1+1))&&(Count_1 != 0))
- {
- Res[Count_2] += Res_Save[Count_2 - 1];
- //printf( "Res[%d] %lf \n" , Count_2 ,Res[Count_2]);
- }
- else
- {
- Res[Count_2] += Res_Save[Count_2] + Res_Save[Count_2 - 1];
- //printf( "Res[%d] %lf \n" , Count_2 ,Res[Count_2]);
- }
- }
- // printf( "Res : ");
- for(Count_Z = 0;Count_Z<= N;Count_Z++)
- {
- Res_Save[Count_Z] = Res[Count_Z] ;
- Res[Count_Z] = 0;
- //printf( "[%d] %lf " ,Count_Z, Res_Save[Count_Z]);
- }
- //printf(" \n" );
- }
- //printf( "Res : ");
- for(Count_Z = 0;Count_Z<= N;Count_Z++)
- {
- *(az+Count_Z) += pow(2,N-Count) * (*(as+Count)) * Res_Save[Count_Z];
- *(bz+Count_Z) += (*(bs+Count)) * Res_Save[Count_Z];
- //printf( " %lf " ,*(bz+Count_Z));
- }
- //printf(" \n" );
- }
- for(Count_Z = N;Count_Z >= 0;Count_Z--)
- {
- *(bz+Count_Z) = (*(bz+Count_Z))/(*(az+0));
- *(az+Count_Z) = (*(az+Count_Z))/(*(az+0));
- }
- //------------------------display---------------------------------//
- printf("bz = [" );
- for(Count_Z= 0;Count_Z <= N ;Count_Z++)
- {
- printf("%lf ", *(bz+Count_Z));
- }
- printf(" ] \n" );
- printf("az = [" );
- for(Count_Z= 0;Count_Z <= N ;Count_Z++)
- {
- printf("%lf ", *(az+Count_Z));
- }
- printf(" ] \n" );
- printf("--------------------------------------------------------\n" );
- return (int) 1;
- }
- int main(void)
- {
- int count;
- struct DESIGN_SPECIFICATION IIR_Filter;
- IIR_Filter.Cotoff = (double)(pi/2); //[red]
- IIR_Filter.Stopband = (double)((pi*3)/4); //[red]
- IIR_Filter.Stopband_attenuation = 30; //[dB]
- int N;
- IIR_Filter.Cotoff = 2 * tan((IIR_Filter.Cotoff)/2); //[red/sec]
- IIR_Filter.Stopband = 2 * tan((IIR_Filter.Stopband)/2); //[red/sec]
- N = Buttord(IIR_Filter.Cotoff,
- IIR_Filter.Stopband,
- IIR_Filter.Stopband_attenuation);
- printf("N: %d \n" ,N);
- printf("--------------------------------------------------------\n" );
- double as[N+1] , bs[N+1];
- Butter(N,
- IIR_Filter.Cotoff,
- as,
- bs);
- double az[N+1] , bz[N+1];
- Bilinear(N,
- as,bs,
- az,bz);
- printf("Finish \n" );
- return (int)0;
- }
3.间接设计实现的IIR滤波器的性能
3.1设计指标
3.2程序执行结果
其频率响应如下所示。博客地址:http://blog.csdn.net/thnh169/
转载于:https://www.cnblogs.com/karl-wu/articles/4363409.html
[数字信号处理]IIR滤波器的间接设计(C代码)(转)相关推荐
- [数字信号处理]IIR滤波器基础
1.IIR滤波器构造 之前在介绍FIR滤波器的时候,我们提到过,IIR滤波器的单位冲击响应是无限的!用差分方程来表达一个滤波器,应该是下式这个样子的. 这个式子是N次差分方程的表达式.我们明显可以看出 ...
- 数字信号处理——FIR滤波器设计
一.线性相位. 并非所有的FIR滤波器都具有线性相位,只有当FIR滤波器的系数是对称的(包括偶对称和奇对称),才具有线性相位. 根据FIR滤波器的阶数(M=点数N-1)及其对称特性,可以分成以下四种类 ...
- 数字信号处理(七)FIR数字滤波器的设计
文章目录 FIR滤波器 线性相位FIR滤波器的条件及特点 线性相位FIR滤波器 线性相位条件 线性相位FIR滤波器幅度特性 类型1:h(n)=h(N-n-1),N=奇数 类型2:h(n)=h(N-n- ...
- Matlab | 数字信号处理:用窗函数法设计FIR数字滤波器
========================================== 博主github:https://github.com/MichaelBeechan 博主CSDN:https:/ ...
- 【数字信号处理】Python离散信号卷积的代码实现/时域直接法/列表法/信号与系统
Python离散信号卷积的代码实现(时域直接法) 1.卷积 卷积是一种积分变换的数学方法,在许多方面得到了广泛应用.卷积是两个变量在某范围内相乘后求和的结果.如果卷积的变量是离散的数组/数列,则卷积的 ...
- IIR滤波器设计(调用MATLAB IIR函数来实现)
转载请注明文章来源 – http://blog.csdn.net/v_hyx ,请勿用于任何商业用途 对于滤波器设计,以前虽然学过相关的理论(现代数字信号处理和DSP设计),但一直不求甚解,也没用过. ...
- 数字信号处理(六)IIR数字滤波器的设计
文章目录 数字滤波器 数字滤波器技术指标 数字低通滤波器的幅频响应曲线 IIR滤波器设计方法 IIR滤波器的函数模型设计法(间接法) 模拟低通滤波器的技术指标 模拟滤波器原型介绍 1.巴特沃斯模拟低通 ...
- FPGA数字信号处理(六)直接型IIR滤波器Verilog设计
该篇是FPGA数字信号处理的第六篇,2-5篇介绍了DSP系统中极其常用的FIR滤波器.本文将简单介绍另一种数字滤波器--IIR滤波器的原理,详细介绍使用Verilog HDL设计直接型IIR滤波器的方 ...
- matlab对图像信号进行频谱分析及滤波,数字信号处理课程设计---应用 Matlab对信号进行频谱分析及滤波...
数字信号处理课程设计---应用 Matlab对信号进行频谱分析及滤波 课课 程程 设设 计 (论文) 报计 (论文) 报 告告 书书 课程名称课程名称 数字信号处理 题题 目目 应用Matlab 对信 ...
最新文章
- windows下备份mysql 数据库
- 如何查看linux的资源,Linux系统资源查看(示例代码)
- 数据结构:列表(双向链表)的了解与示例
- 高中英语计算机辅助教学例子,计算机辅助教学在英语听力中的运用
- mongodb性能 mysql_MySQL和MongoDB的性能测试
- 丰田pcs可以关闭吗_论安全性能,广汽丰田TNGA车型如何?
- leetcode 326 [easy]--- Power of Three
- oracle临时表空间可以删除吗,Oracle临时表空间的增删改查
- 网络其他计算机无法访问,win7局域网别人无法访问我的电脑是为什么 win7其他电脑无法访问我的电脑如何修复...
- ctf工具整理-持续更新
- 百年铁树要开花,贾跃亭要还钱了?
- 利用Python去除图片水印,真的一点都不难!
- LeetCode笔记05:最长公共前缀
- com.monotype.android.font.ktoppo,Zawgyi Myanmar Fonts Free
- mybais-plus出现Invalid bound statement (not found)的解决方案
- pc端ui图片尺寸_pc端常用电脑屏幕 ((响应式PC端媒体查询)电脑屏幕分辨率尺寸大全)...
- 【CSP-J/S】复赛注意事项
- 反编译工具jad使用方法
- L1-044 稳赢分数 15(c++)
- UI设计网站 | 常用的UI设计网站大集合