1.IIR滤波器的直接设计原理
不利用模拟滤波器，直接进行数字滤波器的设计的方法，称为直接设计。回忆之前所说的IIR滤波器的直接设计，我们首先设计了巴特沃斯模拟滤波器，然后进行双线性变换，得到数字滤波器。我所使用的是巴特沃斯低通滤波器作为原型滤波器，其振幅特性如下所示。
首先，我们先把数字滤波器的指标，根据下式转为模拟滤波器的指标。
然后根据模拟滤波器的设计指标，计算次数N，然后计算极点，最后选择出稳定的极点，计算模拟滤波器的传递函数。然后，我们通过拉普拉斯反变换，再使用差分替代微分，得到了S平面-Z平面的对应关系，将其转为数字滤波器。
现在，我们将这个过程，稍微调整一下。
同样的，这个振幅特性具有巴特沃斯低通模拟滤波器的一些特性。函数为单调递减的，通带与阻带都没有纹波。这个新的滤波器，称为巴特沃斯低通数字滤波器。
与之前一样的，我们的首先先确定所需要的滤波器的次数N，根据次数计算出极点，选择稳定的极点来计算出传递函数，然后就可以得到滤波器的系数了。
首先是次数的计算，次数的计算，还是要根据阻带衰减指标，根据下式计算。
将巴特沃斯低通数字滤波器的振幅特性带入，我们就可以计算出所需要的滤波器的次数N。
下一步，我们是需要计算极点。这里与之前的间接设计不同，我们需要做一些变形。首先，将其振幅特性做平方，变为
然后，把
部分稍微做下变形，
部分稍微做下变形，
然后，将
替换为Z，带入，则得到了如下式子
替换为Z，带入，则得到了如下式子
通过这个式子，就可以很方便的计算极点与零点。很容易的能看出，这个滤波器的零点是-1，并且为N重极点（这里是振幅特性的平方所以不是2N）。此时，分母多项式为
由于一步解开很麻烦，我们先将这个式子，关于
解开，我们可以得到如下解，
解开，我们可以得到如下解，
再解z，然后再我们就可以得到极点了，如下所示。
这里所求得的极点为2N个，为了所设计的滤波器是稳定的，我们需要选择出稳定的极点。在Z平面内，其摸小于1的极点，就是稳定的极点。或者来说，只要滤波器的所有极点均在Z平面的单位圆里，那么，这个滤波器就是稳定的。（题外话，由于FIR滤波器的极点全是0，故FIR滤波器无需考虑稳定性。因为FIR滤波器必定稳定。）
最后，将极点带入传递函数的公式，就可以得到滤波器的系数了。
这里的K可以根据下式去求，
到此，我们就设计出了一个IIR数字滤波器。
首先，还是次数的计算。代码如下
- N = Ceil(0.5*( log10 ( pow (10, IIR_Filter.Stopband_attenuation/10) - 1) /
- log10 (IIR_Filter.Stopband/IIR_Filter.Cotoff)));
然后，是为了计算极点
，我们先计算
。这里，将指数函数换为了三角函数。
- poles_1.Real_part = (0.5)*Cotoff*cos((k+dk)*(pi/N));
- poles_1.Imag_Part= (0.5)*Cotoff*sin((k+dk)*(pi/N));
计算出
之后，我们计算极点
，然后选择稳定的极点。
之后，我们计算极点
，然后选择稳定的极点。
- for(k = 0;k <= ((2*N)-1) ; k++)
- {
- poles_1.Real_part = (0.5)*Cotoff*cos((k+dk)*(pi/N));
- poles_1.Imag_Part= (0.5)*Cotoff*sin((k+dk)*(pi/N));
- poles_2.Real_part = 1 - poles_1.Real_part ;
- poles_2.Imag_Part= -poles_1.Imag_Part;
- poles_1.Real_part = poles_1.Real_part + 1;
- poles_1.Real_part = poles_1.Real_part;
- Complex_Division(poles_1,poles_2,
- &poles[count].Real_part,
- &poles[count].Imag_Part);
- if(Complex_Abs(poles[count])<1)
- {
- poles[count].Real_part = -poles[count].Real_part;
- poles[count].Imag_Part= -poles[count].Imag_Part;
- count++;
- if (count == N) break;
- }
- }
- int Complex_Division(COMPLEX a,COMPLEX b,
- double *Res_Real,double *Res_Imag)
- {
- *(Res_Real) = ((a.Real_part)*(b.Real_part) + (a.Imag_Part)*(b.Imag_Part))/
- ((b.Real_part)*(b.Real_part) + (b.Imag_Part)*(b.Imag_Part));
- *(Res_Imag)= ((a.Real_part)*(b.Imag_Part) - (a.Imag_Part)*(b.Real_part))/
- ((b.Real_part)*(b.Real_part) + (b.Imag_Part)*(b.Imag_Part));
- return (int)1;
- }
- double Complex_Abs(COMPLEX a)
- {
- return (double)(sqrt((a.Real_part)*(a.Real_part) + (a.Imag_Part)*(a.Imag_Part)));
- }
- double K_z = 0.0;
- for(count = 0;count <= N;count++) {K_z += *(az+count);}
- K_z = (K_z/pow ((double)2,N));
2.IIR滤波器的直接设计代码（C语言）
- #include <stdio.h>
- #include <math.h>
- #include <malloc.h>
- #include <string.h>
- #define pi ((double)3.1415926)
- typedef struct
- {
- double Real_part;
- double Imag_Part;
- } COMPLEX;
- struct DESIGN_SPECIFICATION
- {
- double Cotoff;
- double Stopband;
- double Stopband_attenuation;
- };
- int Ceil(double input)
- {
- if(input != (double)((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 Complex_Division(COMPLEX a,COMPLEX b,
- double *Res_Real,double *Res_Imag)
- {
- *(Res_Real) = ((a.Real_part)*(b.Real_part) + (a.Imag_Part)*(b.Imag_Part))/
- ((b.Real_part)*(b.Real_part) + (b.Imag_Part)*(b.Imag_Part));
- *(Res_Imag)= ((a.Real_part)*(b.Imag_Part) - (a.Imag_Part)*(b.Real_part))/
- ((b.Real_part)*(b.Real_part) + (b.Imag_Part)*(b.Imag_Part));
- return (int)1;
- }
- double Complex_Abs(COMPLEX a)
- {
- return (double)(sqrt((a.Real_part)*(a.Real_part) + (a.Imag_Part)*(a.Imag_Part)));
- }
- double IIRFilter (double *a, int Lenth_a,
- double *b, int Lenth_b,
- double Input_Data,
- double *Memory_Buffer)
- {
- int Count;
- double Output_Data = 0;
- int Memory_Lenth = 0;
- if(Lenth_a >= Lenth_b) Memory_Lenth = Lenth_a;
- else Memory_Lenth = Lenth_b;
- Output_Data += (*a) * Input_Data; //a(0)*x(n)
- for(Count = 1; Count < Lenth_a ;Count++)
- {
- Output_Data -= (*(a + Count)) *
- (*(Memory_Buffer + (Memory_Lenth - 1) - Count));
- }
- //------------------------save data--------------------------//
- *(Memory_Buffer + Memory_Lenth - 1) = Output_Data;
- Output_Data = 0;
- //----------------------------------------------------------//
- for(Count = 0; Count < Lenth_b ;Count++)
- {
- Output_Data += (*(b + Count)) *
- (*(Memory_Buffer + (Memory_Lenth - 1) - Count));
- }
- //------------------------move data--------------------------//
- for(Count = 0 ; Count < Memory_Lenth -1 ; Count++)
- {
- *(Memory_Buffer + Count) = *(Memory_Buffer + Count + 1);
- }
- *(Memory_Buffer + Memory_Lenth - 1) = 0;
- //-----------------------------------------------------------//
- return (double)Output_Data;
- }
- int Direct( double Cotoff,
- double Stopband,
- double Stopband_attenuation,
- int N,
- double *az,double *bz)
- {
- printf("Wc = %lf [rad/sec] \n" ,Cotoff);
- printf("Ws = %lf [rad/sec] \n" ,Stopband);
- printf("As = %lf [dB] \n" ,Stopband_attenuation);
- printf("--------------------------------------------------------\n" );
- printf("N: %d \n" ,N);
- printf("--------------------------------------------------------\n" );
- COMPLEX poles[N],poles_1,poles_2;
- double dk = 0;
- int k = 0;
- int count = 0,count_1 = 0;;
- if((N%2) == 0) dk = 0.5;
- else dk = 0;
- for(k = 0;k <= ((2*N)-1) ; k++)
- {
- poles_1.Real_part = (0.5)*Cotoff*cos((k+dk)*(pi/N));
- poles_1.Imag_Part= (0.5)*Cotoff*sin((k+dk)*(pi/N));
- poles_2.Real_part = 1 - poles_1.Real_part ;
- poles_2.Imag_Part= -poles_1.Imag_Part;
- poles_1.Real_part = poles_1.Real_part + 1;
- poles_1.Real_part = poles_1.Real_part;
- Complex_Division(poles_1,poles_2,
- &poles[count].Real_part,
- &poles[count].Imag_Part);
- if(Complex_Abs(poles[count])<1)
- {
- poles[count].Real_part = -poles[count].Real_part;
- poles[count].Imag_Part= -poles[count].Imag_Part;
- 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" );
- COMPLEX Res[N+1],Res_Save[N+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));
- Res_Save[count].Real_part += Res[count - 1].Real_part;
- Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;
- }
- }
- for(count = 0;count <= N;count++)
- {
- Res[count].Real_part = Res_Save[count].Real_part;
- Res[count].Imag_Part= Res_Save[count].Imag_Part;
- *(az + N - count) = Res[count].Real_part;
- }
- }
- double K_z = 0.0;
- for(count = 0;count <= N;count++) {K_z += *(az+count);}
- K_z = (K_z/pow ((double)2,N));
- printf("K = %lf \n" , K_z);
- for(count = 0;count <= N;count++)
- {
- Res[count].Real_part = 0;
- Res[count].Imag_Part= 0;
- Res_Save[count].Real_part = 0;
- Res_Save[count].Imag_Part= 0;
- }
- COMPLEX zero;
- zero.Real_part = 1;
- zero.Imag_Part = 0;
- Res[0].Real_part = 1;
- Res[0].Imag_Part= 0;
- 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], zero,
- &(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],zero,
- &(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;
- }
- }
- for(count = 0;count <= N;count++)
- {
- Res[count].Real_part = Res_Save[count].Real_part;
- Res[count].Imag_Part= Res_Save[count].Imag_Part;
- *(bz + N - count) = Res[count].Real_part;
- }
- }
- for(count = 0;count <= N;count++)
- {
- *(bz + N - count) = *(bz + N - count) * K_z;
- }
- //------------------------display---------------------------------//
- printf("bz = [" );
- for(count= 0;count <= N ;count++)
- {
- printf("%lf ", *(bz+count));
- }
- printf(" ] \n" );
- printf("az = [" );
- for(count= 0;count <= N ;count++)
- {
- printf("%lf ", *(az+count));
- }
- printf(" ] \n" );
- printf("--------------------------------------------------------\n" );
- return (int)1;
- }
- int main(void)
- {
- int count;
- struct DESIGN_SPECIFICATION IIR_Filter;
- IIR_Filter.Cotoff = (double)(pi/4); //[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 = Ceil(0.5*( log10 ( pow (10, IIR_Filter.Stopband_attenuation/10) - 1) /
- log10 (IIR_Filter.Stopband/IIR_Filter.Cotoff)));
- double az[N+1] , bz[N+1];
- Direct(IIR_Filter.Cotoff,
- IIR_Filter.Stopband,
- IIR_Filter.Stopband_attenuation,
- N,
- az,bz);
- double *Memory_Buffer;
- Memory_Buffer = (double *) malloc(sizeof(double)*(N+1));
- memset(Memory_Buffer,
- 0,
- sizeof(double)*(N+1));
- FILE* Input_Data;
- FILE* Output_Data;
- double Input = 0 ;
- double Output = 0;
- Input_Data = fopen("input.dat","r");
- Output_Data = fopen("output.txt","w");
- while(1)
- {
- if(fscanf(Input_Data, "%lf", &Input) == EOF) break;
- Output = IIRFilter( az, (N+1),
- bz, (N+1),
- Input,
- Memory_Buffer );
- fprintf(Output_Data,"%lf,",Output);
- }
- printf("Finish \n" );
- return (int)0;
- }
3.IIR滤波器的直接设计(例)
3.1 设计指标
3.2 程序执行结果
根据之前所说的原理，使用前面的程序，我们可以得到满足设计指标的IIR滤波器的参数。程序在gcc下编译通过，执行结果如下
其理想的输出为0.5kHz正弦信号，实际输出为下。
其中，红色的----是Matlab计算的输出，粉色的o是用Ｃ语言计算的输出，蓝色的线是理想输出（也就是混合前的0.5kHz信号）。
4601
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
