我正在從事一個專案，該專案將計算正弦波作為控制回路的輸入。正弦波的頻率為 280Hz，控制回圈每 30us 運行一次，所有內容都用 C 語言撰寫，用于 Arm Cortex M7。目前我們只是在做：

double time;
void control_loop() {
    time  = 30e-6;
    double sine = sin(2 * M_PI * 280 * time);
    ...
}

出現兩個問題/問題：

跑久了time會變大。突然之間，正弦函式的計算時間急劇增加（見圖）。為什么是這樣？這些功能通常是如何實作的？有沒有辦法避免這種情況（沒有明顯的精度損失），因為速度對我們來說是一個重要因素？我們正在使用 math.h (Arm GCC) 中的 sin。
一般如何處理時間：當長時間運行時，變數time不可避免地會達到雙精度的極限。即使使用計數器也time = counter * 30e-6; 只能改善這一點，但不能解決它。由于我當然不是第一個想要長時間生成正弦波的人，因此必須有一些想法/論文/......關于如何快速準確地實作這一點。

uj5u.com熱心網友回復：

不是將正弦計算為時間的函式，而是保持正弦/余弦對并通過復數乘法將其推進。這不需要任何三角函式或查找表；只有四個乘法和偶爾的重新規范化：

static const double dx = cos(2 * M_PI * 280 * 30e-6);
static const double dy = sin(2 * M_PI * 280 * 30e-6);
double x = 1, y = 0; // complex x   iy
int counter = 0;

void control_loop() {
    double xx = dx*x - dy*y;
    double yy = dx*y   dy*x;
    x = xx, y = yy;

    // renormalize once in a while, based on
    // https://www.gamedev.net/forums/topic.asp?topic_id=278849
    if((counter   & 0xff) == 0) {
        double d = 1 - (x*x   y*y - 1)/2;
        x *= d, y *= d;
    }

    double sine = y; // this is your sine
}

如果需要，可以通過重新計算dx、來調整頻率dy。

解釋

有些人要求在評論中解釋為什么會這樣。最簡單的解釋是使用角和三角恒等式：

xx = cos((n 1)*a) = cos(n*a)*cos(a) - sin(n*a)*sin(a) = x*dx - y*dy
yy = sin((n 1)*a) = sin(n*a)*cos(a)   cos(n*a)*sin(a) = y*dx   x*dy

正確性遵循歸納法。

如果我們根據歐拉公式將這些正弦/余弦對視為復數，這本質上就是De Moivre公式。

一種更有洞察力的方法可能是從幾何角度看待它。復數乘法exp(ia)相當于按a弧度旋轉。因此，通過重復乘以dx idy = exp(ia)，我們1 0i沿著單位圓遞增地旋轉我們的起點。在y根據再次歐拉公式坐標，是電流相位的正弦值。

正常化

雖然相位在每次迭代中繼續前進，但由于舍入誤差，幅度（又名范數）會x iy偏離1。然而，我們對生成幅度正弦感興趣1，因此我們需要歸一化x iy以補償數字漂移。當然，直接的方法是將其除以自己的規范：

double d = 1/sqrt(x*x   y*y);
x *= d, y *= d;

這需要計算倒數平方根。即使我們每 X 次迭代只標準化一次，避免它仍然很酷。幸運的|x iy|是已經接近1，因此我們只需要稍微調整一下就可以避免它。展開d大約的運算式1（一階泰勒近似），我們得到代碼中的公式：

d = 1 - (x*x   y*y - 1)/2

TODO：要完全理解這種近似的有效性，需要證明它補償舍入誤差的速度比它們累積的速度更快——從而得到需要應用它的頻率的界限。

uj5u.com熱心網友回復：

如何生成一個可愛的正弦。

DAC 是 12 位，所以你只有 4096 級。每個周期發送超過 4096 個樣本是沒有意義的。在現實生活中，您將需要更少的樣本來生成高質量的波形。

使用查找表創建 C 檔案（使用您的 PC）。將輸出重定向到檔案 ( https://helpdeskgeek.com/how-to/redirect-output-from-command-line-to-text-file/ )。

#define STEP   ((2*M_PI) / 4096.0)

int main(void)
{
    double alpha = 0;

    printf("#include <stdint.h>\nconst uint16_t sine[4096] = {\n");
    for(int x = 0; x < 4096 / 16; x  )
    {
        for(int y = 0; y < 16; y  )
        {
            printf("%d, ", (int)(4095 * (sin(alpha)   1.0) / 2.0));
            alpha  = STEP;
        }
        printf("\n");
    }
    printf("};\n");
}

https://godbolt.org/z/e899d98oW

配置定時器每秒觸發4096*280=1146880次溢位。設定定時器以生成 DAC 觸發事件。對于 180MHz 定時器時鐘，它不會精確，頻率將為 279.906449045Hz。如果您需要更高的精度，請更改樣本數以匹配您的計時器頻率或/和更改計時器時鐘頻率（H7 計時器可運行至 480MHz）
將 DAC 配置為使用 DMA，并在觸發事件時將步驟 1 中創建的查找表中的值傳輸到 DAC。
使用示波器欣賞美麗的正弦波。請注意，您的微控制器內核根本不會加載。您將擁有它用于其他任務。如果你想改變周期簡單地重新配置定時器。您可以根據需要每秒執行多次。要重新配置定時器，請使用定時器 DMA 突發模式 - 這將在更新事件時自動重新加載 PSC 和 ARR 暫存器，而不會干擾生成的波形。

我知道這是高級 STM32 編程，它需要暫存器級編程。我用它在我們的設備中生成復雜的波形。

這是正確的做法。沒有控制回路，沒有計算，沒有核心負載。

uj5u.com熱心網友回復：

該函式可以改寫為

double n;
void control_loop() {
    n  = 1;
    double sine = sin(2 * M_PI * 280 * 30e-6 * n);
    ...
}

這與問題中的代碼完全相同，具有完全相同的問題。但現在可以簡化：

280 * 30e-6 = 280 * 30 / 1000000 = 21 / 2500 = 8.4e-3

這意味著當n達到 2500 時，您已經輸出了 21 個周期的正弦波。這意味著您可以n重新設定為 0。生成的代碼是：

int n;
void control_loop() {
    n  = 1;
    if (n == 2500)
        n = 0;
    double sine = sin(2 * M_PI * 8.4e-3 * n);
    ...
}

只要您的代碼可以無問題地運行 21 個周期，它就會永遠運行而不會出現問題。

uj5u.com熱心網友回復：

我對現有的答案感到震驚。你發現的第一個問題很容易解決，當你解決第一個問題時，下一個問題就神奇地消失了。

您需要對數學有基本的了解才能了解它是如何作業的。回想一下，sin(x 2pi)在sin(x)數學上只是。當您的sin(float)實作切換到另一種演算法時，您看到的時間會大幅增加，而您確實想避免這種情況。

請記住，float它只有 6 位有效數字。100000.0f*M_PI x將這 6 位數字用于100000.0f*M_PI，因此沒有任何剩余x。

因此，最簡單的解決方案是跟蹤x自己。在t=0您初始化x為0.0f. 每 30 us，您增加x = M_PI * 280 * 30e-06;. 這個公式中沒有出現時間！最后，如果x>2*M_PI，你遞減x-=2*M_PI;(Since sin(x)==sin(x-2*pi)

您現在有一個x很好地保持在0to范圍內的3.1417，其中sin速度很快并且 6 位精度都很有用。

uj5u.com熱心網友回復：

正如一些評論中所指出的，time價值隨著時間的推移不斷增長。這帶來了兩個問題：

該sin函式可能必須在內部執行模數才能使內部值進入支持的范圍。
由于增加了更高的數字，時間的解析度會隨著值的增加而變得越來越差。

進行以下更改應該可以提高性能：

double time;
void control_loop() {
    time  = 30.0e-6;
    if((1.0/280.0) < time)
    {
        time -= 1.0/280.0;
    }
    
    double sine = sin(2 * M_PI * 280 * time);
    ...
}

請注意，一旦進行此更改，您將不再擁有時間變數。

uj5u.com熱心網友回復：

使用查找表。您在與 Eugene Sh. 的討論中的評論：

與正弦頻率（如 280.1Hz）的小偏差就可以了。

在這種情況下，控制間隔為 30 μs，如果您有一個包含 119 個樣本的表格，您一遍又一遍地重復，您將得到 280.112 Hz 的正弦波。由于您有一個 12 位 DAC，如果您將其直接輸出到 DAC，您只需要 119 * 2 = 238 個位元組來存盤它。如果您將其用作進一步計算的輸入，就像您在評論中提到的那樣，您可以將其存盤為float或double根據需要。在帶有嵌入式靜態 RAM 的 MCU 上，從記憶體加載最多只需要幾個周期。

uj5u.com熱心網友回復：

我想直接解決您代碼中的嵌入式編程問題 - @0____________ 的答案是在微控制器上執行此操作的正確方法，我不會重述相同的內容。

代表時間變數應該永遠是浮點數。如果您的增量不是 2 的冪，則錯誤將始終累積。即使是這樣，最終您的增量將小于最小增量并且計時器將停止。始終使用整數表示時間。你可以選擇一個足夠大的整數來忽略翻轉——一個代表毫秒的無符號 32 位整數將需要 50 天來翻轉，而一個無符號 64 位整數將需要超過 5 億年。
在不關心信號相位的情況下生成任何周期信號不需要時間變數。相反，您可以保留一個在周期結束時重置為 0 的內部計數器。（當您將 DMA 與查找表一起使用時，這正是您所做的 - 計數器是 DMA 控制器的下一個讀取指標。）
每當您在微控制器中使用超越函式（例如正弦）時，您的第一個想法應該是“我可以為此使用查找表嗎？” 您無法享受現代作業系統在 4GHz 多核處理器上優化調整負載的樂趣。您經常處理單個執行緒，該執行緒將等待您的 200MHz uC 使 FPU 退出待機狀態并執行近似演算法。超驗函式的成本很高。LUT 也有成本，但如果您不斷地使用該功能，您很有可能會更喜歡 LUT 的權衡。

uj5u.com熱心網友回復：

如果您有幾千位元組的可用記憶體，則可以使用查找表完全消除此問題。

采樣周期為 30 μs 時，2500 個樣本的總持續時間為 75 ms。這恰好等于 280 Hz 下 21 個周期的持續時間。

我還沒有測驗或編譯以下代碼，但它至少應該演示該方法：

double sin2500() {
  static double *table = NULL;
  static int n = 2499;
  if (!table) {
    table = malloc(2500 * sizeof(double));
    for (int i=0; i<2500; i  ) table[i] = sin(2 * M_PI * 280 * i * 30e-06);
  }
  n = (n 1) % 2500;
  return table[n];
}

uj5u.com熱心網友回復：

其他人基于模的概念的變體怎么樣：

int t = 0;
int divisor = 1000000;
void control_loop() {
    t  = 30 * 280;
    if (t > divisor) t -= divisor;
    double sine = sin(2 * M_PI * t / (double)divisor));
    ...
}

它以整數計算模，然后不會導致舍入誤差。

uj5u.com熱心網友回復：

我認為可以使用模數，因為它sin()是周期性的。然后你不必擔心這些問題。

double time = 0;
long unsigned int timesteps = 0;
double sine;
void controll_loop()
{
  timesteps  ;
  time  = 30e-6;
  if( time > 1 )
  {
    time -= 1
  }
  sine = sin( 2 * M_PI * 280 * time );
  ...
}

uj5u.com熱心網友回復：

有一種替代方法可以計算增加一些非常小的角度的一系列正弦（和余弦）值。它本質上是計算一個圓的 X 和 Y 坐標，然后將 Y 值除以某個常數以生成正弦，然后將 X 值除以相同的常數以生成余弦。
如果您滿足于生成一個“非常圓的橢圓”，您可以使用以下 hack，這歸功于 1960 年代的 Marvin Minsky。它比計算正弦和余弦快得多，盡管它在系列中引入了非常小的誤差。以下是 Hakmem 檔案第 149 項的摘錄。概述了明斯基圓演算法。

專案 149（明斯基）：圓演算法這是在點繪圖顯示幕上繪制幾乎圓形的優雅方法：

NEW X = OLD X - epsilon * OLD Y
NEW Y = OLD Y   epsilon * NEW(!) X

這將形成一個以原點為中心的非常圓的橢圓，其大小由初始點決定。epsilon 決定回圈點的角速度，對偏心率有輕微影響。如果 epsilon 是 2 的冪，那么我們甚至不需要乘法，更不用說平方根、正弦和余弦了！“圓”將非常穩定，因為這些點很快就會變成周期性的。

當我試圖在顯示黑客中保存一個暫存器時，錯誤地發明了圓演算法！Ben Gurley 只使用了大約六到七個指令就實作了驚人的顯示技巧，這是一個了不起的奇跡。但它基本上是面向行的。我突然想到擁有曲線會令人興奮，我試圖用最少的指令獲得曲線顯示技巧。

這是 hakmem 的鏈接：http ://inwap.com/pdp10/hbaker/hakmem/hacks.html

轉載請註明出處，本文鏈接：https://www.uj5u.com/qukuanlian/338880.html

標籤：C 表现时间精确三角学

上一篇：如何在一行上讀取多個浮點數，然后將它們添加到陣列中？

下一篇：結構體中字符指標的記憶體分配

C中的無限正弦生成

解釋

正常化