Dsp相关问题思考笔记(一)

本文探讨了使用泰勒展开式对三角函数sinx进行近似计算的方法,特别是在工程应用中,通过控制展开阶数m确保信噪比不低于40dB。分析表明,当x取2π时,至少需要10项;若减小x的值,例如取π/2,可以将m降至3。通过编程验证,实际信噪比满足要求,证实了理论计算的准确性。

问题

三角函数 sinxsinxsinx可以通过泰勒展开写成形式:sinx=∑n=0+∞(−1)nx2n+1(2n+1)!sinx=\sum_{n=0}^{+\infty}{(-1)^{n}\frac{x^{2n+1}}{(2n+1)!}}sinx=n=0+(1)n(2n+1)!x2n+1工程上可以根据此公式,进行三角函数的近似计算,从而输出一个正选波。但是展开的阶数不可能取到无限项,应在满足信噪比的情况下,取到尽可能小的有限项。即用前m项进行近似:sinx=∑n=0m−1(−1)nx2n+1(2n+1)!+errorsinx=\sum_{n=0}^{m-1}{(-1)^{n}\frac{x^{2n+1}}{(2n+1)!}}+errorsinx=n=0m1(1)n(2n+1)!x2n+1+error现在考虑,在满足信噪比不低于40dB的前提下,满足的m最小值。

思考过程

1基本

error实际上是前一项代数式的高阶无穷小,具体为形式:error=x2m+1(2m+1)!−x2m+3(2m+3)!+……error=\frac{x^{2m+1}}{(2m+1)!}-\frac{x^{2m+3}}{(2m+3)!}+……error=(2m+1)!x2m+1(2m+3)!x2m+3+由于高阶无穷小可以忽略掉,可以近似认为,误差由第m+1项决定,即error≈x2m+1(2m+1)!error\approx\frac{x^{2m+1}}{(2m+1)!}error(2m+1)!x2m+1则信号的信噪比有SNR=signalpowernosiepower=12[x2m+1(2m+1)!]2≥104SNR=\frac{signal\quad power}{nosie\quad power}=\frac{\frac{1}{2}}{[\frac{x^{2m+1}}{(2m+1)!}]^{2}}\ge10^{4}SNR=nosiepowersignalpower=[(2m+1)!x2m+1]221104考虑最恶劣情况,则xxx取区间内的最大值,则信噪比最坏。正常输出sinxsinxsinx时,为避免累计误差,都会利用其周期性,即xxx的最大值为2π2\pi2π,则公式变为12[2π2m+1(2m+1)!]2≥104\frac{\frac{1}{2}}{[\frac{2\pi^{2m+1}}{(2m+1)!}]^{2}}\ge10^{4}[(2m+1)!2π2m+1]221104此公式为一个非线性方程,但可以知道,随着阶数mmm的增大,信噪比是不断提升的,求解非线性方程的经典算法是二分法。编写求解该非线性方程的C代码。

/*
Function:solve a nonlinear equation with dichotomy method.
Notes:this solve is limited to a integer.
*/
#include "stdio.h"
#include "math.h"

#define PI	3.1415926535897932384626433832795

double fact(int n)
{

     //	return (n==0) || (n==1) ? 1 : n* fact(n-1);

	double result=1;
	if(n==0)
		return 1;
	else
	{
		while(n!=1)
		{
			result*=n;
			n--;
		}
		return result;
	}
}
void main(void)
{
	double	expression,
			search_begin,
			search_end,
			temp;
	int m;
	search_begin=0;
	search_end=20;
	while(1<(search_end-search_begin))
	{
		m=(int)((search_begin+search_end)/2);
		temp=pow(2*PI,2*m+1)/fact(2*m+1);
		expression=0.5/(temp*temp)-1e4;
		if(expression<0)
			search_begin=m;
		else
			search_end=m;
	}
	if(expression<0)
		m++;
	printf("The minimum m is %d\r\n",m);
}

输出结果:

The minimum m is 10

2 进阶

sinxsinxsinxxxx取到2π2\pi2π时,满足信噪比最少需要取到10项。如果想要在保证信噪比前提下进一步减少项数,则要让xxx取值进一步减小。可以想到,最少可以利用三角函数的14\frac{1}{4}41个周期,通过平移,反转来还原整个周期。即xxx取到π2\frac{\pi}{2}2π,计算此时的项数。

The minimum m is 3

通过减少xxx的取值,将项数从10降低到3项,极大的降低了运算量,大大的提高了效率。

3 结果验证

现在通过编程,用两项泰勒展开,去输出一个实际的正选波,求取信噪比,验证上述推理过程。sinx≈x−x33!+x55!sinx\approx x-\frac{x^3}{3!}+\frac{x^5}{5!}sinxx3!x3+5!x5利用泰勒2项泰勒展开产生正选波的代码如下:

/*
Function: Generate sin sequences use taylor expansion.
*/
/*
Function: Generate sin sequences use taylor expansion.
*/
#include "stdio.h"
#include "stdint.h"

#define PI	3.1415926535897932384626433832795
#define SEQ_LENTN	10000
void main(void)
{
	int32_t		fs,
				f0,
				Amp,
				pha,
				wav_val,
				n;
	int8_t		sign;		
	double		pha0,
				var_tmp;
	FILE		*fp = NULL;
 
	fp = fopen("sinwav.txt", "w+");
	fs=8000;
	f0=330;
	Amp=10000;
	pha0=PI/6;
	pha=0;
	sign=1;
	for(n=0;n<SEQ_LENTN;n++)
	{
		var_tmp=2*PI*pha/fs+pha0;
		if(var_tmp>PI)
		{
			var_tmp-=PI;
			sign=-1;
		}
		else
			sign=1;
		if(var_tmp>PI/2)
			var_tmp=PI-var_tmp;

		wav_val=sign*Amp*(var_tmp-var_tmp*var_tmp*var_tmp/6+\
								var_tmp*var_tmp*var_tmp*var_tmp*var_tmp/120);
		
		pha+=f0;
		if(pha>fs)pha-=fs;
		fprintf(fp, "%d\n",wav_val);
	}
	fclose(fp);
}

产生的正选波如图所示
泰勒方法产生sin

现在编写代码,求解实际的信噪比,跟理论值相比较。

for(n=0;n<SEQ_LENTN;n++)
	{
		var_tmp=2*PI*pha/fs+pha0;
		wav_val_ideal=Amp*sin(var_tmp);
		if(var_tmp>PI)
		{
			var_tmp-=PI;
			sign=-1;
		}
		else
			sign=1;
		if(var_tmp>PI/2)
			var_tmp=PI-var_tmp;

		wav_val_actual=sign*Amp*(var_tmp-var_tmp*var_tmp*var_tmp/6+\
								var_tmp*var_tmp*var_tmp*var_tmp*var_tmp/120);
		signal_power+=wav_val_ideal*wav_val_ideal;
		nosie_power+=(wav_val_ideal-wav_val_actual)*(wav_val_ideal-wav_val_actual);
		pha+=f0;
		if(pha>fs)pha-=fs;
	}
	SNR=signal_power/nosie_power;
	SNR=10*log10(SNR);
	printf("SNR is %.2f dB.\r\n",SNR);

输出结果:

SNR is 54.73 dB.

满足大于40dB的信噪比要求,与理论分析是一致的。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值