问题
三角函数 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的信噪比要求,与理论分析是一致的。