→ HuangJC: 上一個會在標題打我 ID的，我記得是 D姐.. 10/28 05:50
→ HuangJC: 兩個果然半斤八兩 10/28 05:50
推 joned: 上一個會叫我去別的板不要在這個是D姐 兩人要在一起? 10/28 07:38
→ st3336: 傅立葉玩出心得了嗎? 10/28 10:19
512 筆的傅立葉陣列完成了
我當初的問題是：輸入 512 筆，輸出也是同一個陣列 512 筆
可是輸出怎麼看？
有人回一句：時域轉頻域
這句我懂，如果我找一位老師也只教我這個，那不夠細節
可是學問本來就要自己做，所以我努力很久終於搞懂了
https://www.ptt.cc/bbs/Math/M.1571655907.A.7B7.html
這邊，發表在數學板，有程式有結果
然後湊答案發現，s[0] 要除以 4，才變成直流成份 1
而 s[1] 要除以 2，才變成振幅 1
為什麼這麼奇怪？但能回答我的人非常少
因為他們都只注重頻率響應，而不需要關注直流成份
更甚且，除以 2 也很少人在做，因為他們要看的常是'特徵頻率在哪裡'
而不是振幅有多大
只有我在用的東西當然只有我研究
https://www.itread01.com/content/1544320336.html
不過這個網址裡有和我相同的結論
他也沒有推導，而是直接給結論，不知是否和我一樣湊答案
反正就執行看看，然後看著結論湊一下就懂了
我接下來的挑戰是要做五萬筆陣列的傅立葉
同事馬上叫我不要做，他說沒人有這個需求，我一定是搞錯需求了
但他後來還是給了我方向：
你的應用只用到FFT.
如果資料太多可以切段作轉換 再將結果加權後相加。
FFT的長度看你需要的解析度決定
因為 你只需要顯示頻譜的大小而已，可以用上述方法進行。
以上教完了。
他搞錯了，誰說我只要大小...
加權要加多少權，他也沒說（應該是他不知道）
這邊我還是得自己研究
但有人給我一個網頁
http://www.fftw.org/
說是最權威的傅立葉函式庫
我同事給我的傅立葉明明很簡短，翻頁兩個 page down 就沒了
怎麼有人大費週章做函式庫？！
也許這種專業的有本事處理五萬筆數據
行的話當然用，反正又不用錢
（軟體工程師有很強烈的分享性格啊，一堆不用錢又強大的東西）
→ DaimlerHuang: 這裡已經有點像個板了.... 10/28 13:54
→ DaimlerHuang: 題外話,你用的傅立葉是快速複利葉,資料點數量 10/28 13:54
→ DaimlerHuang: 必須是2的次方,近6萬筆數字,使用65536輸入 10/28 13:55
→ DaimlerHuang: 輸出實部虛部各65536筆,至於直流部分,我不懂 10/28 13:57
→ DaimlerHuang: 理論上有頻率就不算是直流,直流大概是a[0] 10/28 13:58
數學和程式的差別，在於程式必需再考慮速度和溢位
65536筆很有可能是過大的陣列
如果引起溢位或堆疊不足，那我這程式也不用跑了
直流是 a[0] 沒錯，我程式也付上了，跑了就明白
→ st3336: 其實應該都要除4 因N=4 疊加四倍 1與3共軛 振福2疊加/4=1 10/28 14:27
→ st3336: 你的取樣點就是512 實際上1~256 對稱 257-512 10/28 14:43
→ st3336: 光看FFTW的feature基本上就很猛了 快速轉換/多工是重點 10/28 14:57
你是說，對稱的部份先疊加在一起，然後除一樣多的數？
因為對稱的關係，我是有一半的結果不看的
※ 編輯: HuangJC (111.82.25.183 臺灣), 10/28/2019 15:08:12
→ st3336: 你看s[1]跟s[3]是共軛的 實際上能量屬於同一個頻率 10/28 15:08
→ st3336: 總之在公式推導的結論上 s[1]+s[3] 代表真正的震幅 10/28 15:09
→ st3336: 而且轉換是一個疊加過程 你有n個點就會sigma n次 10/28 15:10
→ st3336: 所以最後要除n 10/28 15:10
→ st3336: 你math版有一個測試512的 你如果會畫圖 你畫1-512 就知道 10/28 15:11
→ st3336: 更精確地說 s[0]是第一點 2-512震幅譜會對稱 10/28 15:12
而且是以 s[256] 為中央向兩側對稱，這個中央我還不太敢用
→ st3336: s[0]就是你所謂的DC level 10/28 15:13
不要講 256 點那麼多，我們來講 16 點 （16hz 取樣）
輸入 s[i] = 1 +
1 * sin(x + 30.0 / 180 * M_PI) +
2 * sin(2 * x + 15.0 / 180 * M_PI) +
3 * sin(3 * x) +
4 * sin(4 * x) +
5 * sin(5 * x) +
6 * sin(6 * x);
因為我會讓 x 走完一個週期, 所以 sin(x) 就是 1hz, sin(6x)就是 6hz
我後來才搞懂，所謂幾 hz 那是我自己的解讀
因為程式裡根本沒含入這件事
但我把 sin(x)解讀為 1hz，這件事定義下來
也就會把輸出時的 s[1] 定義為 1hz 的係數項
上面 30/180*PI 這部份算式，都只是在基頻上轉個角度
這樣在觀察共軛上會更方便
來看輸出
after FFT
s[0]=(15.999999,0.000000) = 15.999999 deg 0.000000
s[1]=(4.000000,-6.928203) = 8.000000 deg -60.000003
s[2]=(4.141104,-15.454813) = 15.999999 deg -75.000001
s[3]=(-0.000000,-24.000000) = 24.000000 deg 89.999999
s[4]=(0.000000,-32.000000) = 32.000000 deg -89.999999
s[5]=(0.000001,-39.999999) = 39.999999 deg -89.999999
s[6]=(0.000001,-47.999999) = 47.999999 deg -89.999999
s[7]=(0.000001,-0.000000) = 0.000001 deg -6.930822
s[8]=(0.000001,0.000000) = 0.000001 deg 0.000000
s[9]=(0.000001,0.000000) = 0.000001 deg 6.930820
s[10]=(0.000001,47.999999) = 47.999999 deg 89.999999
s[11]=(0.000001,39.999999) = 39.999999 deg 89.999999
s[12]=(0.000000,32.000000) = 32.000000 deg 89.999999
s[13]=(-0.000000,24.000000) = 24.000000 deg -89.999999
s[14]=(4.141104,15.454813) = 15.999999 deg 75.000001
s[15]=(4.000000,6.928203) = 8.000000 deg 60.000003
s[0] 是 0hz
s[1] & s[15] 是 1hz
s[2] & s[14] 2hz
s[3] & s[13] 3hz
.
.
.
s[7] & s[9] 7hz
s[8]........ 最好不要用
這就我說不敢用中央項的原因
※ 編輯: HuangJC (111.82.25.183 臺灣), 10/28/2019 15:39:55
我一開始就開了取樣定理一個玩笑
當 s= sin(x); // 1hz
那我用 2hz 去取樣，2hz 有達到兩倍嘛，符合取樣定理
結果取到全是 0...
所以我想是要'高於兩倍'才安全，剛好兩倍用起來還是毛毛的
又如上面 16 項的輸出，中央項是 s[8], 按順序推算是 8hz 的係數
那我們也可以實驗一下含 8hz 的輸入
s[i] = 1 + 1 * sin(x) +
2 * sin(2 * x) +
3 * sin(3 * x) +
4 * sin(4 * x) +
5 * sin(5 * x) +
6 * sin(6 * x) +
7 * sin(7 * x) +
8 * sin(8 * x);
這次不含角度，不觀察共軛，只查看中央項的實用性
s[0]=(16.000000,0.000000) = 16.000000 deg 0.000000
s[1]=(0.000000,-7.999999) = 7.999999 deg -89.999999
s[2]=(0.000000,-15.999999) = 15.999999 deg -90.000000
s[3]=(-0.000000,-23.999999) = 23.999999 deg 89.999999
s[4]=(0.000000,-32.000000) = 32.000000 deg -90.000000
s[5]=(0.000000,-40.000000) = 40.000000 deg -90.000000
s[6]=(-0.000000,-48.000001) = 48.000001 deg 90.000000
s[7]=(0.000000,-56.000000) = 56.000000 deg -90.000000
s[8]=(-0.000000,0.000000) = 0.000000 deg 0.000000 <= 這項， 0
s[9]=(0.000000,56.000000) = 56.000000 deg 90.000000
s[10]=(-0.000000,48.000001) = 48.000001 deg -90.000000
s[11]=(0.000000,40.000000) = 40.000000 deg 90.000000
s[12]=(0.000000,32.000000) = 32.000000 deg 90.000000
s[13]=(-0.000000,23.999999) = 23.999999 deg -89.999999
s[14]=(0.000000,15.999999) = 15.999999 deg 90.000000
s[15]=(0.000000,7.999999) = 7.999999 deg 89.999999
8hz 的振幅應該是 1hz 的 8倍，這裡卻是 0
所以中央項我不敢用
※ 編輯: HuangJC (111.82.25.183 臺灣), 10/28/2019 15:57:47
→ st3336: 並不是中央向不要用 而是不要用到nyquist frequency 10/28 16:22
→ st3336: 如果你硬要他取樣16的話 8當然解不出來 取32就解決了 10/28 16:24
取樣 32hz 時，中央項已經變成 16 了
我不是在談 8hz 解不解得出來，我是在談 '中央項不要用'
所以你還是沒用到中央項啊，不是嘛...
而且別項都有共軛兩組，中央項沒有
真要用你把它當什麼？怎麼解讀呢？
會不會它就是疊加，兩組疊在一起
那既然是共軛，虛數項應該會消失，實數項變兩倍？
→ st3336: 以我處理時序資料為例子 不會遇到這麼低取樣率的 10/28 16:26
→ st3336: 主要還是端看你要處理的資料類型 10/28 16:27
喔～ 你有在處理時序資料？
那太好了，就是有在實戰的戰友
我的問題你應該都可以解了
目前我還在寫 UI，暫時把數學擱置
因為 UI 不好，做對做錯我也不知道..
要再等等才有下一次的回饋..
※ 編輯: HuangJC (111.82.25.183 臺灣), 10/28/2019 17:16:12
→ st3336: 我是沒自己寫拉 現成很多東西都寫好了 主要玩的還在後面 10/28 17:22
→ HuangJC: 我也無法再筆算積分，只想'叫用'已寫好的快速傅立葉 10/28 17:28
→ HuangJC: 但光是使用，這些討論就要懂啊，不然不會解讀傳回值 10/28 17:28
→ HuangJC: 請問你直流成份，各頻率係數是多少？我也只要這程度 10/28 17:28
→ st3336: 你的問題是甚麼 我用matlab帶你的case 結果是一樣的ˊ 10/28 18:07
→ st3336: 頂多差在截斷誤差而已 比如說7.9999 跟 8 的差異 10/28 18:07
我的問題解了啊，就是 s[0] 是什麼，s[1] 是什麼，s[2]是什麼...
之前只知道是頻域，現在知道是直流成份，1hz, 2hz, 3hz...的參數
而且也知道了必需除以 N 倍
---
下一個問題是我要寫出高通，低通數位濾波器
不過這得等我把 UI 完成了
→ st3336: 沒甚麼參數吧 也不一定是所謂 1 2 3HZ 10/28 20:32
→ st3336: 1 2 3端看你的點數/頻率取樣綠 10/28 20:54
→ st3336: 高通低通就是寫transfer function 看你要幾個pole這樣 10/29 00:16
→ st3336: 看你是要頻譜相乘 還是直接轉時域做convolution 10/29 00:17
直接來的，這題完全不要去想前面的傅立葉
因為我要做五萬筆的資料
五萬筆變成頻域，刪掉不要的頻率
再重新從頻域轉回時域
這種做法會瘋掉 XD
理論是理論，我要的是五萬筆刷～ for loop 跑完就做完高通，低通
程式我已經有了，拷來的
但跑起來不太對勁
所以等 UI 弄好再看問題
→ DaimlerHuang: 刪掉不要的頻率?五萬多筆數據怎麼能只宣告[512]??? 10/29 09:42
→ DaimlerHuang: 討論提到取樣,Nyquist卻違背取樣定理? 10/29 09:44
我看不懂你在打什麼，句子必需完整點
→ st3336: 實際上你把濾波器轉時域做convolution會比轉五萬比傅立葉 10/29 11:15
→ st3336: 來的快 10/29 11:15
是的
→ st3336: 你又要高通又要低通 幹嘛不直接帶通 10/29 11:16
因為沒拿到程式 XD
→ st3336: 我這邊不太懂你是五萬筆時間序列 還是一筆五萬點的時間序 10/29 11:18
我也看不懂你的句子
s[time] = v; //儀器所讀到的電壓，time 間隔單位是千分之一秒或百分之一秒
陣列從 0~5萬
這好像你講的兩個都是嘛..
推 st3336: 五萬筆 N點未知的時間序列 vs 一筆 五萬點的時間序列 10/29 13:01
→ st3336: 這裡這樣講 就是一筆 五秒的資料 取樣率1ms 這樣 10/29 13:02
如果是五萬個陣列，這我不會提出來
因為五萬個之間沒有交互關係
我自己 for loop 解決就好
→ st3336: 那就簡單多了 你也只是要做頻譜分析而已吧 10/29 13:20
我的問題對懂的人來說會非常簡單，不過一次談一件就好
DaimlerHuang八成就是混著看文章，所以回了我看不懂的句子
頻譜分析是頻譜分析
帶通是帶通
兩件我都要做，現在只談頻譜分析，而且我已經看懂怎麼做了
帶通等我 UI 完成，試過有想法再來改
→ st3336: 維基跟我想的差不多 帶通就是頻域高通*低通 10/29 13:32
https://imgur.com/a/JTorrBn
程式介面就是這樣
老闆不會管我用高通低通合併起來做帶通
反正做出來就好
這年頭，寫程式當然是先求有再求好 XD
如果帶通其實就是高通加低通
那萬一執行效率也一樣呢？
直接跑一個帶通花的時間和跑一次高通再跑一次低通花的時間一樣
如果是，我沒必要管它們是不是一起的..
所以搞懂弄出來就好
→ st3336: 我說維基百科的 你之前問得我在維基百科都有看到 10/29 13:41
→ HuangJC: 我說我有拿到程式，這些 google 不難 10/29 13:42
→ st3336: 我只是單純覺得帶通 應該就是頻率域的高通*低通 10/29 13:42
→ HuangJC: 但是，google 好幾個不同的地方，拿到的程式一樣 10/29 13:42
→ HuangJC: 內容農場；程式一樣，也一樣都不能跑，沒人解決 10/29 13:43
https://www.facebook.com/notes/%E9%BB%83%E7%91%9E%E6%98%8C/%E9%97%9C%E6%96%BC%E5%BB%BA%E7%AB%8B%E5%B7%A2%E7%8B%80-csplitterwnd/2696877440372540/
https://tinyurl.com/yymrh8jn
我舉另一例子好了
這是我寫的筆記
這例子我 google 到好幾個地方有，他們都錯成一樣的
照著寫程式都不會過
阿共那邊現在是不是很多內容農場啊？都抄來抄去的..
→ st3336: 頻率濾波器應該沒有太多分歧吧 我看都是transfer func 10/29 13:44
→ HuangJC: 等等我附上我研究到一半的東西你就知道我說哪裡沒解決了 10/29 13:45
→ HuangJC: 就像一個傅立葉，你說程式做的和 matlab 一樣 10/29 13:47
→ HuangJC: 一樣沒告訴我必需除以 N 啊，我必需除以 N 才能用 10/29 13:47
→ HuangJC: 而直流成份必需除以 2N，特別不一樣咧 10/29 13:48
→ HuangJC: 所以大部份談數位濾波器的，都是附一個多項式公式 10/29 13:48
→ HuangJC: 然後跑了 matlab 說：看，這式子有用，它就是一階濾波 10/29 13:49
→ st3336: 所以說 你的程式只做到一半 還沒有達到理論推導的結果 10/29 13:49
→ HuangJC: 但他們附的程式其實怪怪的 XDDD 10/29 13:49
→ HuangJC: 我手上的程式輸出會去修改輸入陣列，我看不懂這什麼邏輯 10/29 13:50
→ st3336: 陣列大小吧 要砍一半 10/29 13:51
→ HuangJC: 如果他想表達這叫 IIR, 程式設計技巧不該是這樣的 10/29 13:51
→ st3336: 你得另外開一個同大小的陣列 不然都用 s 原本input都丟了 10/29 13:51
→ HuangJC: 因這種改變輸入陣列的做法催毀原值，那我要另外保留原值? 10/29 13:52
→ HuangJC: 低通我做了有用，高通沒用，我準備留三套陣列了 10/29 13:52
→ HuangJC: 所以那程式有必要修正.. 10/29 13:52
→ st3336: 其實5萬點 應該也只是一秒內的事情 10/29 14:01
手上的程式有做到啊
其實我是競爭者，我不做的話目前已經有程式了
→ DaimlerHuang: ms間隔取樣5萬筆,所以收了50秒?1秒1000個資料 10/29 14:16
→ DaimlerHuang: 台灣交流電60Hz,一個週期有1000/60個資料 10/29 14:16
→ DaimlerHuang: 你需要對完整5萬個資料做FFT轉換? 還是 10/29 14:17
只需要對其中一個週期做FFT轉換?(其它週期視做相同
交流電應該是你舉的例，前面我舉 1+sin(x) 也是我為了了解傅利葉而做的測試
我實際的資料不是週期波
因此我認為要做五萬個資料做 FFT 轉換
我同事認為不可能，是因為他覺得還沒聽過這種需求
關於需求的部份我現在也不是很清楚，我只知道有一串數字是五萬多筆
這五萬筆是十秒內還是一秒內，我也還在問
不過這不影響傅立葉不是嘛
反正這五萬筆要餵進去
反正輸出陣列的頻率解析度，是根據輸入陣列的頻率來定義
現在或許可以賭：還真沒有五萬筆一次塞進去的需求
（這就是在猜老闆要什麼了）
其實或許真的沒需求吧！
因為塞進去愈多筆的差別是可以分析出愈低頻的資料
了不起就少了些超低頻的頻譜值
但我已經不想管老闆要不要了
身為工程師，我突然任性的想：為什麼我拿到的程式不能塞五萬筆？
而我同事則是說：別鬧了，快點賺錢，五萬筆沒人在要求的
→ st3336: 我是覺得用octave就好了 10/29 14:18
我是一個工程師，我的工作是產出一支程式能做
而不是教客戶：你去用別人的東西，把錢付給別人就好
就算 octave 能做，我也得挖出 octave 的程式包起來（它是 open source 吧）
變成執行檔交差
而不能教客戶說：來，裝套 octave,我教你，它很簡單
那我不成了 octave 的推廣業務？
不過如果你說的是：
必需要有開發工具，大家都推 matlab，我看你用 octave 當工具好了
我同意，有人說要教我
→ HuangJC: 所以我的工作其實就是讀檔，數學，繪圖 10/29 16:21
→ HuangJC: 這三個動作去操作 octave 是可以完成，但檔案有它的格式 10/29 16:22
→ HuangJC: 這個讀檔就是我程式的價值；本來就很簡單，已完成 10/29 16:22
→ HuangJC: 我是在補數學這一塊；對懂數學的人來說，它很簡單 10/29 16:22
→ st3336: 如果是要製作產品那就另當別論了 我一直以為你只是要研究 10/29 17:17
→ st3336: 其實也不用糾結五萬筆數值 一筆時間序列本來就是可變 10/29 17:33
→ st3336: 只要dimention開夠大 就不會出錯 10/29 17:35
→ st3336: 你的結果只要比對現有程式的結果 答案一樣就能用了 10/29 17:36
→ st3336: 你能用C寫GUI 我還蠻佩服的 10/29 17:38
UI 完成了階段任務，來做帶通吧～
上面有誤會澄清的東西我消一消，然後開新的一篇回文
※ 編輯: HuangJC (111.83.189.192 臺灣), 10/29/2019 18:04:40
→ HuangJC: 如果可以用工具我還談什麼五萬筆無法解傅立葉? 10/29 18:05
→ HuangJC: 所以我很前面就說了：數學和程式的差別，在於程式必需思 10/29 18:05
→ HuangJC: 考速度和溢位。所以現在，程式還是我自己面對，沒錯吧.. 10/29 18:06