Re: [問題] 還有人記得傅立葉變換嗎？

[HuangJC (吹笛牧童)]

※ 引述《isaacting (2312312)》之銘言：
: 其實網路上有很多資源可以使用
: 1. FFTW
: ttp://www.fftw.org/
: 就一般開發者的狀況，不可能有人寫得比這個library更好更快
: 除非你是專業的軟體工程師，而且還懂一堆硬體的optimization instructions set
: 2. Octave
: 他其實就是免費的Matlab啦, 對於fft這個指令也有內建
: Octave也是使用 fftw來進行FFT的運算
: 如果單純是要開發應用程式的話 ，會建議使用 fftw，這會讓你省去一堆麻煩
: 然後你也可以用 Octave來進行驗證
: 當然，如果你想知道如何使用fftw的話，歡迎站內信
: 我一般是用Linux在做訊號處理的開發

我在 win 上開發，已經成功用上了

1。速度：快逾 20 倍

我手上的資料是五萬六千點取樣
同事的雖然陣列只開最大 512，不過擴充不難

擴充完後測，他的程式 430 ns
這個 DLL 15 ns
快了 20倍以上


2。支援任意 size,不必 2的 N次方

因為他那是 FFT，必需接受 2的N次方的陣列
所以我都必需準備 size = 65536, 在五萬六千點之後補 0
這個根本不精確！
而這個 DLL 強大，可以使用 arbitrary input size
我試著處理 s = 1 + sin(x) + 2 * sin(2*x) 這樣的訊號源來做測試
其中 sin(x) 就是 1hz, 假設我用 10hz 取樣
(按照取樣定理，我這訊號最大是 2hz，所以我必需用 4hz 以上取樣，10hz 很足夠了）
取 10 個點，用 FFT 卻必需準備大小為 16 的陣列
算完根本誤差極大
若取 10 個點，用這個 DLL，它可以只用 size=10 的陣列，果然結果就精準算出


3. 32bit 版本可以用， 64 bit 程式會當掉，
   不知是否和我用模擬器而非真正的 win 有關
（我的開發環境是 Mac OS + 模擬器跑 win）

其他我還在研究

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 223.138.199.136 (臺灣)
※ 文章網址: https://www.ptt.cc/bbs/Math/M.1573674884.A.F78.html
※ 編輯: HuangJC (223.138.199.136 臺灣), 11/14/2019 03:56:05

p = fftw_plan_dft_1d(FFT_N, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);

這是照抄手冊寫出的程式
我搞不懂為什麼要繞這樣一圈來寫
先建立 plan, 再執行 plan
它提供了什麼彈性？
有必要嗎？


※ 編輯: HuangJC (223.138.199.136 臺灣), 11/14/2019 04:11:01
→ HuangJC     : 好吧，或許有，因為我要連做 20次，所以用到了 XD
→ hsnuyi      : 你有修過DSP跟數值分析/穩定度吧?
沒有，
https://zh.wikipedia.org/wiki/%E6%95%B0%E5%80%BC%E7%A8%B3%E5%AE%9A%E6%80%A7
是指這個嗎？

我是學過'用逼進法求解，有可能發散'，不知是否指這回事
但還是不明白這裡為什麼要提這件事

→ recorriendo : 是誤差?還是頻域取樣的問題? #1SDMTp7K (Math)
好，來談這個問題
首先我雖然是理科的學生，但畢業很久，都還老師了
而且我沒學到的科目，用的可能不是精確的字眼
對我來說，誤差或是頻域取樣的問題
就是同一個問題 :P

那一篇我早看過了，頻率洩漏本來就會發生，那還能指望什麼
再來，從邏輯上我也要說
如果我想表達兩個不同的波型
但我餵給程式的是完全相同的資料
那我能指望傳回的結果有何不同？

舉 1hz 的 sin 波取樣 10點為例
如果我塞入 16 點的陣列裡，後面補 0

和另一種圖形
前面 10/16 是 sin 波，後面 6/16 是 0
但這前面 sin 波部份是一秒！
我的波就是真的長這樣嘛
然後我做 16 點取樣

這樣傳入的陣列資料，不是一模模一樣樣？
所以事實上我可以解釋成
"我正在用一種新的波型來計算，希望算到相近的答案"
答案不可能一樣，波型都變了

--------------

我們來講另一個可能好了，內插法，這同事教我的
方法是：假設輸入波型以 10hz 取樣
（今天我們要談的就是取樣設備有限制，它就是運作在 10hz
  這也是沒法子的事，如果可以運作在 16hz 那大家都沒事，事情很簡單）

好，那我就可以取到 10 個點的值
接著，我不要填 10 個值，再補 6 個 0
我用內插法把這 10 個值變成 16 個值，如何？

程式我直接列出

#define M_PI            3.14159265358979323846
#define SAMPLE_RATE     10      //      hz
#define FFT_N           16

double sample[SAMPLE_RATE];

void test()
{
    fftw_complex *in;
    fftw_plan p;

    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * FFT_N);
    memset(in, 0, sizeof(fftw_complex) * FFT_N);

    //fill data
    double  deltaDeg = 2 * M_PI / SAMPLE_RATE;
    int     i;

    //這裡是 10 個點的取樣
    for(i = 0; i < SAMPLE_RATE; i++) {
        double x = deltaDeg * i;
        sample[i] = sin(x)；  //很簡單，我只是輸入 sin 波
    }

    //這裡是把 10 點個變成 16 個點的內插法
    for(i = 0; i < FFT_N; i++) {
        double x = (double)i / FFT_N * 10;
        int x1 = (int)x;
        int x2 = x1 + 1;
        in[i][0] = sample[x1] * (x - x1) + sample[x2] * (x2 - x);
        in[i][1] = 0;
    }

    p = fftw_plan_dft_1d(FFT_N, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);

    //read data
    TRACE("output:\n");
    for(i = 0; i < FFT_N; i++) {
        double real = in[i][0];
        double imag = in[i][1];
        if ( i < FFT_N / 2 )
            TRACE("s[%d]=(%f,%f) = %f freq %f\n",
            i, real, imag, sqrt(real * real + imag * imag), (double) i );
    }

    fftw_free(in);
}

※ 編輯: HuangJC (223.138.199.136 臺灣), 11/16/2019 02:47:27

好，請注意到這是用內插法把 10 點取樣擴充成 16 點
完全不是以往學的 補 0
那我們可以讀到什麼結果？

※ 編輯: HuangJC (223.138.199.136 臺灣), 11/16/2019 02:51:09

output:
s[0]=(0.000000,0.000000) = 0.000000 freq 0.000000
s[1]=(1.175571,-7.012620) = 7.110471 freq 1.000000
s[2]=(0.000000,-0.000000) = 0.000000 freq 2.000000
s[3]=(1.175571,-0.376322) = 1.234336 freq 3.000000
s[4]=(0.000000,0.000000) = 0.000000 freq 4.000000
s[5]=(1.175571,0.685697) = 1.360936 freq 5.000000
s[6]=(0.000000,-0.000000) = 0.000000 freq 6.000000
s[7]=(1.175571,1.657851) = 2.032348 freq 7.000000

這結果挺好的
所以，不要再補 0 啦
用內插法好了 XD

還是我附上補 0 的結果？

output:
s[0]=(0.000000,0.000000) = 0.000000 freq 0.000000
s[1]=(4.367883,-1.809236) = 4.727762 freq 0.625000
s[2]=(-2.883839,-2.883839) = 4.078364 freq 1.250000
s[3]=(-0.201906,0.487443) = 0.527605 freq 1.875000
s[4]=(-0.726543,0.000000) = 0.726543 freq 2.500000
s[5]=(-0.072232,-0.174384) = 0.188752 freq 3.125000
s[6]=(-0.193845,0.193845) = 0.274138 freq 3.750000
s[7]=(-0.289519,-0.119923) = 0.313373 freq 4.375000

有覺得比較好嗎？
因為頻率間隔的計算是 10hz/16 = 0.625
所以不會精確的有 1hz 這個讀值，判讀變得很麻煩
必需從最接近的 0.625hz & 1.25 hz 去內插

※ 編輯: HuangJC (223.138.199.136 臺灣), 11/16/2019 02:52:37

其實我還有個大膽的想法
既然 FFT 是運作在週期波上，那我自己把取樣到的值重覆 N 次不就好了？
因為 10hz, 我取樣到 10 個值嘛

那我重覆三次，不就有 30個值？
把它塞入 32個值的陣列裡，不是補的 0 就比較少了？
本來 10 塞入 16 要補 6 個零，有效資料是 10/16 = 63%
現在 30 塞入 32 只補 2 個零，有效資料是 30/32 = 94%

有沒有想過這樣會比較精確？

這是結果：
output:
s[0]=(0.000000,0.000000) = 0.000000 freq 0.000000
s[1]=(0.130241,-0.654765) = 0.667592 freq 0.312500
s[2]=(0.749410,-1.809236) = 1.958303 freq 0.625000
s[3]=(8.080340,-12.093083) = 14.544228 freq 0.937500
s[4]=(-2.883839,2.883839) = 4.078364 freq 1.250000
s[5]=(-1.603337,1.071315) = 1.928317 freq 1.562500
s[6]=(-1.176792,0.487443) = 1.273751 freq 1.875000
s[7]=(-0.920980,0.183194) = 0.939023 freq 2.187500
s[8]=(-0.726543,0.000000) = 0.726543 freq 2.500000
s[9]=(-0.563101,-0.112008) = 0.574133 freq 2.812500
s[10]=(-0.421000,-0.174384) = 0.455687 freq 3.125000
s[11]=(-0.297790,-0.198977) = 0.358149 freq 3.437500
s[12]=(-0.193845,-0.193845) = 0.274138 freq 3.750000
s[13]=(-0.110592,-0.165513) = 0.199060 freq 4.062500
s[14]=(-0.049674,-0.119923) = 0.129803 freq 4.375000
s[15]=(-0.012499,-0.062838) = 0.064069 freq 4.687500

看到最接近的 0.93 那個頻率，這結果還不錯 XD
但是無論如何，不用湊 2的 N次方，還是快樂多了

只有 10個取樣點，那就用 10 個值的陣列
結果會是

output:
s[0]=(0.000000,0.000000) = 0.000000 freq 0.000000
s[1]=(-0.000000,-5.000000) = 5.000000 freq 1.000000
s[2]=(0.000000,-0.000000) = 0.000000 freq 2.000000
s[3]=(0.000000,-0.000000) = 0.000000 freq 3.000000
s[4]=(0.000000,-0.000000) = 0.000000 freq 4.000000

有沒有，精準的有 1hz 那個讀值，而且沒有頻率洩漏
所以我真喜歡這個副程式 XD
麻省理工就是神～

本來快速傅立葉就是為了實用而產生的，就是快
相對的慢速傅立葉（其實沒這個名字；但相對就慢啊）應該慢多了
結果用這東西可以輸入任意 size 的資料
而且還是比其他用快速傅立葉的副程式快了二十倍
又快又好用

※ 編輯: HuangJC (223.138.199.136 臺灣), 11/16/2019 03:12:29
→ HuangJC     : 我能用，我同事不能用，因為沒有他的平台副程式..
→ HuangJC     : 對了，我的內插還只是線性一階內插..
→ HuangJC     : 如果有專屬的 sin 波內插，那可漂亮了
→ recorriendo : 一樣的輸入 輸出卻不同 就是你程式設計者本身的問題
→ recorriendo : 這不是什麼誤差 而是你不願去了解數學細節 曲解輸出
→ recorriendo : 的意義  更不是頻率洩漏 那完全是另一個問題 跟補不
→ recorriendo : 補零無關
→ recorriendo : 你可以直接補到128點,256點等等 看看得出來結果怎麼
→ recorriendo : 變化 或許有助了解其中原理
→ wohtp       : 對一個isolated, well-defined wave packet來說，
→ wohtp       : zero padding沒什麼差，反正你真的繼續取樣下去也是
→ wohtp       : 繼續拿到一串零。
→ wohtp       : 除此之外，zero padding當然是在惡搞輸入訊號，最後
→ wohtp       : 你會拿到一樣的結果才有鬼
> 反正你真的繼續取樣下去也是繼續拿到一串零。

才不是咧，我在想良好的表達都是要附圖的
最近我有想把圖畫出來講
（反之沒附圖而回應我的人，覺得我能聽懂？聽不懂的話溝通障礙要怪在我身上嗎？）
我不像熟 matlab 的人可以輕輕鬆鬆產出畫面
我都必需一行一行建制完整程式

有時我會想，我已經有自己的解決方案，可以停手了
但既然有人陪我討論，那這熱情值得鼓勵
我想法子生圖就是了

> 除此之外，zero padding當然是在惡搞輸入訊號，最後你會拿到一樣的結果才有鬼

怪的是這就是我要表達的啊，所以你應該懂我在說什麼...
那回過頭來說吧，竟然有人說錯的是我，因為程式是我寫的
Orz.. 而且還不只一個
我相信是溝通上出了問題，唉..

※ 編輯: HuangJC (223.141.233.39 臺灣), 11/20/2019 02:12:04
→ wohtp       : 如果你的波形不是一個小包包就不能補零
→ wohtp       : 差別在哪你自己都已經舉例出來了
→ HuangJC     : 小包包是指？我聽不懂