MATLAB中有这个函数
说明
语法
[z,p,k]=butter(n,Wn)
[z,p,k] = butter(n,Wn,'ftype')
[b,a]=butter(n,Wn)
[b,a]=butter(n,Wn,'ftype')
[A,B,C,D]=butter(n,Wn)
[A,B,C,D] = butter(n,Wn,'ftype')
[z,p,k]=butter(n,Wn,'s')
[z,p,k] = butter(n,Wn,'ftype','s')
[b,a]=butter(n,Wn,'s')
[b,a]=butter(n,Wn,'ftype','s')
[A,B,C,D]=butter(n,Wn,'s')
[A,B,C,D] = butter(n,Wn,'ftype','s')
描述
butter 用来设计低通、带通、高通、和带阻数字和模拟的巴特沃斯滤波器。巴特沃斯滤波器的特征是通带内幅度响应最大平坦，且整体上是单调的。
巴特沃斯滤波器牺牲了在通带和阻带内的单调衰减陡度。除非需要巴特沃斯滤波器的平滑性，椭圆或切比雪夫滤波器可以用更小的滤波器阶数获得更陡峭的衰减特性。
数字域
[z,p,k] = butter(n,Wn) 设计一个阶数为n，归一化截止频率为Wn的低通数字巴特沃斯滤波器。此函数用n列的向量z和p返回零点和极点，以及用标量k返回增益。
[z,p,k] = butter(n,Wn,'ftype') 设计一个高通、低通或带阻滤波器，字符串'ftype'取值是：
'high' 用于设计归一化截止频率为Wn的高通数字滤波器
'low' 用于设计归一化截止频率为Wn的低通数字滤波器
'stop' 用于设计阶数为2*n的带阻数字滤波器，Wn应该是有两个元素的向量Wn=[w1 w2]。阻带是w1 <ω <w2.
截止频率 是幅度响应为处的的频率。对巴特沃斯滤波器，归一化截止频率Wn必须是介于0和1之间的数，这里的1对应于尼奎斯特频率，即每秒π弧度(π rad/s)。
如果Wn是含有两个元素的向量，Wn=[w1 w2]，butter 返回阶数为 2*n的数字带通滤波器，通频带为w1 <ω <w2.
通过返回不同数量的输出参数，butter 直接地得到其它的滤波器实现。要获得传输函数形式，使用如下所示的两个输出参数。
注意 参考下面的限制 了解关于影响形成传输函数的数值问题。
[b,a] = butter(n,Wn) 设计一个阶为n，归一化截止频率为Wn的数字低通巴特沃斯滤波器。它返回滤波器系数在长度为n+1的行向量b和a中，这两个向量包含z的降幂系数。
[b,a] = butter(n,Wn,'ftype') 设计一个高通、低通或带阻滤波器，字符串'ftype' 是上面描述的'high'、'low'、或 'stop'。
要获得状态空间形式，使用下面所示的4个输出参数：
[A,B,C,D] = butter(n,Wn) 或
[A,B,C,D] = butter(n,Wn,'ftype') 其中 A、 B,、C,、和D 是
并且u是输入, x是状态向量, y 是输出。
模拟域
[z,p,k] = butter(n,Wn,'s') 设计一个阶n，截止角频率为Wn rad/s的模拟低通巴特沃斯滤波器。它返回零点和极点在长n或2*n的列向量z和p中，标量k返回增益。butter的截止角频率Wn必须大于0 rad/s。
如果Wn是有两个元素w1<w2的向量， butter(n,Wn,'s') 返回阶 2*n 带通模拟滤波器，其通带是w1 <ω <w2。
[z,p,k] = butter(n,Wn,'ftype','s') 通过使用上面描述的ftype 值可以设计一个高通、低通或带阻滤波器。
只要返回不同数量的输出参数，butter 可以直接地获得其它的模拟滤波器实现。要获得传输函数形式，使用如下所示的两个输出参数：
[b,a] = butter(n,Wn,'s') 设计一个阶n、截止角频率为Wn rad/s的模拟低通巴特沃斯滤波器。它返回滤波器的系数在长n+1的行向量b和a中，这两个向量包含下面这个传输函数中s的降幂系数：
[b,a] = butter(n,Wn,'ftype','s') 通过设置上面描述的ftype 值，可以设计一个高通、低通或带阻滤波器。
要获得状态空间形式，使用下面的四个参数：
[A,B,C,D] = butter(n,Wn,'s') 或
[A,B,C,D] = butter(n,Wn,'ftype','s') 其中A、 B、 C、和D 是
并且u 作为输入, x 是状态向量, y 是输出。
举例
高通滤波器
对于1000Hz的采样，设计一个9阶高通巴特沃斯滤波器，截止频率300Hz，相应的归一化值为0.6：
[z,p,k] = butter(9,300/500,'high')
[sos,g] = zp2sos(z,p,k) % 转换为二次分式表示形式
Hd = dfilt.df2tsos(sos,g) % 创建dfilt对象
h = fvtool(Hd) % 绘制幅度响应
set(h,'Analysis','freq') % 显示频率响应
#define pi 3.14159265358979384626double h(double n)
{
double x
x=1.0/(pi*(n-64))
x*=sin(pi*(n-64))-sin(3*pi/4*(n-64))
x*=0.54-0.46*cos(2*pi*n/128)
x*=r128(n)
return x
}
double r128(double n)
{
double x
/*函数R128的公式你没给出，你自己完成吧*/
return x
}
/// 巴特沃斯滤波器, 带增益Mat BoostButterworthFilter( const int &row, const int &col, const float &cutoff, const int &n, const float &boost )const
{
assert( row>1 && col>1 )
assert( cutoff>0 && cutoff<0.5 )
assert( n>=1 )
if ( boost>=1 )
{
return (1-1/boost)*( BHPF( row, col, cutoff, n ) ) + 1/boost
}
else
{
return (1-boost)*BLPF( row, col, cutoff, n) + boost
}
}
/// 巴特沃斯高通滤波器
Mat CIllumNorm::BHPF( const int &row, const int &col, const float &cutoff, const int &n )
{
return 1-BLPF( row, col, cutoff, n )
}
/// 巴特沃斯低通滤波器
/*
cutoff is the cutoff frequency of the filter 0 - 0.5
n is the order of the filter, the higher n is the sharper
1
f = --------------------
2n
1.0 + (w/cutoff)
*/
Mat BLPF( const int &row, const int &col, const float &cutoff, const int &n )
{
assert( row>1 && col>1 )
assert( cutoff>0 && cutoff<0.5 )
assert( n>=1 )
Mat X = Mat::zeros( row, col, CV_32F )
Mat Y = Mat::zeros( row, col, CV_32F )
if ( col%2 )
{
int t = -(col-1)/2
for ( int j=0 j<col j++ )
{
for ( int i=0 i<row i++ )
{
X.at<float>( i, j ) = ((float)t)/(col-1)
}
t++
}
}
else
{
int t = -col/2
for ( int j=0 j<col j++ )
{
for ( int i=0 i<row i++ )
{
X.at<float>( i, j ) = ((float)t)/col
}
t++
}
}
if ( row%2 )
{
int t = -(row-1)/2
for ( int i=0 i<row i++ )
{
for ( int j=0 j<col j++ )
{
Y.at<float>( i, j ) = ((float)t)/(row-1)
}
t++
}
}
else
{
int t = -row/2
for ( int i=0 i<row i++ )
{
for ( int j=0 j<col j++ )
{
Y.at<float>( i, j ) = ((float)t)/row
}
t++
}
}
Mat H = Mat::zeros( row, col, CV_32F )
/// 计算频域的巴特沃斯分类器
for ( int i=0 i<row i++ )
{
for ( int j=0 j<col j++ )
{
float x = X.at<float>( i, j )
float y = Y.at<float>( i, j )
float radius = sqrtf( x*x + y*y )
/// 1.0 ./ (1.0 + (radius ./ cutoff).^(2*n))
H.at<float>( i, j ) = 1.0 / ( 1.0+std::pow( radius/cutoff, 2*n ) )
}
}
/// shift 将中心点移到左上角
H = CenterShift(H)
return H
}
Mat CenterShift( const Mat &_src )
{
Mat src = Mat_<float>(_src)
Mat dst = Mat::zeros( src.rows, src.cols, src.type() )
int N = src.rows
int D = src.cols
int *rowIndex = new int[N]
int *colIndex =new int[D]
memset( rowIndex, 0, sizeof(rowIndex[0])*N )
memset( colIndex, 0, sizeof(colIndex[0])*D )
/// 行 顺序
int begin = N/2
for ( int i=0 i<N i++ )
{
rowIndex[i] = begin
begin++
if ( begin>=N )
{
begin = 0
}
}
/// 列 顺序
begin = D/2
for ( int i=0 i<D i++ )
{
colIndex[i] = begin
begin++
if ( begin>=D )
{
begin = 0
}
}
/// 重新排序
for ( int i=0 i<N i++ )
{
for ( int j=0 j<D j++ )
{
dst.at<float>( i, j ) = src.at<float>( rowIndex[i], colIndex[j] )
}
}
/// 释放
delete []rowIndex
delete []colIndex
rowIndex = NULL
colIndex = NULL
return dst
}
