Dsp相關問題思考筆記（一）

問題

三角函數 s i n x sinx sinx可以通過泰勒展開寫成形式: s i n x = ∑ n = 0 + ∞ ( − 1 ) n x 2 n + 1 ( 2 n + 1 ) ! sinx=\sum_{n=0}^{+\infty}{(-1)^{n}\frac{x^{2n+1}}{(2n+1)!}} sinx=n=0∑+∞​(−1)n(2n+1)!x2n+1​工程上可以根據此公式，進行三角函數的近似計算，進而輸出一個正選波。但是展開的階數不可能取到無限項，應在滿足信噪比的情況下，取到盡可能小的有限項。即用前m項進行近似： s i n x = ∑ n = 0 m − 1 ( − 1 ) n x 2 n + 1 ( 2 n + 1 ) ! + e r r o r sinx=\sum_{n=0}^{m-1}{(-1)^{n}\frac{x^{2n+1}}{(2n+1)!}}+error sinx=n=0∑m−1​(−1)n(2n+1)!x2n+1​+error現在考慮，在滿足信噪比不低于40dB的前提下，滿足的m最小值。

思考過程

1基本

error實際上是前一項代數式的高階無窮小，具體為形式： e r r o r = x 2 m + 1 ( 2 m + 1 ) ! − x 2 m + 3 ( 2 m + 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項決定，即 e r r o r ≈ x 2 m + 1 ( 2 m + 1 ) ! error\approx\frac{x^{2m+1}}{(2m+1)!} error≈(2m+1)!x2m+1​則信号的信噪比有 S N R = s i g n a l p o w e r n o s i e p o w e r = 1 2 [ x 2 m + 1 ( 2 m + 1 ) ! ] 2 ≥ 1 0 4 SNR=\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​]221​​≥104考慮最惡劣情況，則 x x x取區間内的最大值，則信噪比最壞。正常輸出 s i n x sinx sinx時，為避免累計誤差，都會利用其周期性，即 x x x的最大值為 2 π 2\pi 2π，則公式變為 1 2 [ 2 π 2 m + 1 ( 2 m + 1 ) ! ] 2 ≥ 1 0 4 \frac{\frac{1}{2}}{[\frac{2\pi^{2m+1}}{(2m+1)!}]^{2}}\ge10^{4} [(2m+1)!2π2m+1​]221​​≥104此公式為一個非線性方程，但可以知道，随着階數 m m m的增大，信噪比是不斷提升的，求解非線性方程的經典算法是二分法。編寫求解該非線性方程的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 進階

s i n x sinx sinx中 x x x取到 2 π 2\pi 2π時，滿足信噪比最少需要取到10項。如果想要在保證信噪比前提下進一步減少項數，則要讓 x x x取值進一步減小。可以想到，最少可以利用三角函數的 1 4 \frac{1}{4} 41​個周期，通過平移，反轉來還原整個周期。即 x x x取到 π 2 \frac{\pi}{2} 2π​,計算此時的項數。

The minimum m is 3

通過減少 x x x的取值，将項數從10降低到3項，極大的降低了運算量，大大的提高了效率。

3 結果驗證

現在通過程式設計，用兩項泰勒展開，去輸出一個實際的正選波，求取信噪比，驗證上述推理過程。 s i n x ≈ x − x 3 3 ! + x 5 5 ! sinx\approx x-\frac{x^3}{3!}+\frac{x^5}{5!} sinx≈x−3!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);
}
           

産生的正選波如圖所示

現在編寫代碼，求解實際的信噪比，跟理論值相比較。

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的信噪比要求，與理論分析是一緻的。