看板 Math 關於我們 聯絡資訊
: 好,來談這個問題 : 首先我雖然是理科的學生,但畢業很久,都還老師了 : 而且我沒學到的科目,用的可能不是精確的字眼 : 對我來說,誤差或是頻域取樣的問題 : 就是同一個問題 :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的值就對了 不懂人家給你的 11/20 11:56
recorriendo : 是什麼就指說有誤差是你本來自己的問題 11/20 11:57
recorriendo : 給你一公斤的東西你說台斤秤上顯示1.666斤這叫誤差? 11/20 11:58
我聽得懂你的意思,所以那我們就別用'誤差'這個詞 但還是得有一個詞來溝通吧? 請問那個詞是什麼? 謝謝
recorriendo : 還有頻率洩漏跟你用幾點做FFT本來就無關 要減少洩漏 11/20 11:59
既然不是頻率洩漏,那也就不用減少洩漏啊!
recorriendo : 請去搜尋window function 11/20 12:00
recorriendo : 專有名詞不是你愛怎樣定義就隨你定義 11/20 12:01
j0958322080 : 反正他高興就好了,他應該認為自己懂了 11/20 12:12
話不是這樣說,你們也說我錯,我同事也說我錯 我就像乒乓球一樣被踢來踢去 要是說,有才華的人總是有脾氣 我忍受完脾氣就有收獲,那也就值了 千萬不要什麼都得不到 那我當然認為自己已經懂了,止血 這難道不是對我最有利的做法? 上面大哥丟出一句'window function',那就是啦 有點收獲我就很感激了 (還有後面一篇回文;我再想想怎麼去吸收)
sunev : 這個板主要是討論數學,程式方面的限制請洽專板 11/20 12:57
戰文時才會拿這個當擋箭牌,如果專版能討論得更深入,那我會去專板 現在講究 PI 型人才(跨領域,站兩隻腳的意思) 數學當然為軟體服務 板上如果認為我洗板,那我可以回文時附加在自己文章之後,不增加文章 如果認為這不算數學,我覺得是數學啊... 不過我如果有問題,接下來會去程式專板,謝謝
sin55688 : 補零後輸入怎會一樣呢? 完全看不懂 11/20 14:36
是不一樣,所以你其實看懂了 前面有人說一樣,我才看不懂 訊號輸入 => 主程式取樣 => FFT 副程式 => FFT 結果 好,這問題其實是這樣的 訊號輸入是 sin(x) 波 經過取樣,送入 FFT 副程式,產出 FFT 的結果 其中 FFT 副程式是個黑盒子,廠商寫的,我們不檢討它 所謂的一樣,所謂的不一樣,這是溝通字眼的問題 因為我已經附圖附程式了,所以我認為我的字眼不該有誤會 ※ 編輯: HuangJC (223.141.233.39 臺灣), 11/20/2019 18:31:59