matlab buttord函数 C实现

最新推荐文章于 2024-05-17 03:41:42 发布
最新推荐文章于 2024-05-17 03:41:42 发布
阅读量778 收藏 7
点赞数 10
分类专栏： 数字滤波器 文章标签： matlab c语言 信号处理
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/CSharp_Beginner/article/details/135739389
版权
本文介绍了用于滤波器设计的buttord函数，该函数根据给定的参数如通带波纹、阻带衰减和类型，计算所需滤波器的最小阶数，并确定相应的频率数组。函数还包括lclfminbnd和bscost辅助函数，确保了与Matlab等工具的兼容性。

摘要生成于 C知道 ，由 DeepSeek-R1 满血版支持， 前往体验 >

buttord

/************************************************************************
 * w是频率范围数组，由ws和wp组成
 * nw为w数组长度,低通或高通长度为2，带通或带阻长度为4
 * rp为通带波纹
 * rs为阻带衰减
 * order为所求最小阶数
 * wn为所求最小阶数对应的频率数组
*/
void buttord(double* w,int nw,double rp,double rs,int* order,double* wn)
{
	int type=0;
	double* wp=(double*)malloc(sizeof(double)*nw/2);
	double* ws=(double*)malloc(sizeof(double)*nw/2);
	double* WA=(double*)malloc(sizeof(double)*nw/2);
	double* WN=(double*)malloc(sizeof(double)*nw/2);
	double minWA, W0;

	for(int i=0;i<nw/2;i++)
	{
		wp[i]=w[i];
		ws[i]=w[nw/2+i];
	}

	if(nw==2)
	{
		type=0;
	}
	else if(nw==4)
	{
		type=2;
	}
	if(wp[0]<ws[0])
	{
		type+=1;
	}
	else
	{
		type+=2;
	}


	for(int i=0;i<nw/2;i++)
	{
		wp[i]=tan(PI*wp[i]/2.0);
		ws[i]=tan(PI*ws[i]/2.0);
	}

	if (type == 1)
	{
		for (int i = 0; i < nw / 2; i++)
		{
			WA[i] = ws[i] / wp[i];
		}
	}
	else if(type==2)
	{
		for (int i = 0; i < nw / 2; i++)
		{
			WA[i] = wp[i] / ws[i];
		}
	}
	else if (type == 3)
	{
		wp[0] = lclfminbnd(wp[0], ws[0] - 1e-12,0, wp, ws, rs, rp, 0);
		wp[1] = lclfminbnd(ws[1] + 1e-12, wp[1],1, wp, ws, rs, rp, 0);		
	}
	else if(type==4)
	{
		for (int i = 0; i < nw / 2; i++)
		{
			WA[i] = (ws[i] * ws[i] - wp[0] * wp[1]) / (ws[i] * (wp[0] - wp[1]));
		}
	}

	minWA=INFINITY;
	for (int i = 0; i < nw / 2; i++)
	{
		double absValue = fabs(WA[i]);
		minWA = fmin(minWA, absValue);
	}

	*order =  ceil(log10((pow(10, 0.1 * fabs(rs)) - 1.0)/(pow(10, 0.1 * fabs(rp)) - 1.0)) / (2.0 * log10(minWA)));

	W0 = minWA / (pow((pow(10, .1 * fabs(rs)) - 1), (1 / (2 * fabs(*order)))));

	if (type == 1) {  // low
        WN[0] = W0 * wp[0];
    } else if (type == 2) {  // high
        WN[0] = wp[0] / W0;
    } else if (type == 3) {  // stop
        WN[0] = fabs(((wp[1] - wp[0]) + sqrt(pow((wp[1] - wp[0]), 2) + 4 * pow(W0, 2) * wp[0] * wp[1])) / (2 * W0));
        WN[1] = fabs(((wp[1] - wp[0]) - sqrt(pow((wp[1] - wp[0]), 2) + 4 * pow(W0, 2) * wp[0] * wp[1])) / (2 * W0));
        // Sorting WN in ascending order
        if (WN[0] > WN[1]) {
            double temp = WN[0];
            WN[0] = WN[1];
            WN[1] = temp;
        }
    } else if (type == 4) {  // pass
        double W0Array[2] = {-W0, W0};
        WN[0] = -W0Array[0] * (wp[1] - wp[0]) / 2 + sqrt(pow(W0Array[0], 2) / 4 * pow((wp[1] - wp[0]), 2) + wp[0] * wp[1]);
        WN[1] = -W0Array[1] * (wp[1] - wp[0]) / 2 + sqrt(pow(W0Array[1], 2) / 4 * pow((wp[1] - wp[0]), 2) + wp[0] * wp[1]);
        // Sorting WN in ascending order
        if (WN[0] > WN[1]) {
            double temp = WN[0];
            WN[0] = WN[1];
            WN[1] = temp;
        }
    }
	for (int i = 0; i < nw / 2; i++) {
        wn[i]=(2/PI)*atan(WN[i]);
    }

	free(wp);
	free(ws);
	free(WA);
	free(WN);
}

其中lclfminbnd代码如下：


double lclfminbnd(double ax, double bx,int ind,double* wp,double* ws,double rs,double rp,int type) 
{
	// Initialization
    double seps = sqrt(EPS);
    double c = 0.5 * (3.0 - sqrt(5.0));
    double a = ax, b = bx;
    double v = a + c * (b - a);
    double w = v, xf = v;
    double d = 0.0, e = 0.0;
    double x = xf;
    double fx = bscost(x,ind,wp,ws,2,rs,rp,type);
    int num = 1;
    double fv = fx, fw = fx;
    double xm = 0.5 * (a + b);
	double tol=1.0e-4;
    double tol1 = seps * fabs(xf) + tol / 3.0;
    double tol2 = 2.0 * tol1;

	// Main loop
	while (fabs(xf - xm) > (tol2 - 0.5 * (b - a)))
	{
		int gs = 1;
		if (fabs(e) > tol1)
		{
			gs = 0;
			// Parabolic interpolation
			double r = (xf - w) * (fx - fv);
			double q = (xf - v) * (fx - fw);
			double p = (xf - v) * q - (xf - w) * r;
			q = 2.0 * (q - r);
			if (q > 0.0)
			{
				p = -p;
			}
			q = fabs(q);
			double temp = e;
			e = d;
			if ((fabs(p) < fabs(0.5 * q * temp)) && (p > q * (a - xf)) && (p < q * (b - xf)))
			{
				d = p / q;
				x = xf + d;
				// Parabolic interpolation step
				if ((x - a) < tol2 || (b - x) < tol2)
				{
					int si = (xm - xf > 0) - (xm - xf < 0);
					d = tol1 * si;
				}
			}
			else
			{
				gs = 1;
			}
		}
		if (gs)
		{
			if (xf >= xm)
			{
				e = a - xf;
			}
			else
			{
				e = b - xf;
			}
			d = c * e;			
		}
		double si = d >= 0?1.0:-1.0;
		x = xf + si * fmax(fabs(d), tol1);
		double fu = bscost(x,ind,wp,ws,2,rs,rp,type);
		num++;
		// Update variables
		if (fu <= fx)
		{
			if (x >= xf)
			{
				a = xf;
			}
			else
			{
				b = xf;
			}
			v = w;
			fv = fw;
			w = xf;
			fw = fx;
			xf = x;
			fx = fu;
		}
		else
		{ // fu > fx
			if (x < xf)
			{
				a = x;
			}
			else
			{
				b = x;
			}
			if ((fu <= fw) || (w == xf))
			{
				v = w;
				fv = fw;
				w = x;
				fw = fu;
			}
			else if ((fu <= fv) || (v == xf) || (v == w))
			{
				v = x;
				fv = fu;
			}
		}
		xm = 0.5 * (a + b);
		tol1 = seps * fabs(xf) + tol / 3.0;
		tol2 = 2.0 * tol1;
	}

	return xf;
}

又涉及到bscost函数实现：

double bscost(double wp,int ind,double* WP,double* WS,int n,double rs,double rp,int type)
{
    double result;
    double* WA=(double*)malloc(sizeof(double)*n);
    double* tempwp=(double*)malloc(sizeof(double)*n);
    double minWA=INFINITY;
    for (int i = 0; i < n; i++)
    {
        tempwp[i]=WP[i];
    }
    

    tempwp[ind]=wp;
    for (int i = 0; i < n; i++)
    {
       WA[i]=(WS[i]*(tempwp[0]-tempwp[1]))/(WS[i]*WS[i]  - tempwp[0]*tempwp[1]);
    }
    
    for (int i = 0; i < n; i++)
	{
		double absValue = fabs(WA[i]);
		minWA = fmin(minWA, absValue);
	}

    switch (type)
    {
        case 0: // Butter
        {

            result = log10((pow(10, 0.1 * fabs(rs)) - 1.0) / pow(10, 0.1 * fabs(rp)) - 1.0) / (2 * log10(minWA));
            break;
        }
        case 1: // cheby
        {
            result = acosh(sqrt((pow(10, 0.1 * fabs(rs)) - 1.0) / (pow(10, 0.1 * fabs(rp)) - 1.0))) / acosh(minWA);
            break;
        }
        case 2: // ellip
        {
            double epsilon;
            double k1;
            double k;
            // Assuming ellipke function is already defined
            double *ks = (double *)malloc(sizeof(double) * 2);
            double *ks1 = (double *)malloc(sizeof(double) * 2);
            epsilon = sqrt(pow(10, 0.1 * rp) - 1.0);
            k1 = epsilon / sqrt(pow(10, 0.1 * rs) - 1.0);
            k = 1.0 / minWA;
            ks[0] = k * k;
            ks[1] = 1 - k * k;
            ks1[0] = k1 * k1;
            ks1[1] = 1 - k1 * k1;
            double *capk = ellipke(ks, 2);
            double *capk1 = ellipke(ks1, 2);
            result = (capk[0] * capk1[1]) / (capk[1] * capk1[0]);
            free(ks);
            free(ks1);
            free(capk);
            free(capk1);
            break;
        }
        default:
            break;
    }

    free(tempwp);
    free(WA);
    return result;
    
}

以上函数都经过和matlab对应函数实现，验证结果一致。

MATLAB IIR滤波器设计函数buttord与butter(3)
2401_85015040的博客
05-17 913
如果Wp和Ws都是向量，并且由Ws指定的间隔包含由Wp指定的间隔（Ws（1）
butter滤波器是iir吗_MATLAB IIR滤波器设计函数buttord与butter
weixin_39559523的博客
12-20 749
(4)Ws —阻带截止频率标量| 二元向量阻带截止频率，指定为标量或二元向量，值在0到1之间，其中1对应于归一化的奈奎斯特频率πrad / sample。数据类型：单| 双(5)Rp-通带纹波标量通带纹波，指定为标量，单位为dB。数据类型：单| 双(6)Rs-阻带衰减标量阻带衰减，指定为标量，单位为dB。数据类型：单| 双Buttord的阶数预测公式适用于模拟和数字情况下的模拟范围。对于数字情况，...
C++ 实现 matlab 的 buttord函数
hugang的博客
04-16 2322
/ 数字通带截止频率 Hz double digital_passband_cutoff_frequency;// 数字阻带截止频率 Hz double digital_stopband_cutoff_frequency;// 通带最大衰减 db int passband_max_attenuation;// 阻带最小衰减 db int stopband_min_attenuation;// 采样频率 Hz int fs;
MATLAB中buttord的使用
zzy0222的博客
10-28 6万+
MATLAB中buttord的使用 我看网上很多对于buttord函数的使用都是正确的，我也是doc了一些MATLAB里面的buttord函数才知道的。 前半段代码的解释 有一个很关键的问题是信号的采样频率是1kHz，但是在实际归一化的过程中使用的是500Hz。这个500Hz实际上是奎因斯特频率 其实大家看看时域采样定理就可以知道，Nyquist频率实际上是限带信号的原始最高频率，就是那个Wm...
latex不等于_你可能不知道的MATLAB操作#第四话
weixin_39755136的博客
11-23 551
前言回顾 #第三话#介绍了MATLAB内部的绘图命令。比如如何利用函数命令控制数据图的坐标轴、线型、颜色等；如何绘制子图与动图；如何在标题或者坐标轴标签中增加LaTeX公式。此外，小编说一哈下写本系列推送的目的不是为了让大家对MATLAB的用法有全面的了解或者替代书本，因为本人也做不到。出发点是为大家整理最常用的命令和最容易犯的错误，可以指导大家更高效地自学。更多的自学资料...
MATLAB IIR滤波器设计函数buttord与butter(2)
最新发布
2401_85015025的博客
05-17 1160
如果Wp和Ws都是向量，并且由Ws指定的间隔包含由Wp指定的间隔（Ws（1）
matlab中buttord用法_matlab butter函数
weixin_39658019的博客
12-20 5149
匿名用户1级2013-12-01 回答语法[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]...
matlab buttord函数用法
07-28
buttord函数是MATLAB中的一个滤波器设计函数，用于计算Butterworth滤波器的阶数和截止频率。它的语法如下： [n, Wn] = buttord(Wp, Ws, Rp, Rs) 其中，Wp和Ws分别是通带和阻带的归一化频率，Rp和Rs是通带和阻带的...
matlab中buttord用法_MATLAB中buttord的使用
weixin_39762856的博客
12-30 1785
MATLAB中buttord的使用我看网上很多对于buttord函数的使用都不太清楚，我也是doc了一些MATLAB里面的buttord函数才知道的。 前半段代码的解释有一个很关键的问题是信号的采样频率是1kHz，但是在实际归一化的过程中使用的是500Hz。这个500Hz实际上是奎因斯特频率 其实大家看看时域采样定理就可以知道，Nyquist频率实际上是限带信号的原始最高频率，就是那个Wm/2*...
下载matlab后出现函数冲突,Matlab程序运行后出错buttord函数的问题
weixin_36322704的博客
03-22 739
本帖最后由 shuFeeYoung 于 2013-10-6 21:12 编辑错如提示:Error using buttord (line 59)Invalid option for final argument.Error in bianshengqi (line 42)[N,Wn]=buttord(wp1,ws1,Rp,Rs,'x1'); %选择滤波器的最小阶数请指教~怎么就用错了....代...
C语言实现matlab的butter函数(免complex.h)
11-28
代码验证请看我的博客: http://blog.csdn.net/zone53/article/details/78280901
Matlab关于数字信号处理的常用函数集锦
12-24
Matlab关于数字信号处理的常用函数集锦
频域滤波matlab代码_连续时间滤波器Matlab仿真实例
weixin_39997037的博客
12-13 826
使用Matlab对给定的信号和预期需求，设计了几种连续时间滤波器，包括四种最常见的连续滤波器：巴特沃斯(Butterworth)、切比雪夫滤波器I型和II型(chebyshev)、考尔滤波器(Cauer)；验证了滤波器性能，并进行了分析和对比，给出了代码，并将滤波器设计过程常用操作简单封装成函数，简化编程。待处理信号不同音所在的频率这里使用gamme.m函数，生成一段依次从Do、 Re、Mi、.....
matlab sprintf_如何在matlab用带通滤波器进行滤波？
weixin_39521068的博客
11-14 748
以带通滤波器——巴特沃斯滤波器为例matlab给出的代码如下：[n,Wn] = buttord(Wp,Ws,Rp,Rs)是不是显得很不人性化，即便是看了help中的英文解释，仍旧是对输入输出参数一知半解~那么我就用简单的语言来描述一下~首先要先理解一个带通滤波器的幅频特性图中，分为通带、阻带、过渡带三大部分。通带截止频率和阻带截止频率对这三部分进行划分，二者之间的部分即为过渡带。一.如何设计带通滤...
MATLAB IIR滤波器设计函数buttord与butter
芒果爱火锅的博客
04-21 3万+
设计巴特沃斯IIR滤波器可使用butter函数和buttord函数。 第一：使用buttord函数先求得最小阶数和截止频率：[n,Wn]= buttord(Wp,Ws,Rp,Rs)。 理解：上式中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。简单来说，就是在Wp处，通带内波纹系数或者说是通带内达到最大衰减为Rp，如（3db）,而在Ws处，阻带达到最小衰减为R...
MATLAB_数字信号处理_模拟滤波器_设计巴特沃斯滤波器
qq_45908742的博客
12-07 1万+
巴特沃斯滤波器-百度百科巴特沃斯滤波器-维基百科简介：巴特沃斯滤波器是一种模拟滤波器，它在频率响应方面具有特殊的属性。它被设计为具有均匀的幅度响应，即在通带内，它对所有频率的增益是相等的，而在阻带内，它对所有频率的增益都是零。为了实现这种频率响应，巴特沃斯滤波器在设计时采用了一些特殊的方法。具体而言，它是由一组极点和零点组成的，这些极点和零点的位置是根据所需的通带和阻带响应来计算的。在巴特沃斯滤波器中，极点和零点都是共轭对，这意味着它们的实部相等，虚部互为相反数。
数字信号处理|Matlab设计巴特沃斯低通滤波器（冲激响应不变法和双线性变换法）
zzztutu的博客
12-29 3万+
将巴特沃斯低通滤波器从理论到仿真完整实现
使用Matlab设计数字滤波器，从原理到代码
热门推荐
乐观的lishan的博客
12-05 9万+
本文是在Matlab下设计数字滤波器的一些总结，减少了理论的说明，重在代码的实现方法，给出了普遍使用的数字滤波器代码，修改参数即可移植使用。
巴特沃斯滤波（matlab/C++）
weixin_50624597的博客
07-13 1万+
文章目录一、滤波器简介二、用matlab做巴特沃斯低通滤波器2.1基本数据2.2做出原信号的频谱函数2.3做出巴特沃斯低通滤波器2.4用滤波器过滤信号 并得出频谱图2.5对高频的信号的低通滤波三、MATLAB中filter的理解与使用3.1filter概念与基本语法3.2以最简单的 y = filter(b,a,X) 为例3.3可实现差分方程四、巴特沃斯模拟滤波器（低通，高通，带通，带阻）设计-MATLAB实现4.1基础知识介绍4.2函数介绍4.2.1buttord - 求解滤波器的阶数N和3dB截止频率w
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值