关于matlab如何辅助设计滤波器最近大致搞明白一些,这里主要是一些比较传统的数字滤波器设计方法(还有一些直接转换到CCS,VHDL或者C Header的这里没有提到),通过matlab算出相应的传递函数,然后C语言进行仿真测试,最后再改写到所需要用的设备上。
首先谈谈matlab的设计方式,主要两种
1.辅助设计 Filter Designe & analysis Tool,通过辅助工具设计,比较省事
进去观察下有哪些选项可以设置
大部分还是基本看了就知道需要怎么设,然后就是准备导出传递函数
FILE中的EXPORT,然后设置10进制输出
最后就能导出相应的传函
%
% Generated by MATLAB(R) 7.8 and the Signal Processing Toolbox 6.11.
%
% Generated on: 06-Mar-2012 23:13:19
%
% Coefficient Format: Decimal
% Discrete-Time IIR Filter (real)
% -------------------------------
% Filter Structure : Direct-Form II, Second-Order Sections
% Number of Sections : 2
% Stable : Yes
% Linear Phase : No
SOS matrix:
1 -0.3274034970481377 1 1 -1.487740666582311 0.62366222645998592
1 -1.4109006645592876 1 1 -1.5591910986997555 0.90887414412475032
Scale Values:
0.072426444477067295
0.59358927160118591
得出一个SOS矩阵
这个可以通过下图看下怎么转换
对于以零点-极点-增益为参数的传递函数H(Z),可用函数zp2sos将其转换为LX6sos数组。函数sos2ss将二阶滤波器节映射为直接II型状态空间形式,而函数sos2tf将一组 2阶滤波器节转换为传递函数。对于滤波器节为1阶的情况,系数b2i和a2i将被没为零。函数sos2zp将2阶滤波器节转换为零点-极点-增益的形式。
2.通过matlab函数进行实现由于设计的方法挺多了,就不列举了,直接通过函数调用就行了
这里主要是为了方便实现转换为比较容易实现的是级联的设计方式
相关的函数
function [b0,B,A] = dir2cas(b,a);
% 直接型到级联型的型式转换(复数对型)
% ---------------------------------------------------------
% [b0,B,A] = dir2cas(b,a)
% b = 直接型的分子多项式系数
% a = 直接型的分母多项式系数
% b0 = 增益系数
% B = 包含各bk的K乘3维实系数矩阵
% A = 包含各ak的K乘3维实系数矩阵
% compute gain coefficient b0
b0 = b(1); b = b/b0;
a0 = a(1); a = a/a0;
b0 = b0/a0;
%
M = length(b); N = length(a);
if N > M
b = [b zeros(1,N-M)];
elseif M > N
a = [a zeros(1,M-N)]; N = M;
else
NM = 0;
end
%
K = floor(N/2); B = zeros(K,3); A = zeros(K,3);
if K*2 == N;
b = [b 0];
a = [a 0];
end
%
broots = cplxpair(roots(b));
aroots = cplxpair(roots(a));
for i=1:2:2*K
Brow = broots(i:1:i+1,:);
Brow = real(poly(Brow));
B(fix((i+1)/2),:) = Brow;
Arow = aroots(i:1:i+1,:);
Arow = real(poly(Arow));
A(fix((i+1)/2),:) = Arow;
end
获得了相关的传函之后,就可以进行测试了,测试也是两部分
首先通过matlab测试一些相关的幅频特性
%对于有理分式H(jw),MATLAB提供freqs函数处理方法。调用格式为H=freqs(b,a,w)
%其中,b为分子多项式的系数,a为分母多项式系数,w为需计算的频率特性
clear all
clc
w=linspace(0,5,200);
b=[1];
a=[1 2 2 1];
H=freqs(b,a,w);
plot(w,abs(H));
测试下
其次是用C语言大致模拟下,大致也体现了怎么使用传递函数,如下
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <conio.h>
#include <dos.h>
#include <time.h>
#include <io.h>
#include <process.h>
void main()
{
int i;
int x0,x1,x2,y0,y1,y2,q0,q1,q2;
FILE *fp;
x0=0;x1=0;x2=0;
y0=0;y1=0;y2=0;
q0=0;q1=0;q2=0;
system("cls");
fp=fopen("C:\\data0.txt","wt");
printf("Press any key to continue when ready\n");
printf("or Press ESE to Cancel\n");
getch();
for(i=0;i<500;i++)
{
x2=1600;
if(i==200) x2=3200;
if(i==350) x2=800;
y2=(113*x2-107*x1+113*x0+1618*y1-662*y0)/1024;
q2=(252*y2-436*y1+252*y0+1818*q1-913*q0)/1024;
y0=y1;y1=y2;x0=x1;x1=x2;q0=q1;q1=q2;
printf("%d %d\n",x2/16,q2/16);
fprintf(fp,"%d %d\n",x2/16,q2/16);
}
fclose(fp);
system("pause");
}
这里113,107,113,1618,662为传函的参数调理之后,如下
这样大致就设计好了,这个C语言的版本只是一个直流信号的测试,不过大致可以明白,将传函的微分方程改写成差分方程,并进行计算,通过这个思想,可以很方便的将其移植到其他的单片机或者FPGA上,还是比较好实现的。