2010-10-19 89 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/