2010-10-19 87 views
1

我想为octave写一个.oct函数,给定一个正弦波值,在-1和1之间,正弦波周期返回一个周期长度的正弦波矢量,矢量中的最后一个值是给定的正弦波值。我的代码到目前为止是:递归创建一个正弦波给定一个单一的正弦波值和周期

#include <octave/oct.h> 
#include <octave/dColVector.h> 
#include <math.h> 
#define PI 3.14159265 

DEFUN_DLD (sinewave_recreate, args, , "args(0) sinewave value, args(1) is period") 
{ 
octave_value_list retval; 

double sinewave_value = args(0).double_value(); 
double period = args(1).double_value(); 
ColumnVector output_sinewave(period);     
double degrees_inc = 360/period; 
double output_sinewave_degrees; 

output_sinewave_degrees = asin(sinewave_value) * 180/PI; 
output_sinewave(period-1) = sin(output_sinewave_degrees * PI/180); 

for (octave_idx_type ii (1); ii < period; ii++) // Start the loop 
    { 
    output_sinewave_degrees = output_sinewave_degrees - degrees_inc; 

    if (output_sinewave_degrees < 0) 
    { 
    output_sinewave_degrees += 360 ; 
    } 

    output_sinewave(period-1-ii) = sin(output_sinewave_degrees * PI/180); 
    } 

retval(0) = output_sinewave;               

return retval;                   
} 

但给补丁结果。我的意思是,它有时会相当精确地重新创建正弦波,而其他时间则会消失。我只是通过创建一个给定的正弦波来确定这一点,将时间的最后一个值填入函数中,通过时间向后重新创建正弦波,然后比较两者的曲线。显然我做错了什么,但我似乎无法确定什么。

+0

http://pizer.wordpress.com/2010/02/08/fast-digital-sine-oscillator/ – sellibitze 2010-10-19 05:37:00

+0

一般来说,任何给定值都有两个答案,具体取决于正弦波在经过给定值时是上升还是下降。 – 2010-10-20 01:23:53

回答

4

让我们先从一些三角恒等式:

sin(x)^2 + cos(x)^2 == 1 
sin(x+y) == sin(x)*cos(y) + sin(y)*cos(x) 
cos(x+y) == cos(x)*cos(y) - sin(x)*sin(y) 

鉴于目前点x正弦和余弦,我们可以准确地计算出大小d的步骤后的值,预先计算sd = sin(d)cd = cos(d)后:

sin(x+d) = sin(x)*cd + cos(x)*sd 
cos(x+d) = cos(x)*cd - sin(x)*sd 

给定初始正弦值,可以计算初始余弦值:

cos(x) = sqrt(1 - sin(x)^2) 

请注意,有两种可能的解决方案,对应于两个可能的平方根值。还要注意,这些身份中的所有角度均以弧度表示,如果您要返回波形,则d需要为负数。

+1

它[完美](http://ideone.com/OzKu5)。 – Keynslug 2010-10-19 15:01:55

0

你可能会尝试一个更简单的方法来通过。 只是记得,如果

y = sin(x) 

y然后一阶导数将等于

dy/dx = cos(x) 

因此,在计算的每一步你加入的y一些增量等于当前值

dy = cos(x) * dx 

但是这可能会降低您的准确度作为一个副作用。你可以探究它。 HTH。


似乎略有改善方程往往更准确:

dy = cos(x + dx/2) * dx 

this看看。

+0

这避免了反正弦计算,这可能是错误来自的地方。但是,有一种方法可以精确地求解它,而不是通过近似计算,而不必为每个步骤计算余弦。 – 2010-10-19 11:44:57

1

Mike指出cos(x)有两种可能的解决方案让我意识到我需要解决正弦波的相位模糊问题。我的第二个成功的尝试是:

#include <octave/oct.h> 
#include <octave/dColVector.h> 
#include <math.h> 
#define PI 3.14159265 

DEFUN_DLD (sinewave_recreate_3, args, , "args(0) sinewave value, args(1) is period, args(2) is the phase") 
{ 
octave_value_list retval; 

double sinewave_value = args(0).double_value(); 
double period = args(1).double_value(); 
double phase = args(2).double_value(); 
ColumnVector output_sinewave(period);     
double X0 = asin(sinewave_value); 

if (sinewave_value < 0 & phase > 180 & phase < 270) 
{ 
X0 = PI + (0 - X0); 
} 

if (sinewave_value < 0 & phase >= 270) 
{ 
X0 = X0 + 2 * PI; 
} 

if (sinewave_value > 0 & phase > 90) 
{ 
X0 = PI - X0; 
} 

if (sinewave_value > 0 & phase < 0) 
{ 
X0 = X0 + PI/2; 
} 

double dx = PI/180 * (360/period); 

for (octave_idx_type ii (0); ii < period; ii++) // Start the loop 
{ 
output_sinewave(period-1-ii) = sin(X0 - dx * ii); 
} 

retval(0) = output_sinewave;               

return retval;                   
} 

感谢也是由于Keynslug。

1

有一个简单的公式。这里是在Python的例子:

import math 
import numpy as np 

# We are supposing step is equal to 1degree 
T = math.radians(1.0/360.0) 
PrevBeforePrevValue = np.sin(math.radians(49.0))   # y(t-2) 
PrevValue = np.sin(math.radians(50.0))      # y(t-1) 

ValueNowRecursiveFormula = ((2.0*(4.0-T*T))/(4.0+T*T))*PrevValue - PrevBeforePrevValue 

print("From RECURSIVE formula - " + str(ValueNowRecursiveFormula)) 

的细节可以在这里找到: http://howtodoit.com.ua/en/on-the-way-of-developing-recursive-sinewave-generator/