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

[HuangJC (吹笛牧童)]

: 好，來談這個問題
: 首先我雖然是理科的學生，但畢業很久，都還老師了
: 而且我沒學到的科目，用的可能不是精確的字眼
: 對我來說，誤差或是頻域取樣的問題
: 就是同一個問題 :P
:
: 那一篇我早看過了，頻率洩漏本來就會發生，那還能指望什麼
: 再來，從邏輯上我也要說
: 如果我想表達兩個不同的波型
: 但我餵給程式的是完全相同的資料
: 那我能指望傳回的結果有何不同？
:
: 舉 1hz 的 sin 波取樣 10點為例

我用了一個網頁
https://www.desmos.com/calculator
用這個來幫我產生圖形，不用寫程式

現在我們來繪 y=sin(x)

https://imgur.com/a/Sb6znAc
見圖

圖中有兩條曲線，上面那條就是 y = sin(x)
粉紅色的區間裡就是 0 < x < 2PI
如果我另外定義粉紅色區間為 一秒，那這圖形就是 1hz

注意到下面那條曲線，它在粉紅色區間裡和上面那條曲線一模一樣

https://imgur.com/a/XRUXlTv
見圖

這條曲線怎麼畫的呢？
它就是我把 sin(x) 取傅立葉轉換，得出的係數，再寫成反轉換形式重建繪出的
因為兩條曲線可以完全重合
所以證明我的假想無誤

但是要做 FFT，有些條件要設定；取樣頻率啊什麼的
所以我交代一下：我是用 10hz 取樣，在一秒裡取了 10 個點
取完後，塞入 size = 16 的陣列裡，所以後面補了 6 個 0

而這樣事實上我所表達的根本是另一種圖形
這圖形其實是先有一秒的 sin 波，緊接著一段 y=0 的直線
其中 sin 波佔 10/16, y=0 佔 6/10, 加起來就是 16/16, 一整個週期


這上下兩個不同的波型，在 10 hz 取樣裡，根本是一模一樣的輸入值
不然你告訴我要怎麼輸入

input:
s[0]=0.000000
s[1]=0.587785
s[2]=0.951057
s[3]=0.951057
s[4]=0.587785
s[5]=0.000000
s[6]=-0.587785
s[7]=-0.951057
s[8]=-0.951057
s[9]=-0.587785

而既然我們用一樣的輸入，就不該有不同的輸出

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

如果我計算錯誤，那我不可能把這兩個波型繪得這麼像

這裡是我繪圖的一些輸入
https://imgur.com/a/FPpeDMI

可以看到，我把 sin(x) 輸成 sin(x)+5
因為這樣我才把可以把圖往上移，做成兩個波比對
而 0<x<2PI, 或者 0<x<2PI*16/10, 這是一種區間的繪圖工具
滑鼠去按一下前面，就可以 enable or disable 這個區間
方便看出一個週期

比較有趣的是振幅必需調整
\frac{\left(4.367883\cos0.625x+1.809236\sin0.625x-2.883839\cos1.25x+
2.883839\sin1.25x-0.201906\cos1.875x-0.487443\sin1.875x-0.726543\cos2.5x-
0.072232\cos3.125x+0.174384\sin3.125x-0.193845\cos3.75x-0.193845\sin3.75x
-0.289519\cos4.375x+0.119923\sin4.375x\right)}{7.852}

這是直接拷下來的式子，啊..別理會什麼 \,  \frac 之類代碼
總之這串接成一行，直接貼上去就可以繪出圖形
這串就是看懂傳回陣列後把值塞回去做反轉換
但振幅必需除以 7.852
至於為什麼是 7.852，我也不知道
但這樣可以完美重現原波型就是了
（在取樣點都有通過，因為只做 10hz 取樣，所以夠精密了）



https://imgur.com/nRSMvHT

這是另一張圖，標出了取樣點
受限於取樣密度，所以我重建的函數只會在取樣點和原函數重合
淺色區域就是 y=sin(x), 而深色區域是 y=sin(x) 再加一小段 y=0
以我的認知，所謂的補 0 做法
就是這兩條曲線，會輸入一模一樣的值去做 FFT
既然輸入一模一樣，輸出當然也一模一樣

如果我們把兩個波重疊，把圖放大來看
其中 sin(x) 的部份並不能完全重疊
因為後面補了 6 個 0 啊
sin 波部份還要一模一樣也太強人所難了；就接近而已

再來看到解析的頻率怎麼解釋？

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

這段就重貼前面而已
這其實是 1hz 的 sin 波，可是頻率解析度 = 10hz / 16 = 0.625
而 0.625 重建回來的就是這樣，有 0.625hz, 有 1.25 hz
但就是沒有 1hz 這一項
而重建回的波型，也算有準

這若不叫頻率洩露，也不叫誤差，該叫什麼我們再討論一下
而如果要我把輸入陣列補到 size = 128
也麻煩教一下我怎麼補
畢竟我是 10hz 取樣的

如果說把波形重覆 12 次，那就是 120筆資料，再補 8 個 0
如果說 8 個點不要補零，而是繼續向第 13次重覆的波型取樣
那也行！ XD
但我得說，我認為這是另一個波型
這是怎樣的波型我也可以不厭其煩的畫出來給你
但真的就是另一個波型

所以這明明就是用另一個波型在求 sin(x) 的 FFT 轉換
在我認為，會一模一樣才奇怪
只能說是一種輔助
而愈是接近答案的輔助，愈好
（但這又不叫誤差，所以不能說'誤差愈小'；該怎麼表達？）

當然如果用 16 hz 取樣，那什麼事都沒有
因為這只是個 1hz 的 sin 波，用兩倍以上取樣就已足夠
而且會很精準的重建回來（因為太單純，沒其他頻率了）
可是今天的難題就是：我就偏偏 10hz 取樣，這零要怎麼補？

OK，我認為就是會有這種事

至於說快速傅立葉並沒規定輸入陣列是 2的N次方
或許在數學上是的
但在程式設計上，並不是
因為副程式是使用別人寫好，現成的
人家已經給了這個限制

當我跟同事說輸入陣列不必是 2的N次方時，我被同事噓回來，換他說我不懂
拜託不要這樣，兩邊各有各的堅持
其實，誰會花這麼多時間，寫程式繪圖和你討論？
我們採用這種副程式，已經是業界的共識了
（除非用 fftw, 也就是傅立葉轉換，但可以任意 size，不必 2的 N次方
  老實說它為什麼無此限制，速度又快，我們也不知道
  但討論怎麼補 0這個題目時，就先假設不知道 fftw 吧..）
同事只丟給我一句：這副程式是 chip 廠商提供的，內建加速
那就得用啊..

我們剩下的只有主程式可以變化，要檢討都在主程式檢討
那麼，不這麼補 0 ，不然怎麼補？
不叫做誤差，不然叫什麼..


再來附上內插法的圖形

https://imgur.com/a/wykbB3t

上面是 sin(x), 下面是用 10 hz 取樣，取得 10 個點
然後再內插變成 16 個點

附上輸入陣列的程式

#define SAMPLE_RATE     10      //      hz
#define FFT_N           16

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

        TRACE("input:\n");
        for(i = 0; i < SAMPLE_RATE; i++) {
                double x = deltaDeg * i;
                sample[i] = sin(x);
        }

        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;
        }

這個方法不用補 0，還原的波型也不漂亮
但是主頻率抓得好

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


→ recorriendo : 搞懂人家給你的是0.625hz的值就對了 不懂人家給你的
→ recorriendo : 是什麼就指說有誤差是你本來自己的問題
→ recorriendo : 給你一公斤的東西你說台斤秤上顯示1.666斤這叫誤差?
我聽得懂你的意思，所以那我們就別用'誤差'這個詞
但還是得有一個詞來溝通吧？
請問那個詞是什麼？
謝謝

→ recorriendo : 還有頻率洩漏跟你用幾點做FFT本來就無關 要減少洩漏
既然不是頻率洩漏，那也就不用減少洩漏啊！

→ recorriendo : 請去搜尋window function
→ recorriendo : 專有名詞不是你愛怎樣定義就隨你定義
→ j0958322080 : 反正他高興就好了，他應該認為自己懂了
話不是這樣說，你們也說我錯，我同事也說我錯
我就像乒乓球一樣被踢來踢去
要是說，有才華的人總是有脾氣
我忍受完脾氣就有收獲，那也就值了
千萬不要什麼都得不到
那我當然認為自己已經懂了，止血
這難道不是對我最有利的做法？
上面大哥丟出一句'window function',那就是啦
有點收獲我就很感激了
（還有後面一篇回文；我再想想怎麼去吸收）

推 sunev       : 這個板主要是討論數學，程式方面的限制請洽專板
戰文時才會拿這個當擋箭牌，如果專版能討論得更深入，那我會去專板
現在講究 PI 型人才（跨領域，站兩隻腳的意思）
數學當然為軟體服務
板上如果認為我洗板，那我可以回文時附加在自己文章之後，不增加文章
如果認為這不算數學，我覺得是數學啊...
不過我如果有問題，接下來會去程式專板，謝謝

推 sin55688    : 補零後輸入怎會一樣呢?  完全看不懂
是不一樣，所以你其實看懂了
前面有人說一樣，我才看不懂

訊號輸入 => 主程式取樣 => FFT 副程式 => FFT 結果


好，這問題其實是這樣的
訊號輸入是 sin(x) 波
經過取樣，送入 FFT 副程式，產出 FFT 的結果

其中 FFT 副程式是個黑盒子，廠商寫的，我們不檢討它
所謂的一樣，所謂的不一樣，這是溝通字眼的問題
因為我已經附圖附程式了，所以我認為我的字眼不該有誤會

※ 編輯: HuangJC (223.141.233.39 臺灣), 11/20/2019 18:31:59