Re: [問題] 還有人記得傅立葉變換嗎?
※ 引述《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
→
11/14 04:16,
4年前
, 1F
11/14 04:16, 1F
→
11/14 13:59,
4年前
, 2F
11/14 13:59, 2F
沒有,
https://zh.wikipedia.org/wiki/%E6%95%B0%E5%80%BC%E7%A8%B3%E5%AE%9A%E6%80%A7
是指這個嗎?
我是學過'用逼進法求解,有可能發散',不知是否指這回事
但還是不明白這裡為什麼要提這件事
→
11/15 17:33,
4年前
, 3F
11/15 17:33, 3F
好,來談這個問題
首先我雖然是理科的學生,但畢業很久,都還老師了
而且我沒學到的科目,用的可能不是精確的字眼
對我來說,誤差或是頻域取樣的問題
就是同一個問題 :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
→
11/16 03:13,
4年前
, 4F
11/16 03:13, 4F
→
11/16 03:17,
4年前
, 5F
11/16 03:17, 5F
→
11/16 03:18,
4年前
, 6F
11/16 03:18, 6F
→
11/16 16:10,
4年前
, 7F
11/16 16:10, 7F
→
11/16 16:12,
4年前
, 8F
11/16 16:12, 8F
→
11/16 16:13,
4年前
, 9F
11/16 16:13, 9F
→
11/16 16:13,
4年前
, 10F
11/16 16:13, 10F
→
11/16 16:17,
4年前
, 11F
11/16 16:17, 11F
→
11/16 16:18,
4年前
, 12F
11/16 16:18, 12F
→
11/18 15:24,
4年前
, 13F
11/18 15:24, 13F
→
11/18 15:25,
4年前
, 14F
11/18 15:25, 14F
→
11/18 15:25,
4年前
, 15F
11/18 15:25, 15F
→
11/18 15:26,
4年前
, 16F
11/18 15:26, 16F
→
11/18 15:27,
4年前
, 17F
11/18 15:27, 17F
> 反正你真的繼續取樣下去也是繼續拿到一串零。
才不是咧,我在想良好的表達都是要附圖的
最近我有想把圖畫出來講
(反之沒附圖而回應我的人,覺得我能聽懂?聽不懂的話溝通障礙要怪在我身上嗎?)
我不像熟 matlab 的人可以輕輕鬆鬆產出畫面
我都必需一行一行建制完整程式
有時我會想,我已經有自己的解決方案,可以停手了
但既然有人陪我討論,那這熱情值得鼓勵
我想法子生圖就是了
> 除此之外,zero padding當然是在惡搞輸入訊號,最後你會拿到一樣的結果才有鬼
怪的是這就是我要表達的啊,所以你應該懂我在說什麼...
那回過頭來說吧,竟然有人說錯的是我,因為程式是我寫的
Orz.. 而且還不只一個
我相信是溝通上出了問題,唉..
※ 編輯: HuangJC (223.141.233.39 臺灣), 11/20/2019 02:12:04
→
11/20 02:12,
4年前
, 18F
11/20 02:12, 18F
→
11/20 02:14,
4年前
, 19F
11/20 02:14, 19F
→
11/20 02:29,
4年前
, 20F
11/20 02:29, 20F
討論串 (同標題文章)